$title Test CNS model with side constraints (cns14, SEQ=975)

$onText
  The square nonlinear system
    y = sqr(x)
    y = cos(x)
  has two solutions:
     x = +/-  0.82413
     y = sqr(x)

  In this test, we use a side constraint (not a variable bound but an actual
  constraint) to exclude one of the two solutions and find the other.
  We add two little wrinkles that don't change the basic idea:
    1. cos(x) is rejected by some global solvers, so we approximate
    2. we add a little bit of distortion so the NZ pattern is less trivial

  Contributor: Steven Dirkse, December 2024
$offText

$if not set TESTTOL $set TESTTOL 1e-6
scalar tol / %TESTTOL% /;
$if not set SLOWOK $set SLOWOK 0
scalar slowOK 'slow solves are OK: just abort.noerror in this case' / %SLOWOK% /;

* for small x this is close to cos(x)
$macro CLS(x) [1 - sqr(x)/2 + power(x,4)/24 - power(x,6)/720]

set j / 1 * 10 /;

parameter
  w(j) 'add a little wobble'
  xhat(j) 'positive x solution' /
    1  0.822198565229547
    2  0.81937320292214
    3  0.816547848660956
    4  0.813722946671731
    5  0.810898928514749
    6  0.808076212972153
    7  0.80525518578422
    8  0.802434859791175
    9  0.799524644509739
    10 0.79048366658928
  /
  yhat(j) 'y solution' /
    1  0.676010396454883
    2  0.671372397049558
    3  0.666750363401474
    4  0.66214502182321
    5  0.657557067471169
    6  0.652987164533831
    7  0.64843591397624
    8  0.643901704196337
    9  0.639239657178362
    10 0.624864424433268
  /
  ;
w(j) = 1.0 + 0.1 * ( [(ord(j)-1) / (card(j)-1)] - 0.5);

variables x(j), y(j);

equations
  sq(j)
  co(j)
  pl(j) 'positive x selected by linear constraint'
  nl(j) 'negative x selected by linear constraint'
  pn(j) 'positive x selected by nonlinear constraint'
  nn(j) 'negative x selected by nonlinear constraint'
  ;

sq(j)..  sqr(x(j)) =E= y(j);
co(j)..  CLS(x(j)*w(j)) - 0.05 * CLS(x(j+1)) =E= y(j);
pl(j)..  x(j) + y(j) =G= 1;
nl(j)..  x(j) =L= y(j) - 1;
pn(j)..  1 - sqr(x(j)-1) =G= y(j);
nn(j)..  y(j) =L= 1 - sqr(x(j)+1);

$ifThen set MCPSOL
model m / sq, co /;
x.l(j) = 0.8;
solve m using mcp;
parameter xhat(j), yhat(j);
xhat(j) = x.l(j);
yhat(j) = y.l(j);
execute_unload 'cnsSol.gdx', w, xhat, yhat;
$exit
$endif


* different equation orderings may test different permutations in the solver
model cnsplA / sq, co, pl /;
model cnsplB / sq, pl, co /;
model cnsplC / pl, co, sq /;

model cnsnlA / sq, co, nl /;
model cnsnlB / sq, nl, co /;
model cnsnlC / nl, co, sq /;

model cnspnA / sq, co, pn /;
model cnspnB / sq, pn, co /;
model cnspnC / pn, co, sq /;

model cnsnnA / sq, co, nn /;
model cnsnnB / sq, nn, co /;
model cnsnnC / nn, co, sq /;


$macro CHECKMODEL(m) abort$[(m.modelstat <> %modelStat.optimal%) and (m.modelstat <> %modelStat.solvedUnique%) and (m.modelstat <> %modelStat.solved%)] 'unexpected .modelstat', m.modelstat
$macro CHECKSOLVE(m) abort$[m.solvestat <> %solveStat.normalCompletion%] 'unexpected .solvestat', m.solvestat
$macro CHECKINFES(m) abort$[m.numinfes <> 0] 'unexpected .numinfes', m.numinfes 
$macro CHECKDEPND(m) abort$[m.numdepnd <> 0] 'unexpected .numdepnd', m.numdepnd 
$macro CHECKSLOW(m)  abort.noerror$[slowOK and %solvestat.ResourceInterrupt% = m.solvestat] 'Solve too slow';
$macro CHECKSTAT(m)  CHECKSLOW(m);  CHECKMODEL(m);  CHECKSOLVE(m);  CHECKINFES(m);  CHECKDEPND(m)

