Description
This model demonstrates how to implement a generic simple outer approximation algorithm for convex mixed-binary non-linear programs. This model uses two represenations of the model instance: the index model instance (with original names) and a vector/scalar model instance. The latter allows to deal with the Jacobian and similar objects in a generic n by m matrix fashion. The modifications (OA cuts, solution elemination cuts, fixing of variables, etc) for the subproblems are done by GAMS source rewriting in the scalar model instance space. The only custom part that needs to be done by the user is to identify the non-linear constraints and the binary variables in the original indexed model. Contributor: Michael Bussieck, November 2012 START Problem specific part
Small Model of Type : GAMS
Category : GAMS EMP library
Main file : oa.gms
$title Outer Approximation for Convex Minimization Problem with Binary Variables (OA,SEQ=95)
$onText
This model demonstrates how to implement a generic simple outer approximation algorithm
for convex mixed-binary non-linear programs.
This model uses two represenations of the model instance: the index
model instance (with original names) and a vector/scalar model
instance. The latter allows to deal with the Jacobian and similar
objects in a generic n by m matrix fashion.
The modifications (OA cuts, solution elemination cuts, fixing of
variables, etc) for the subproblems are done by GAMS source rewriting
in the scalar model instance space.
The only custom part that needs to be done by the user is to identify
the non-linear constraints and the binary variables in the original
indexed model.
Contributor: Michael Bussieck, November 2012
$offText
* START Problem specific part
$if not set model $set model fuel
$call gamslib -q %model%
$goTo %model%
$label fuel
$set binvar status_vm
$set nlequ costfn_em, oileq_em
$set binvarc binvar(j) = sum(status_vm(j,u1),1);
$set nlequc nlequ(i) = 1$costfn_em(i) + sum(oileq_em(i,u1),1);
$set solve ucom min cost
$goTo cont
$label procsel
$set binvar y1_vm, y2_vm, y3_vm
$set nlequ inout2_em, inout3_em
$set binvarc binvar(j) = 1$y1_vm(j) + 1$y2_vm(j) + 1$y3_vm(j);
$set nlequc nlequ(i) = 1$inout2_em(i) + 1$inout3_em(i);
$set solve process min pr
$goTo cont
$label batchdes
$set binvar y_vm
$set nlequ time_em, obj_em
$set binvarc binvar(j) = sum(y_vm(j,u1,u2), 1);
$set nlequc nlequ(i) = 1$time_em(i) + 1$obj_em(i);
$set solve batch min cost
$goTo cont
* STOP Problem specific part
* Generic OA algorithm
$label cont
$onEcho > convert.op2
dumpgdx jacobian.gdx
objvar
$offEcho
$onEcho > convert.opt
gams
gmsinsert
objvar
dictmap
$offEcho
* Assume last solve statement is the MINLP to be solved
$call gams %model% minlp=convert optfile=1 lo=%gams.lo%
$if errorlevel 1 $abort problems running %model%
Set i scalar equations
j scalar variables
nlequ(i) scalar non-linear equations
binvar(j) scalar binary variables;
* Original binary variables and non-linear equation maps
Set %binvar%, %nlequ%;
$gdxIn dictmap
$load i j %binvar% %nlequ%
Alias (*,u1,u2,u3);
%binvarc%;
%nlequc%;
Scalar
sstat solve status
mstat model status
obj objective function;
Variables
x(j) scalar variables;
Equations
e(i) scalar equations;
file fx /'insert.gms'/; fx.pw=4095;
put fx 'solve m using rminlp min objvar'
/ 'scalar mstat,sstat,obj;'
/ 'mstat=m.modelstat; sstat=m.solvestat; obj=m.objval;'
/ 'execute_unload "rminlp", mstat, sstat, obj;'
/ 'option minlp=convert; m.optfile=2; solve m using minlp min objvar'
/ '$stop';
putclose;
execute 'gams gams lo=%gams.lo% u1=insert';
abort$errorlevel 'problems solving rminlp';
execute_load 'rminlp', mstat, sstat, obj;
if (not (mstat=%modelStat.optimal% or mstat=%modelStat.locallyOptimal%),
abort 'rminlp not solved to (local) optimality';
);
* Check if integer feasible
execute_load 'jacobian', x;
if (sum(binvar(j), abs(x.l(j)-round(x.l(j))))<1e-6,
abort.noerror 'rminlp integer feasible';
);
Set iter OA iterations / i1*i10 /
oacutpmarg(iter,i) OA cut positive marginal indicator
oacut(iter,i) OA cut
mipcut(iter) MIP cut
mipcut1(iter,j) MIP cut binary variables with level 1 in cut;
Parameter
lb lower bound
ub upper bound
A(i,j) Jacobian
oacutcoef(iter,i,j) OA cut coefficients
oacutrhs(iter,i) OA cut rhs;
lb = obj;
ub = inf;
mipcut(iter) = no; mipcut1(iter,j) = no;
alias(iter,it);
$set fmt 18:15
loop(iter,
execute_load 'jacobian', A, e, x;
* Build outer approximation cut:
loop(nlequ(i)$e.m(i),
oacutcoef(iter,i,j)$A(i,j) = A(i,j);
oacutrhs(iter,i) = sum(j$A(i,j), A(i,j)*x.l(j)) - e.l(i);
oacut(iter,i) = yes;
oacutpmarg(iter,i) = e.m(i)>0;
);
* Solve MIP
put fx;
loop(oacut(it,i),
put / 'equation oa_', it.tl:0 '_' i.tl:0 '; '
'variable boa_', it.tl:0 '_' i.tl:0 ' / lo ' e.lo(i) ', up ' e.up(i) ' /;'
/ 'oa_' it.tl:0 '_' i.tl:0 '.. ';
loop(j$oacutcoef(it,i,j), put$(oacutcoef(it,i,j)>0) '+';
put oacutcoef(it,i,j):%fmt% '*' j.tl:0);
if(oacutpmarg(it,i), put ' =g= '; else put ' =l= ');
put oacutrhs(it,i):%fmt% ' + boa_', it.tl:0 '_' i.tl:0 ';';
);
loop(mipcut(it),
put / 'equation mc_', it.tl:0 ';'
/ 'mc_' it.tl:0 '.. ';
loop(binvar(j), if (mipcut1(it,j), put '+' j.tl:0; else put '-' j.tl:0));
put ' =l= ' (sum(mipcut1(it,binvar),1)-1):0:0 ';';
);
put 'model mp / all';
loop(nlequ(i), put '-' i.tl:0); put '/;'
/ 'solve mp using mip min objvar;'
/ 'scalar mstat,sstat,obj;'
/ 'mstat=mp.modelstat; sstat=mp.solvestat; obj=mp.objval;'
/ 'Parameter binX(*);'
loop(binvar(j), put / 'binX("' j.tl:0 '") = ' j.tl:0 '.l;');
put / 'execute_unload "mip", mstat, sstat, obj, binX;'
/ '$stop';
putclose;
execute 'gams gams optcr=0 lo=%gams.lo% u1=insert';
abort$errorlevel 'problems solving mip';
execute_load 'mip', mstat, sstat, obj, x.l=binX;
if (not (mstat=%modelStat.optimal% or mstat=%modelStat.integerSolution%),
if (mstat=%modelStat.infeasible% or
mstat=%modelStat.integerInfeasible% or
mstat=%modelStat.infeasibleNoSolution%,
abort.noerror 'MIP infeasible', lb, ub;
else
abort 'mip not solved to optimality or infeasibility');
);
* lower bound update
lb = obj;
abort.noerror$(lb>ub) 'crossover', lb, ub;
* Add MIP cut
mipcut(iter) = yes;
mipcut1(iter,binvar(j)) = x.l(j)>0.5;
* Solve fixed NLP
loop(binvar(j), put fx / j.tl:0 '.fx = ' round(x.l(j)):0:0 ';');
put / 'solve m using rminlp min objvar'
/ 'scalar mstat,sstat,obj;'
/ 'mstat=m.modelstat; sstat=m.solvestat; obj=m.objval;'
/ 'execute_unload "rminlp", mstat, sstat, obj;'
/ 'option minlp=convert; m.optfile=2; solve m using minlp min objvar'
/ '$stop';
putclose;
execute 'gams gams lo=%gams.lo% u1=insert';
abort$errorlevel 'problems solving rminlp';
execute_load 'rminlp', mstat, sstat, obj;
if (mstat=%modelStat.optimal% or mstat=%modelStat.locallyOptimal%,
ub = min(ub, obj);
);
* Convergence
abort.noerror$((ub - lb) < 1e-4) 'converged', lb, ub;
);
display 'out of iterations', lb, ub;