$title 'Test of SOS2 variables'  (SOS2A,SEQ=345)

$onText

Do linear interpolation of the points
   ord(I), f(I)
(using SOS2 vars w(I)) to define a function f(x).
By bounding w(I) below we can take out a little interval from the
domain of f.  We minimize the distance between f(x) and Fbar,
with full domain and with the restriction.

Contributor: Steve Dirkse
$offText

$if not set TESTTOL $set TESTTOL 1e-6
scalar mchecks / 0 /;
$if not %MIPMCHECKS% == 0 mchecks = 1;
scalar Fbar / 1.3 /;

set I / 1 * 3 /;

parameter f(I) /
1       1
2       2
3       3
/;

scalar wLo 'lower bound on w(1)' / -1 /;

sos2 variables w(I);
positive variables fplus, fminus;
* the optimization forces
* fplus  = min(0,fx-Fbar)
* fminus = min(0,Fbar-fx)

variables
  obj,
  x,
  fx;

w.lo(I) = 0;

equations
  wsum,
  xdef,
  fxdef,
  objdef,
  defwLo,
  gapplus,
  gapminus;

wsum..  1  =e= sum {I, w(I)};
xdef..  x  =e= sum {I, w(I)*ord(I)};
fxdef..        fx =e= sum {I, w(I)*f(I)};
gapplus..  fplus  =g= fx - Fbar;
gapminus.. fminus =g= Fbar - fx;
defwLo..   w('1') =g= wLo;
objdef..   obj =e= fplus + fminus;

model m / all /;
m.optcr = 0;

scalars
  tol / %TESTTOL% /,
  obj1    'objective of first solve'  / 0 /
  obj2    'objective of second solve' / .1 /
  fplus1    / 0 /
  fminus1   / 0 /
  fplus2    / 0 /
  fminus2   / .1 /
  fx1_L     / 1.3 /
  fx2_L     / 1.2 /
  x1_L      / 1.3 /
  x2_L      / 1.2 /
  ;
parameters
  w1_L(I) /
    1       .7
    2       .3
    3       0
  /
  w2_L(I) /
    1       .8
    2       .2
    3       0
  /
  ;
solve m using MIP minimizing obj;

if {(m.solvestat = %solveStat.capabilityProblems%),
  abort$(m.modelstat <> %modelStat.noSolutionReturned%)             'bad modelstat';
else
  abort$( m.solvestat <> %solveStat.normalCompletion% or m.modelstat <> %modelStat.optimal%) 'wrong status codes';
  abort$( abs(m.objval - obj1) > tol)           'Wrong m.objval';
  abort$( abs(obj.l - obj1) > tol)              'Wrong obj.l';
  abort$( smax(I,abs(w.l(I) - w1_L(I))) > tol)  'Wrong w.l';
  abort$( abs(fplus.l - fplus1) > tol)          'Wrong fplus.l';
  abort$( abs(fminus.l - fminus1) > tol)        'Wrong fminus.l';
  abort$( abs(fx.l - fx1_L) > tol)              'Wrong fx.l';
  abort$( abs(x.l - x1_L) > tol)                'Wrong x.l';
  if {mchecks,
    abort$( smax(I,abs(-wsum.m - ord(I) * xdef.m - f(I) * fxdef.m + defwLo.m$sameas(I,'1') + w.m(I))) > tol)  'Wrong w.m';
    abort$( abs(gapplus.m - objdef.m + fplus.m) > tol)           'Wrong fplus.m';
    abort$( abs(gapminus.m - objdef.m + fminus.m) > tol)         'Wrong fminus.m';
    abort$( abs(fxdef.m - gapplus.m + gapminus.m + fx.m) > tol)  'Wrong fx.m';
  };
};

wLo = .8;
solve m using MIP minimizing obj;
if {(m.solvestat = %solveStat.capabilityProblems%),
  abort$(m.modelstat <> %modelStat.noSolutionReturned%)             'bad modelstat';
else
  abort$( m.solvestat <> %solveStat.normalCompletion% or m.modelstat <> %modelStat.optimal%) 'wrong status codes';
  abort$( abs(m.objval - obj2) > tol)           'Wrong m.objval';
  abort$( abs(obj.l - obj2) > tol)              'Wrong obj.l';
  abort$( smax(I,abs(w.l(I) - w2_L(I))) > tol)  'Wrong w.l';
  abort$( abs(fplus.l - fplus2) > tol)          'Wrong fplus.l';
  abort$( abs(fminus.l - fminus2) > tol)        'Wrong fminus.l';
  abort$( abs(fx.l - fx2_L) > tol)              'Wrong fx.l';
  abort$( abs(x.l - x2_L) > tol)                'Wrong x.l';
  if {mchecks,
    abort$( smax(I,abs(-wsum.m - ord(I) * xdef.m - f(I) * fxdef.m + defwLo.m$sameas(I,'1') + w.m(I))) > tol)  'Wrong w.m';
    abort$( abs(gapplus.m - objdef.m + fplus.m) > tol)           'Wrong fplus.m';
    abort$( abs(gapminus.m - objdef.m + fminus.m) > tol)         'Wrong fminus.m';
    abort$( abs(fxdef.m - gapplus.m + gapminus.m + fx.m) > tol)  'Wrong fx.m';
  };
};