$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 > convertd.op2 jacobian objvar $offecho $onecho > convertd.opt gams gmsinsert objvar dictmap $offecho * Assume last solve statement is the MINLP to be solved $call gams %model% minlp=convertd 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=convertd; 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.Locally Optimal%), 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.Integer Solution%), if (mstat=%ModelStat.Infeasible% or mstat=%ModelStat.Integer Infeasible% or mstat=%ModelStat.Infeasible No Solution%, 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=convertd; 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.Locally Optimal%, ub = min(ub, obj); ); * Convergence abort.noerror$((ub - lb) < 1e-4) 'converged', lb, ub; ); display 'out of iterations', lb, ub;