relief.gms : Relief Mission

**Description**

Huts or villages are located at 20 location on a 10 x 10 grid. The problem is to find the two best drop locations for relief packages such that the total distance traveled by villagers is minimized. We look at this problem as a plant location problem and trace the efficiency frontier of multiple drops. The optimal solution for two drops is D.8 and E.2. Another interesting question may be: how many other drop locations are there that have a total distance traveled with, say, 3 percent of the best location.

**Reference**

- Toczek, J, Relief Mission. ORMS Today 37, 4 (2010), 14.

**Large Model of Type :** MIP

**Category :** GAMS Model library

**Main file :** relief.gms

```
$title Relief Mission (RELIEF,SEQ=353)
$onText
Huts or villages are located at 20 location on a 10 x 10 grid. The problem
is to find the two best drop locations for relief packages such that the total
distance traveled by villagers is minimized.
We look at this problem as a plant location problem and trace the efficiency
frontier of multiple drops. The optimal solution for two drops is D.8 and E.2.
Another interesting question may be: how many other drop locations are there that
have a total distance traveled with, say, 3 percent of the best location.
Toczek J, Relief Mission, ORMS Today 37, 4 (2010), page 14
Keywords: mixed integer linear programming, plant location problem
$offText
$eolCom //
Set
r 'rows' / A*J /
c 'columns' / 1*10 /
d 'drops' / d1*d20 /;
Table dem(r,c) 'relief demands'
1 2 3 4 5 6 7 8 9 10
A 1
B 1 1 1
C 1 1 1 1 1
D 1 1 1
E 1 1
F 1
G 1
H 1 1
I
J 1 1;
Alias (r,rr), (c,cc);
Set huts(r,c) 'hut locations';
Parameter
dis(r,c,r,c) 'eucledian distances'
numdrops 'number of drops'
maxdem 'maximum drop demand';
huts(r,c) = dem(r,c);
maxdem = sum(huts, dem(huts));
dis(huts(r,c),rr,cc) = edist(r.pos-rr.pos,c.pos-cc.pos);
Variable
drop(r,c) 'drop locations'
walk(r,c,r,c) 'distances walked from huts to nearest drop zone'
total 'total distance walked';
Binary Variable drop;
Positive Variable walk;
Equation
demand
supply
deftotal
defnumdrop;
demand(huts).. sum((r,c), walk(huts,r,c)) =e= 1;
supply(r,c).. sum(huts, walk(huts,r,c)) =l= drop(r,c)*maxdem;
deftotal.. total =e= sum((huts,r,c), dis(huts,r,c)*walk(huts,r,c));
defnumdrop.. sum((r,c), drop(r,c)) =e= numdrops;
Model m / all /;
* --------------------------
* get best two drop solution
* --------------------------
numdrops = 2;
option limRow = 0, limCol = 0, optCr = 0, solveLink = 2, resLim = 10;
m.solPrint = 0;
solve m min total using mip;
display drop.l;
* ------------------------
* trace efficient frontier
* ------------------------
Parameter QDrep 'quick and dirty report';
loop(d,
numdrops = d.pos;
solve m min total using mip;
m.solPrint = 2;
// we use round() because a MIP code (CPLEX) may return fractional values
QDrep(r,c,d)$round(drop.l(r,c)) = sum ( huts, dis(huts,r,c)*walk.l(huts,r,c)) + eps*huts(r,c)*drop.l(r,c);
QDrep('max','dist',d) = smax((huts,r,c), dis(huts,r,c)*walk.l(huts,r,c));
QDrep('tot','dist',d) = sum ((huts,r,c), dis(huts,r,c)*walk.l(huts,r,c));
QDrep('CPU','used',d) = m.resUsd;
);
display QDrep;
* ----------------------------------------------------------------
* find all drop points within x percent of best for two drop points
*
* To exclude the n'th integer solution we can write:
* cut(n).. sum(i, abs(x(i) - xsol(i,n)) =g= 1;
* simulating abs() will give
* cut(n).. sum((r,c), cutval(n,r,c)*drop(r,c)) =l= 1;
* where cutval() contains solutions to be excluded
* ----------------------------------------------------------------
Set
nn 'max number of close solutions' / n1*n50 /
n(nn) 'dynamic set';
Parameter
limit 'total distance for drop points'
objval(nn) 'total miles traveled'
cutval(nn,r,c) 'all possible solutions for cut generation';
numdrops = 2;
solve m min total using mip;
limit = 1.03*total.l;
Equation cut(nn) 'known solutions to be eliminated';
cut(n).. sum((r,c), cutval(n,r,c)*drop(r,c)) =l= 1;
Model mm / m, cut /;
n(nn) = no; // clear set of cuts
mm.solveStat = %solveStat.NormalCompletion%;
mm.modelStat = %modelStat.Optimal%;
mm.solPrint = 0;
loop(nn$(mm.solveStat = %solveStat.NormalCompletion% and
mm.modelStat = %modelStat.Optimal% and
total.l < limit),
n(nn) = yes; // add element to set
cutval(nn,r,c) = round(drop.l(r,c));
solve mm min total using mip;
mm.solPrint = 0;
objval(nn) = total.l;
);
option cutval:0:1:1;
display objval, cutval;
Parameter hits;
hits(r,c) = sum(nn, cutval(nn,r,c));
display hits;
```