$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;