$title Non-transitive Dice Design (DICE,SEQ=176)

$onText
Probabilistic dice - an example of a non-transitive relation.
We want to design a set of dice with an integer number on each face
such that on average dice1 beats dice2, and dice2 on average beats
dice3 etc, but diceN has to beat dice1.

MIP codes behave very erratically on such a problem and slight
reformulations can result in dramatic changes in performance.
For example, making the variable fval integer can change performance
so we leave fval free and do a second solve to get integer values.

Gardner, M, Scientific American.

Robert A Bosch, Mindsharpener, Optima, MP Society Newsletter, Vol 70,
June 2003, page 8-9

Robert A Bosch, Monochromatic Squares, Optima, MP Society Newsletter,
Vol 71, March 2004, page 6-7

Keywords: mixed integer linear programming, dice designment, mathematics,
          nontransitive dice
$offText

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

Scalar
   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), (dice,dicep);

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

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) + (fup - flo + 1)*(1 - comp(dice,f,fp)) =g= fval(dice++1,fp) + 1;

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

Model xdice / all /;

$if set nosolve $exit

xdice.resLim = 30;

* 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%] "unexpected solver status";
abort$[(xdice.modelstat <> %modelStat.optimal%)
   and (xdice.modelstat <> %modelStat.integerSolution%)] "unexpected model status";

* the dominance relationship is good after this solve but there is no reason
* fval must be integer: we can fix that with a second solve of a related model

variable z 'sum of face values';
equation zdef;
zdef..   sum{(dice,f), fval(dice,f)} =e= z;
model facecomp / zdef, eq3, eq4 /;

comp.fx(dice,f,fp) = round(comp.l(dice,f,fp));
solve facecomp using rMIP min z;
abort$[facecomp.solvestat <> %solveStat.normalCompletion%] "unexpected solver status";
abort$[facecomp.modelstat <> %modelStat.optimal%]          "unexpected model status";

* rounding will give us integer face values
* we verify we get the expected win totals
parameters
  fv(dice,f)      'integer face values'
  win(dice,f,fp)  'tally wins using the face values'
  winTot(dice)    'computed win totals'
  ;
fv(dice,f) = round(fval.l(dice,f));
win(dice,f,fp) = [fv(dice,f) >= fv(dice++1,fp) + 1];
winTot(dice) = sum{(f,fp), win(dice,f,fp)};
abort$[smin{dice, round(winTot(dice) - wnx.L)} < 0] 'where did the wins go??', winTot, wnx.L;

option  fv:0;
display wn, wnx.L, fv;

Parameter rep1 'chance of winning against next';
rep1(dice,f) = 100*sum(fp, win(dice,f,fp))/card(f);
rep1(dice,'chance') = sum(f, rep1(dice,f))/card(f);

option  rep1:0;
display rep1;