$title Two dimensional scenario GUSS Example (GUSS2DIM,SEQ=423)

$onText
This model is based on the nonconvex mathopt1. Local solvers fail to
find the global optimum. Here we implement a search strategy and partition
the x,y space into equal sized boxes. For each box we start a local optimization
and this way get the global optimum.

This example demonstrates the use of GUSS with a two dimensional scenario
index. The set sample(s1,s2) represent all the two dimensional boxes described
by the lower (x/yLo) and upper (x/yUp) bounds of the boxes. The x,y solution
is stored in parameters sx,sy.

Bussieck, M R, Ferris, M C, and Lohmann, T, GUSS: Solving Collections of Data
Related Models within GAMS. In Kallrath, J, Ed, Algebraic Modeling Systems: Modeling
and Solving Real World Optimization Problems. Springer, Berlin Heidelberg, 2012, pp. 35-56.

Keywords: nonlinear programming, mathematics, global optimization, GUSS, multi-start
$offText

$eolCom //
Variable x, y, obj;

x.lo = -10; y.lo = -15;    // lower bounds
x.l  =   8; y.l  = -14;    // initial value
x.up =  20; y.up =  20;    // upper bounds

Equation objdef, eqs, ineqs;

objdef.. obj =e= 10*sqr(sqr(x) - y) + sqr(x - 1);

eqs..    x   =e= x*y;

ineqs..  3*x + 4*y =l= 25;

Models m / all /;

* NUMBOX is the number of boxes to cover the [x.lo,x.up]
$if not set NUMBOX $set NUMBOX 20
set s / s1*s%NUMBOX% /; alias (s,s1,s2);
set sample(s1,s2) / #s1.#s2 /;
Parameter xLo(s1,s2), xUp(s1,s2), yLo(s1,s2), yUp(s1,s2) 'Boxes that cover the [x.lo,x.up]*[y.lo,y.up] space';
xLo(s,s2) = x.lo + (x.up-x.lo)*(s.pos-1)/%NUMBOX%;
xUp(s,s2) = xLo(s,s2) + (x.up-x.lo)/%NUMBOX%;

yLo(s1,s) = y.lo + (y.up-y.lo)*(s.pos-1)/%NUMBOX%;
yUp(s1,s) = yLo(s1,s) + (y.up-y.lo)/%NUMBOX%;

Set mattrib 'model attributes like objval and modelstat' / system.GUSSModelAttributes /;

Parameter
   sx(s1,s2), sy(s1,s2)  'collectors for x and y'
   srep(s1,s2, mattrib)  'model attibutes like modelstat etc'
   o                     'GUSS options' / SkipBaseCase 1, LogOption 1 /;

Set dict / sample.scenario.''
           o.     opt     .srep
           x.     lower   .xLo
           x.     upper   .xUp
           x.     level   .sx
           y.     lower   .yLo
           y.     upper   .yUp
           y.     level   .sy   /;

solve m minimizing obj using nlp scenario dict;

* Only consider the samples that have been solved to at least feasibility (filter the infeasible ones)
sample(s1,s2) = srep(s1,s2,'ModelStat')=%modelStat.optimal% or
                srep(s1,s2,'ModelStat')=%modelStat.locallyOptimal% or
                srep(s1,s2,'ModelStat')=%modelStat.IntermediateNonoptimal%;

* Get the global minimum from all feasible samples and assign x.l and y.l with the corresponding sample level values
obj.l = smin(sample, srep(sample,'ObjVal'));

* Since there might be more than one sample with the global optimum we allow the singleton to be assigned multiple samples but take the first
singleton set bestSample(s1,s2); option strictSingleton=0; // See UG_GamsCall.html#GAMSAOstrictSingleton
bestSample(sample) = abs(srep(sample,'ObjVal')-obj.l)<1e-6;
x.l = sx(bestSample); y.l = sy(bestSample);
display bestSample, obj.l, x.l, y.l;

*Parameter sobj(s1,s2) 'collector for obj';
*sobj(sample)  = srep(sample,'ObjVal');
*option sobj:0:1:1; display sobj;

* Some test that we are close to the global minimum
abort$(abs(obj.l - 0)>1e-3) 'Not the global minimum', obj.l;
abort$(abs(x.l   - 1)>1e-3) 'Not the correct x',        x.l;
abort$(abs(y.l   - 1)>1e-3) 'Not the correct y',        y.l;