sipres.gms : Global optimization of semi-infinite programs via restriction of the right-hand side

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 /;);