$macro CHECKXP  abort$[smax{j, abs(x.l(j)-xhat(j))} > tol] 'x.l not within tolerance', x.l, xhat, tol
$macro CHECKXN  abort$[smax{j, abs(x.l(j)+xhat(j))} > tol] 'x.l not within tolerance', x.l, xhat, tol
$macro CHECKY   abort$[smax{j, abs(y.l(j)-yhat(j))} > tol] 'y.l not within tolerance', y.l, yhat, tol
$macro CHECKCO  abort$[smax{j, abs(co.l(j)-co.up(j))} > tol]  'co.l not within tolerance', co.l, co.up, tol
$macro CHECKSQ  abort$[smax{j, abs(sq.l(j))} > tol]  'sq.l not within tolerance', sq.l, tol
$macro CHECKPL  abort$[smax{j, pl.lo(j)-pl.l(j)} > 0]  'equ pl not satisfied', pl.lo, pl.l
$macro CHECKNL  abort$[smax{j, nl.l(j)-nl.up(j)} > 0]  'equ nl not satisfied', nl.l, nl.up
$macro CHECKPN  abort$[smax{j, pn.lo(j)-pn.l(j)} > 0]  'equ pn not satisfied', pn.lo, pn.l
$macro CHECKNN  abort$[smax{j, nn.l(j)-nn.up(j)} > 0]  'equ nn not satisfied', nn.l, nn.up

$macro CHECKALL_P(m) CHECKSTAT(m);  CHECKXP;  CHECKY;  CHECKCO;  CHECKSQ
$macro CHECKALL_N(m) CHECKSTAT(m);  CHECKXN;  CHECKY;  CHECKCO;  CHECKSQ

* marginals are left untouched when the solution is returned: we check at model end
x.m(j) = NA;
y.m(j) = 7.0;
pl.m(j) = NA;
pn.m(j) = 0.5;
co.m(j) = 8.0;


x.l(j) = 0.7;  y.l(j) = 0.0;
solve cnsplA using cns;
CHECKALL_P(cnsplA);  CHECKPL;

x.l(j) = 0.7;  y.l(j) = 0.0;
solve cnsplB using cns;
CHECKALL_P(cnsplB);  CHECKPL;

x.l(j) = 0.7;  y.l(j) = 0.0;
solve cnsplC using cns;
CHECKALL_P(cnsplC);  CHECKPL;


x.l(j) = -0.7;  y.l(j) = 0.0;
solve cnsnlA using cns;
CHECKALL_N(cnsnlA);  CHECKNL;

x.l(j) = -0.7;  y.l(j) = 0.0;
solve cnsnlB using cns;
CHECKALL_N(cnsnlB);  CHECKNL;

x.l(j) = -0.7;  y.l(j) = 0.0;
solve cnsnlC using cns;
CHECKALL_N(cnsnlC);  CHECKNL;


x.l(j) = 0.7;  y.l(j) = 0.0;
solve cnspnA using cns;
CHECKALL_P(cnspnA);  CHECKPN;

x.l(j) = 0.7;  y.l(j) = 0.0;
solve cnspnB using cns;
CHECKALL_P(cnspnB);  CHECKPN;

x.l(j) = 0.7;  y.l(j) = 0.0;
solve cnspnC using cns;
CHECKALL_P(cnspnC);  CHECKPN;


x.l(j) = -0.7;  y.l(j) = 0.0;
solve cnsnnA using cns;
CHECKALL_N(cnsnnA);  CHECKNN;

x.l(j) = -0.7;  y.l(j) = 0.0;
solve cnsnnB using cns;
CHECKALL_N(cnsnnB);  CHECKNN;

x.l(j) = -0.7;  y.l(j) = 0.0;
solve cnsnnC using cns;
CHECKALL_N(cnsnnC);  CHECKNN;


abort$[smax{j, 1$(x.m(j) <> NA )} > 0] 'x.m was reset', x.m;
abort$[smax{j, 1$(y.m(j) <> 7.0)} > 0] 'y.m was reset', y.m;
abort$[smax{j, 1$(pl.m(j) <> NA )} > 0] 'pl.m was reset', pl.m;
abort$[smax{j, 1$(pn.m(j) <> 0.5)} > 0] 'pn.m was reset', pn.m;
abort$[smax{j, 1$(co.m(j) <> 8.0)} > 0] 'co.m was reset', co.m;