Description
This model provides an algorithm for the global solution of a
semi-infinite program (SIP) without convexity assumptions.
It terminates finitely with a guaranteed feasible point, and a
certificate of eps^f-optimality. The only assumptions are continuity
of the functions and existence of an eps^f-optimal SIP-Slater point.
The lower and upper bounds are obtained by the solution of two nonconvex
nonlinear programs (NLP) each, thus shifting the nonconvexity to the
global NLP solver.
This example provides a solution to the SIP
min{(x_1,x_2) in [-1000,1000]^2} 1/3*sqr(x_1) + sqr(x_2) + 1/2*x_1;
s.t.
sqr(1 - sqr(x_1)*sqr(y)) - x_1*sqr(y) - sqr(x_2) + x_2 <= 0 for all y in [0,1]
Model has been contributed by Alexander Mitsos, June 2009
Please use the citation below in derived work.
Large Model of Type : NLP
Category : GAMS Model library
Main file : sipres.gms
$title Global Optimization of semi-infinite Programs via Restriction of the right-hand side (SIPRES,SEQ=372)
$onText
This model provides an algorithm for the global solution of a
semi-infinite program (SIP) without convexity assumptions.
It terminates finitely with a guaranteed feasible point, and a
certificate of eps^f-optimality. The only assumptions are continuity
of the functions and existence of an eps^f-optimal SIP-Slater point.
The lower and upper bounds are obtained by the solution of two nonconvex
nonlinear programs (NLP) each, thus shifting the nonconvexity to the
global NLP solver.
This example provides a solution to the SIP
min{(x_1,x_2) in [-1000,1000]^2} 1/3*sqr(x_1) + sqr(x_2) + 1/2*x_1;
s.t.
sqr(1 - sqr(x_1)*sqr(y)) - x_1*sqr(y) - sqr(x_2) + x_2 <= 0 for all y in [0,1]
Model has been contributed by Alexander Mitsos, June 2009
Please use the citation below in derived work.
Mitsos, A, Global Optimization of Semi-Infinite Programs via
Restriction of the Right Hand Side. Optimization 60, 7 (2010).
http://doi.org/10.1080/02331934.2010.527970
Keywords: nonlinear programming, global optimization, semi-infinite programming,
mathematics
$offText
* Test only with deterministic global codes
$ifI %gams.nlp%==antigone $goTo cont
$ifI %gams.nlp%==baron $goTo cont
$ifI %gams.nlp%==lindoglobal $goTo cont
$ifI %gams.nlp%==scip $goTo cont
$log %gams.nlp% does not provide deterministic global optimization
$exit
$label cont
$onEmpty
Parameter
myoptR 'relative termination criteria for SIP' / 1e-3 /
myoptA 'absolute termination criteria for SIP' / 1e-3 /
epsilon 'initial value for epsilon' / 1 /
epsfactor 'division factor of epsilon' / 2 /;
Set
disc / 1*30 /
dlbd(disc) 'discretization of lower-bounding' / /
dubd(disc) 'discretization of upper-bounding' / /;
* dimensionality >= expected number of iterations
Scalar maxiter;
maxiter = card(disc);
Set xset, yset;
Parameter
ylbd(yset,disc) 'lower bound discretization'
yubd(yset,disc) 'upper bound discretizationl';
Variable
x(xset) 'variable for upper bounding and lower bounding problem'
y(yset) 'variables for lower-level'
obj 'upper-level objective function'
vio 'constraint violation for lower-level program';
Model lowmod, ubdmod, lbdmod;
* Problem spefic definitions
Equation
eqobj 'upper-level objective function'
eqconlow 'SIP constraint'
eqconlbd(disc) 'discretized constraints for lower bounding problem'
eqconubd(disc) 'discretized constraints for upper bounding problem';
Set
xset / 1*2 /
yset / 1*1 /;
$macro resetBoundsOnX x.lo(xset) = -1000; x.up(xset) = 1000
resetBoundsOnX;
y.lo(yset) = 0;
y.up(yset) = 1;
eqobj.. obj =e= 1/3*sqr[x('1')] + sqr[x('2')] + 1/2*x('1');
$macro eqcon(y) sqr{1 - sqr[x('1')]*sqr[y]} - x('1')*sqr[y] - sqr[x('2')] + x('2')
eqconlow.. eqcon(y('1')) =e= vio;
eqconlbd(dlbd).. eqcon(ylbd('1',dlbd)) =l= 0;
eqconubd(dubd).. eqcon(yubd('1',dubd)) =l= -epsilon;
* the three subproblems
Model
lowmod / eqconlow /
ubdmod / eqobj, eqconubd /
lbdmod / eqobj, eqconlbd /;
* end of problem specific declarations
File output where results go / 'ex2res.txt' /;
put output;
lowmod.holdFixed = 1;
option solPrint = off, limRow = 0, limCol = 0, optCa = 1e-6, optCr = 1e-4;
Scalar
LBD 'lower bound' / -1e10 /
UBD 'upper bound' / +1e10 /
UBD_LBD 'absolute gap'
iter 'iteration counter' / 0 /
cpulbd / 0 /
cpulbdlow / 0 /
cpuubd / 0 /
cpuubdlow / 0 /
cputot;
Parameter xstar(xset) 'solution';
* dummy initialization just in case
ylbd(yset,dlbd) = 6e66;
yubd(yset,dubd) = 6e66;
$eolCom //
$macro putsol(v) loop(v&set, put ' ' v.l(v&set):20:15)
$macro putstat(m) put 'm solveStat=' m.solveStat ' modelStat=' m.modelStat
loop(disc,
iter = iter + 1;
if(iter > maxiter,
put 'was not expecting iter=' iter /;
abort 'too many iterations';
);
put '####' / 'running at iter=' iter:5:0 ' LBD=' LBD:20:15 ' UBD=' UBD:20:15 /;
solve lbdmod minimizing obj using nlp; // lower bounding problem
cpulbd = cpulbd + lbdmod.resUsd;
putstat(lbdmod);
if(lbdmod.solveStat <> %solveStat.normalCompletion% or
(lbdmod.modelStat <> %modelStat.optimal% and
lbdmod.modelStat <> %modelStat.feasibleSolution% and
lbdmod.modelStat <> %modelStat.locallyOptimal%),
put ' unexpected failed lbdmod' /;
if(LBD < lbdmod.objEst, LBD = lbdmod.objEst;);
else
LBD = lbdmod.objEst;
);
put ' lbd obj=' obj.l:20:15 ' best possible=' lbdmod.objEst:20:15,' x='; putsol(x);
put /;
* lower-level problem for fixed parameters (assuming all good so far)
x.fx(xset) = x.l(xset);
solve lowmod maximizing vio using nlp;
resetBoundsOnX;
cpulbdlow = cpulbdlow + lowmod.resUsd;
putstat(lowmod);
if(lowmod.solveStat = %solveStat.normalCompletion% and
(lowmod.modelStat = %modelStat.optimal% or
lowmod.modelStat = %modelStat.feasibleSolution% or
lowmod.modelStat = %modelStat.locallyOptimal%),
put ' normal vio=' vio.l:20:15 ' best possible=' lowmod.objEst:20:15 ' y='; putsol(y);
put /;
if(lowmod.objEst <= 0,
put 'lower bounding problem found feasible point but will not use' /;
else
dlbd(disc) = yes;
ylbd(yset,disc) = y.l(yset);
);
);
* upper bounding problem
solve ubdmod minimizing obj using nlp;
cpuubd = cpuubd + ubdmod.resUsd;
putstat(ubdmod);
if(ubdmod.solveStat <> %solveStat.normalCompletion%,
put 'abnormal' /;
abort 'failed ubdmod';
elseIf(ubdmod.modelStat = %modelStat.infeasible% or
ubdmod.modelStat = %modelStat.infeasibleNoSolution%),
put ' infeasible epsilon=' epsilon:20:15 /;
epsilon = epsilon/epsfactor;
elseIf(ubdmod.modelStat <> %modelStat.optimal% and
ubdmod.modelStat <> %modelStat.feasibleSolution% and
ubdmod.modelStat <> %modelStat.locallyOptimal%),
put 'unexpected' /;
abort "unexpected ubdmod.modelstat";
else
put ' optimal candidate obj=' obj.l:20:15 ' best possible=' ubdmod.objEst:20:15 ' x='; putsol(x);
put ' epsilon=' epsilon:20:15 /;
* lower level problem for fixed parameters
x.fx(xset) = x.l(xset);
solve lowmod maximizing vio using nlp;
resetBoundsOnX;
cpuubdlow = cpuubdlow + lowmod.resUsd;
putstat(lowmod);
if(lowmod.solveStat <> %solveStat.normalCompletion% or
(lowmod.modelStat <> %modelStat.optimal% and
lowmod.modelStat <> %modelStat.feasibleSolution% and
lowmod.modelStat <> %modelStat.locallyOptimal%),
put 'unexpected' /;
abort 'failed lowmod';
else
put ' normal vio=' vio.l:20:15 ' best possible=' lowmod.objEst:20:15 ' y='; putsol(y);
put /;
if(lowmod.objEst <= 0,
epsilon = epsilon/epsfactor;
if(obj.l < UBD,
UBD = obj.l;
xstar(xset) = x.l(xset));
else
* add point to discretization
dubd(disc) = yes;
yubd(yset,disc) = y.l(yset));
);
);
UBD_LBD = UBD - LBD;
* check that lower bound is below upper bound (otherwise sth is wrong)
if(UBD_LBD < 0,
put 'UBD=' UBD:20:15 ' should be larger than LBD=' LBD:20:15 /;
abort 'unexpected UBD < LBD';
);
* check if termination has occured
if(UBD_LBD < myoptA or LBD > UBD - myoptR*abs(UBD),
put 'converged with UBD=' UBD:10:7 ' LBD=' LBD:10:7 ' difference=' UBD_LBD:10:7 ' xstar=';
loop(xset, put ' ' xstar(xset):20:15;);
put /;
break;
);
);
* printout discretization for upper and lower bounding problem
put / '##' / 'discretization for lbd' /;
loop(disc,
put dlbd(disc);
if(dlbd(disc), loop(yset, put ' ' ylbd(yset,disc):20:15););
put /;
);
put / '##' / 'discretization for ubd' /;
loop(disc,
put dubd(disc);
if(dubd(disc), loop(yset, put ' ' yubd(yset,disc):20:15););
put /;
);
* printout CPU time and solution
cputot = cpulbd + cpulbdlow + cpuubd + cpuubdlow;
put / '##' / 'Solution'
/ cputot:10:2
/ cpulbd:10:2
/ cpulbdlow:10:2
/ cpuubd:10:2
/ cpuubdlow:10:2
/ LBD:10:5
/ UBD:10:5 /;
loop(xset, put xstar(xset):20:15 /;);