Description
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.
Small Model of Type : MIP
Category : GAMS Model library
Main file : trnsindic.gms
$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;