\$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, coordinating relief \$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;