$Title StripMine

$Ontext

This is a model for problem 14 in Model Building in Mathematical Programming
by H. Paul Williams.

$Offtext

Sets
   d        "Depths of the mine"			/1*4/
   i        "One coordinate for a parcel of land"	/1*4/
   j        "One coordinate for a parcel of land"	/1*4/
   n(i,j,d) "Neighbors to be removed"
   e(i,j,d) "Indicates if (i,j,d) actually exists as a block" ;

n(i,j,d)=yes$((ord(d) ne card(d)) and (ord(i) le card(d)+1-ord(d))
	and (ord(j) le card(d)+1-ord(d)));
e(i,j,d)=yes$((ord(i) le card(d)+1-ord(d)) and (ord(j) le card(d)+1-ord(d)));

Table val(i,j,d)  "Value for block i,j,d in percent metal"

                    1    2    3    4

          1.1     .75  2.0   5.0  6.0
          1.2     1.0  3.0  12.0 
          1.3     1.5  4.0
          1.4     1.5

          2.1     .75  2.0   4.0
          2.2     1.0  3.0   6.0
          2.3     2.0  4.0
          2.4     1.5

          3.1      .5   .5
          3.2     .75  1.0
          3.3     1.5  2.0
          3.4     1.5

          4.1     .25
          4.2     .5
          4.3     .75
          4.4     .75 ;

Parameter cost(d)
	/ 1 3000
          2 6000
          3 8000
          4 10000/;

Scalar R "Revenue from a 100% block of ore" /200000/;

* The Unimodularity allows x to be continuous and still get binary solutions

Variables
   x(i,j,d)  "binary var indicating whether plot i,j at depth d is mined" 
   profit    "profit variable";

Positive Variables x;
Free Variables profit;

x.up(i,j,d) = 1;

Equations
   Prs1  "Make sure i,j     entry of higher depth is removed"
   Prs2  "Make sure i+1,j   entry of higher depth is removed"
   Prs3  "Make sure i,j+1   entry of higher depth is removed"
   Prs4  "Make sure i+1,j+1 entry of higher depth is removed" 
   Prof  "Total Anticipated Profit" ;

Prs1(i,j,d)$n(i,j,d).. x(i,j,d+1) =l= x(i,j,d);
Prs2(i,j,d)$n(i,j,d).. x(i,j,d+1) =l= x(i+1,j,d);
Prs3(i,j,d)$n(i,j,d).. x(i,j,d+1) =l= x(i,j+1,d);
Prs4(i,j,d)$n(i,j,d).. x(i,j,d+1) =l= x(i+1,j+1,d);

Prof.. profit =e= sum((i,j,d)$e(i,j,d), x(i,j,d)*(val(i,j,d)/100)*R)
	              - sum((i,j,d)$e(i,j,d), x(i,j,d)*cost(d));

Model StripMine /all/;

Solve StripMine using lp maximizing profit;

Display x.l,profit.l;