$title Test of =g= indicator constraints (INDIC01,SEQ=541)

$onText
Test =g= indicator constraints.  This model is based on DICE from the
GAMS model library.

Contributed by Steve Dirkse, November 2011
$offText

sets    f    faces on a dice  / face1*face6 /
        dice number of dice   / dice1*dice3 / ;

scalars flo  lowest face value  / 1 /
        fup  highest face value
        wn   wins needed - possible bound ;

fup = card(dice) * card(f);

wn = floor(0.5 * sqr(card(f))) + 1;

alias(f,fp);
alias(dice,dicep);

variables  wnx               'number of wins'
           fval(dice,f)      'face value on dice - may be fractional'
           comp(dice,f,fp)   'one implies f beats fp' ;
binary variable comp;

fval.lo(dice,f) = flo;
fval.up(dice,f) = fup;

fval.fx("dice1","face1") = flo;
* can cheat a little and fix values to speed things up
* fval.fx("dice1","face2") = 10;

equation eq1(dice)       'count the wins'
         eq3(dice,f,fp)  'definition of non-transitive relation'
         eq4(dice,f)     'different face values for a single dice';

eq1(dice).. sum((f,fp), comp(dice,f,fp)) =e= wnx;

eq3(dice,f,fp) .. fval(dice,f) =g= fval(dice++1,fp) + 1;

eq4(dice,f-1) .. fval(dice,f-1) + 1 =l= fval(dice,f);

model xdice /all /;

$onEcho > copt.o99
indic eq3(dice,f,fp)$comp(dice,f,fp) 1
$offEcho
$onEcho > cplex.o99
indic eq3(dice,f,fp)$comp(dice,f,fp) 1
$offEcho
$onEcho > gurobi.o99
indic eq3(dice,f,fp)$comp(dice,f,fp) 1
$offEcho
$onEcho > gurobiold.o99
indic eq3(dice,f,fp)$comp(dice,f,fp) 1
$offEcho
$onEcho > xpress.o99
indic eq3(dice,f,fp)$comp(dice,f,fp) 1
$offEcho
$onEcho > xpressold.o99
indic eq3(dice,f,fp)$comp(dice,f,fp) 1
$offEcho
$echo gams/indicatorfile = "cplex.o99" > scip.o99
xdice.optfile = 99;
* all obj values are integer, so any fairly small optcr should work
xdice.optcr = 1e-3;

* if you doubt that this formulation admits solutions with fractional fval,
* uncomment the variable fixing below
* fval.fx("dice3","face1") = 3.75;

solve xdice using mip max wnx;
abort$[xdice.solvestat <> %solveStat.normalCompletion%] "wrong solver status";
abort$[(xdice.modelstat <> %modelStat.optimal%) and (xdice.modelstat <> %modelStat.integerSolution%)] "wrong model status";

* clean and check
parameters
  chk1(dice)
  nWins(dice)    'win count, dice over dice++1'
  v(dice,f)      'face value of dice - will be integer'
  chk3(dice,f,fp)
  dom(dice,f,fp) 'one implies (dice,f) dominates (dice++1,fp)'
  ;

v(dice,f) = floor(fval.l(dice,f)+1e-3);
dom(dice,f,fp) = round(comp.l(dice,f,fp));
nWins(dice) = sum{(f,fp), dom(dice,f,fp)};

execute_unload 'indic01_sol', v, dom, nWins;

abort$[abs(wnx.l - round(wnx.l)) > 1e-6] 'objective not integer', wnx.l;

chk1(dice) = nWins(dice) - round(wnx.l);
abort$card(chk1) 'bogus objective', nWins, wnx.l, dom;

chk3(dice,f,fp) = abs(dom(dice,f,fp) - comp.l(dice,f,fp));
abort$[smax{(dice,f,fp), chk3(dice,f,fp)} > 1e-6] 'non-integer comp returned', chk3, comp.l;

chk3(dice,f,fp)$dom(dice,f,fp) = [v(dice,f) <= v(dice++1,fp)];
abort$card(chk3) 'bogus comp variable', comp.l, dom, chk3;