$title 'MCP GUSS Test' (SCENSOL10,SEQ=920)

$onText
This model tests various combination of fixing variables in an MCP
including GUSS scenario runs.

The original model (which is also the base case) has more variables
than equation and when given directly to the MCP solver, an error
occurs. The fixing (and unfixing) of some variables make it a proper
MCP in the scenarios.

Contributor: Michael Bussieck, December 2022
$offText

Scalars
    x0 / 0 /
    y0 / 0 /;

Variables x, y, xd, yd;
Equations f, fx, fy;

f..   exp(x + y) =E= exp(2);
fx..  x =E= x0 + xd;
fy..  y =E= y0 + yd;

Model m / f, fx.x, fy.y /;

maxExecError = 1;
solve m using mcp;
execError = 0;
abort$(m.solveStat<>%solveStat.solveProcessingSkipped%) 'Expect base case MCP to fail';

* First scenario
xd.fx = 0;
solve m using mcp;
abort$(m.modelStat<>1 or m.solveStat<>1 or abs(x.l-0)>1e-6 or abs(y.l-2)>1e-6) 'xd.fx failed';

* Second scenario
xd.lo = -inf; xd.up = inf;
yd.fx = 0;
solve m using mcp;
abort$(m.modelStat<>1 or m.solveStat<>1 or abs(x.l-2)>1e-6 or abs(y.l-0)>1e-6) 'yd.fx failed';
yd.lo = -inf; yd.up = inf;

$macro check(case) loop(s, abort$(srep(s,'modelStat')<>1 or srep(s,'solveStat')<>1 or \
                                  abs(xval('s1')-0)>1e-6 or abs(xval('s2')-2)>1e-6 or \
                                  abs(yval('s1')-2)>1e-6 or abs(yval('s2')-0)>1e-6) case,srep,xval,yval)

Set
    s 'scenarios' / s1*s2 /
    mattrib       / system.GUSSModelAttributes /;
Parameters
    xVal(s), yVal(s)
    xFx(s)        / s1 EPS        /
    yFx(s)        /        s2 EPS /
    srep(s, mattrib)  'model attibutes like modelstat etc'
    o(*)              'GUSS options' / skipBaseCase 1, updateType 1 /;

Set dictFX /  s.scenario.''
              o.opt.srep
             xd.fixed.xFx
             yd.fixed.yFx
              x.level.xVal
              y.level.yVal /;
solve m using mcp scenario dictFX;
check('fx');

* Now with bounds
* BaseCase
Parameters
    xdLo(s) / s1  EPS           /
    xdUp(s) / s1  EPS           /
    ydLo(s) /           s2  EPS /
    ydUp(s) /           s2  EPS /;

Set dictLOUP /  s.scenario.''
                o.opt.srep
               xd.lower.xdLo
               xd.upper.xdUp
               yd.lower.ydLo
               yd.upper.ydUp
                x.level.xVal
                y.level.yVal /;

solve m using mcp scenario dictLOUP;
check('baseCase');

* Accumulate
xdLo('s2') = -INF; xdUp('s2') = INF;
o('updateType') = 2;
solve m using mcp scenario dictLOUP;
check('accumulate');

* Zero
ydLo('s1') = -INF; ydUp('s1') = INF;
o('updateType') = 0;
solve m using mcp scenario dictLOUP;
check('zero');