$title simple MPEC unique solution LB matches integer var (MPEC08,SEQ=619)

$onText
  extremely simple MPEC model with some useful features:
   a) it has a unique solution
   b) all matched vars are lower-bounded
   c) contains an integer var

Contributor: Steve Dirkse, Sep 2013
$offText

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

scalars
  tol / %TESTTOL% /
  c1  / -1.5 /
  c2  / -3 /
  ;
positive variables x1, x2;
integer variable y / lo 0,  up 10 /;
free variable z;
equation f1, f2, o;

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

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

scalars r1, r2;
file log / '' /;

* -----------------------------------------------------------------

solve m using rmpec 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;
};

r1 = 2*x1.l +   x2.l + y.l + c1;
r2 =   x1.l + 2*x2.l + y.l + c2;
putclose log
 ' ' /
 '** RMPEC: y = ', y.L:6:2 /
 '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 /
 ' ' / ;

abort$(abs(x1.l-0  ) > tol) 'var x1.l should be 0'  , x1.l;
abort$(abs(x2.l-0.8) > tol) 'var x2.l should be 0.8', x2.l;
abort$(abs( y.l-1.4) > tol) 'var  y.l should be 1.4',  y.l;
abort$(abs( z.l-5.8) > tol) 'var  z.l should be 5.8',  z.l;

abort$(abs(f1.l-0.7) > tol) 'equ f1.l should be 0.7', f1.l;
abort$(abs(f2.l-3.0) > tol) 'equ f2.l should be 3.0', f2.l;
abort$(abs( o.l-0  ) > tol) 'equ  o.l should be 0'  ,  o.l;

* -----------------------------------------------------------------

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.integerSolution%
              ] } ] 'Wrong status codes', m.solvestat, m.modelstat;
};

r1 = 2*x1.l +   x2.l + y.l + c1;
r2 =   x1.l + 2*x2.l + y.l + c2;
putclose log
 ' ' /
 '** MPEC: y = ', y.L:6:2 /
 '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 /
 ' ' / ;

abort$(abs(x1.l-0  ) > tol) 'var x1.l should be 0'  , x1.l;
abort$(abs(x2.l-1.0) > tol) 'var x2.l should be 1.0', x2.l;
abort$(abs( y.l-1.0) > tol) 'var  y.l should be 1.0',  y.l;
abort$(abs( z.l-6.0) > tol) 'var  z.l should be 6.0',  z.l;

abort$(abs(f1.l-0.5) > tol) 'equ f1.l should be 0.5', f1.l;
abort$(abs(f2.l-3.0) > tol) 'equ f2.l should be 3.0', f2.l;
abort$(abs( o.l-0  ) > tol) 'equ  o.l should be 0'  ,  o.l;