$title Fixed Charge Transportation Problem with Indicator Constraints (TRNSINDIC,SEQ=412)

$onText
This problem finds a least cost shipping schedule that meets
requirements at markets and supplies at factories and has a fixed charge
if something flows on an network arc.

The model demonstrates how to formulate some implication constraints
"binaryVariable = 1 => equation eq holds" via indicator constraints.
Usually a bigM forumlation can be used to express this logic, but indicator
constraint can do this without having an appriori finite bigM. This model shows
a couple of variants using the bigM and indicator constraints.

Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
Princeton University Press, Princeton, New Jersey, 1963.

Keywords: mixed integer linear programming, transportation problem, indicator
          constraints, scheduling, GAMS language features
$offText

* Only run Cplex, Gurobi, and SCIP
$ifThenI %system.mip%==cplex
$  echo readFile %system.mip%.indic > cplex.opt
$elseIfI  %system.mip%==gurobi
$  echo readFile %system.mip%.indic > gurobi.opt
$elseIfI  %system.mip%==scip
$  echo gams/indicatorfile="%system.mip%.indic" > scip.opt
$else
$  log *** No indicator constraint available with %system.mip%
$  exit
$endIf

Set
   i 'canning plants' / seattle,  san-diego /
   j 'markets'        / new-york, chicago, topeka /;

Parameter
   a(i) 'capacity of plant i in cases'
        / seattle    350
          san-diego  600 /

   b(j) 'demand at market j in cases'
        / new-york   325
          chicago    300
          topeka     275 /;

Table d(i,j) 'distance in thousands of miles'
              new-york  chicago  topeka
   seattle         2.5      1.7     1.8
   san-diego       2.5      1.8     1.4;

Scalar f 'freight in dollars per case per thousand miles' / 90 /;

Parameter
   c(i,j)       'transport cost in thousands of dollars per case'
   fixcost(i,j) 'fixed cost in thousands of dollars';

c(i,j)       = f*d(i,j)/1000;
fixcost(i,j) = 10*d(i,j)/1000;

Scalar
   minshipping 'minimum shipping of cases' / 100 /
   bigM        'sufficiently large number';

bigM = smax(i, a(i));

Variable
   x(i,j)   'shipment quantities in cases'
   use(i,j) 'is 1 if arc is used in solution'
   z        'total transportation costs in thousands of dollars';

Positive Variable x;
Binary   Variable use;

Equation
   cost         'define objective function'
   supply(i)    'observe supply limit at plant i'
   demand(j)    'satisfy demand at market j'
   minship(i,j) 'ensure minimum shipping'
   maxship(i,j) 'ensure zero shipping if use variable is 0';

cost..         z =e= sum((i,j), c(i,j)*x(i,j) + fixcost(i,j)*use(i,j));

supply(i)..    sum(j, x(i,j)) =l= a(i);

demand(j)..    sum(i, x(i,j)) =g= b(j);

minship(i,j).. x(i,j) =g= minshipping*use(i,j);

maxship(i,j).. x(i,j) =l= bigM*use(i,j);

option limRow = 0, limCol = 0, optCr = 0;

Model bigMModel / all /;

Alias (*,m,sm);

Parameter rep(m,sm,*);
$macro dorep(m,sm) rep('m','sm','modelstat') = m.modelStat; rep('m','sm','solvestat') = m.solveStat; rep('m','sm','obj') = m.objVal

solve bigMModel using mip minimizing z;
dorep(bigMModel,1);

* Now let's build a model for the same problem using indicator constraints
Equation
   iminship(i,j) 'ensure minimum shipping using indicator constraints'
   imaxship(i,j) 'ensure zero shipping if use variable is 0 using indicator constraints';

iminship(i,j).. x(i,j) =g= minshipping;

imaxship(i,j).. x(i,j) =e= 0;

Model indicatorModel / cost, supply, demand, iminship, imaxship /;

File fslv 'Indicator Option file' / %system.mip%.indic /;

putClose fslv 'indic iminship(i,j)$use(i,j) 1' / 'indic imaxship(i,j)$use(i,j) 0';

indicatorModel.optFile = 1;
solve indicatorModel using mip minimizing z;
dorep(indicatorModel,1);

* Let's do the same using indicator option with labels
loop((i,j),
   put fslv 'indic ' iminship.tn(i,j) '$' use.tn(i,j) yes
          / 'indic ' imaxship.tn(i,j) '$' use.tn(i,j) no /;
);
putClose fslv;

solve indicatorModel using mip minimizing z;
dorep(indicatorModel,2);

* Now let's build a model for the same problem that can be used with
* and without indicator constraints. This can become handy when
* debugging a model with indicator constraints

Positive Variable minslack(i,j), maxslack(i,j);

Equation
   xminship(i,j)    'ensure minimum shipping using indicator constraints and bigM'
   xmaxship(i,j)    'ensure zero shipping ig use variable is 0 using indicator constraints and bigM'
   bndminslack(i,j) 'ensure minslack is zero if use variable is 1'
   bndmaxslack(i,j) 'ensure maxslack is zero if use variable is 0';

xminship(i,j)..    x(i,j) =g= minshipping - minslack(i,j);

xmaxship(i,j)..    x(i,j) =e= 0 + maxslack(i,j);

bndminslack(i,j).. minslack(i,j) =l= bigM*(1 - use(i,j));

bndmaxslack(i,j).. maxslack(i,j) =l= bigM*use(i,j);

Model indicatorbigMModel / cost, supply, demand, xminship, xmaxship, bndminslack, bndmaxslack /;

* Let's first solve this without use of indicators
indicatorbigMModel.optFile = 0;
solve indicatorbigMModel using mip minimizing z;
dorep(indicatorbigMModel,1);

* Now we will use indicators and therefore we don't need the slacks
putClose fslv 'indic xminship(i,j)$use(i,j) 1' / 'indic xmaxship(i,j)$use(i,j) 0';

minslack.fx(i,j) = 0;
maxslack.fx(i,j) = 0;
indicatorbigMModel.optFile = 1;
solve indicatorbigMModel using mip minimizing z;
dorep(indicatorbigMModel,2);

* We can also mix and match bigM with indicator constraints
putClose fslv 'indic xminship(i,j)$use(i,j) 1';

minslack.fx(i,j) = 0;
maxslack.up(i,j) = inf;
indicatorbigMModel.optFile = 1;
solve indicatorbigMModel using mip minimizing z;
dorep(indicatorbigMModel,3);

$label done
Set mm(m,sm);
option mm < rep;

abort$(smax(mm,rep(mm,'modelstat')) > 1) 'Some model is not solved to global optimality',  rep;
abort$(smax(mm,rep(mm,'solvestat')) > 1) 'Some model has solvestat not Normal Completion', rep;

Scalar diff;
diff = abs(smax(mm,rep(mm,'obj')) - smin(mm,rep(mm,'obj')));
abort$(abs(smax(mm,rep(mm,'obj')) - smin(mm,rep(mm,'obj'))) > 1e-3) 'We get different objective values', rep;