$title simple MPEC unique solution doubly-bounded matches (MPEC07,SEQ=618)

$onText
  extremely simple MPEC model with some useful features:
   a) it has a unique solution
   b) all matched vars are doubly-bounded

Contributor: Steve Dirkse, Sep 2013
$offText

$if not set TESTTOL $set TESTTOL 1e-4

scalars
  tol / %TESTTOL% /
  c1  / -4 /
  c2  / 8 /
  c3  / 8 /
  ;
free variables
  x1  / lo   -1,   up   1 /
  x2  / lo    1,   up   3 /
  x3  / lo -100,   up 100 /
   y  /            up   0 /
  ;
free variable z;
equation f1, f2, f3, o;

f1 .. 3*x1 +   x2         + y + c1 =N= 0;
f2 ..   x1 + 3*x2 +   x3  + y + c2 =G= 0;
f3 ..          x2 + 3*x3  + y + c3 =N= 0;
o  .. sqr(x1-5) + sqr(x2-2) + sqr(x3-2) + sqr(y) =E= z;

model m / f1.x1, f2.x2, f3.x3, o /;

solve m using mpec min z;
if {(m.solvestat = %solveStat.capabilityProblems%),
  abort$[m.modelstat <> %modelStat.noSolutionReturned%] 'Wrong status codes',
    m.solvestat, m.modelstat;
  abort.noerror 'no solution, no point in checking further';
else
  abort$[not((m.solvestat = %solveStat.normalCompletion%) and (m.modelstat = %modelStat.optimal% or
              m.modelstat = %modelStat.locallyOptimal% or m.modelstat = %modelStat.feasibleSolution%))] 'Wrong status codes',
    m.solvestat, m.modelstat;
};

scalars r1, r2, r3;
r1 = 3*x1.l +   x2.l          + y.l + c1;
r2 =   x1.l + 3*x2.l +   x3.l + y.l + c2;
r3 =            x2.l + 3*x3.l + y.l + c3;

file log / '' /;
putclose log
 ' ' /
 'F1 := ', r1:6:2, '  perp x1 ', x1.lo:8:2, ' <= ', x1.L:6:2 ' <= ', x1.up:6:2 /
 'F2 := ', r2:6:2, '  perp x2 ', x2.lo:8:2, ' <= ', x2.L:6:2 ' <= ', x2.up:6:2 /
 'F3 := ', r3:6:2, '  perp x3 ', x3.lo:8:2, ' <= ', x3.L:6:2 ' <= ', x3.up:6:2 /
 '  y = ', y.L:6:2 /
 ' ' / ;

abort$(abs(x1.l-1.0) > tol) 'var x1.l should be  1.0', x1.l;
abort$(abs(x2.l-1.0) > tol) 'var x2.l should be  1.0', x2.l;
abort$(abs(x3.l+2.5) > tol) 'var x3.l should be -2.5', x3.l;
abort$(abs( y.l+1.5) > tol) 'var  y.l should be -1.5',  y.l;
abort$(abs( z.l-39.5)> tol) 'var  z.l should be 39.5',  z.l;

abort$(abs(f1.l+1.5) > tol) 'equ f1.l should be -1.5', f1.l;
abort$(abs(f2.l-0.0) > tol) 'equ f2.l should be  0.0', f2.l;
abort$(abs(f3.l-0.0) > tol) 'equ f3.l should be  0.0', f3.l;
abort$(abs( o.l-0  ) > tol) 'equ  o.l should be  0'  ,  o.l;