Mining Example

This model is an example of a production problem. In production planning problems, choices must be made regarding what resources to use to produce what products. These problems are common across a broad range of manufacturing and mining industries.

In this example we’ll model and solve a multi-period production planning problem. In this case, the application is to optimize mine production across a number of mines over a five year period. The company has different mines, which can be opened, working, or closed. Each mine varies in the quality of its ore and the amount of ore that can be extracted. The quality of ore required to meet production goals varies each year. The aim is to create an optimal production plan for the next 5 years to maximize the total profit, while considering all production capacity, mine restrictions and mine costs.

Note: you can download the model, implemented in Python, here. More information on this type of model can be found in the fifth edition of Model Building in Mathematical Programming, by H. Paul Williams.


H. Paul Williams, Model Building in Mathematical Programming, fifth edition (Page 261-262, 357)


Problem Description

A mining company needs to create a five-year operating plan for a certain area with four mines in it.

They can only operate a maximum of three mines in this area in any given year. However, even though a mine may not operate in a given year, the company must still must pay royalties on that mine in that year, if there is the expectation it will operate again in a future year. If a mine is not going to be operated again, then it can be permanently closed down and no more royalties need to be paid.

The yearly royalties due for each open mine (operating or not) are as follows:

  Royalties
Mine 1 $5 million
Mine 2 $4 million
Mine 3 $4 million
Mine 4 $5 million

There is a maximum amount of ore that can be extracted from each mine in a given year. These limits are as follows:

  Maximum Production
Mine 1 2 x 106 tons
Mine 2 2.5 x 106 tons
Mine 3 1.3 x 106 tons
Mine 4 3 x 106 tons

Each mine produces a different grade of ore. This grade is measured on a scale such that blending ores together results in a linear combination of the quality requirements. For example, If equal quantities of ore from two different mines were combined, the resulting ore would have a grade that is the average of the grade for each of the two ores that were combined. The grade for each mine’s ore is as follows:

  Ore Quality
Mine 1 1.0
Mine 2 0.7
Mine 3 1.5
Mine 4 0.5

Each year, the ore produced from each operating mine must be combined to produce ore of a certain grade. The yearly objectives for the combined ore are as follows:

  Produced Ore Quality Target
Year 1 0.9
Year 2 0.8
Year 3 1.2
Year 4 0.6
Year 5 1.0

The final blended ore sells for $10/ton. Revenues and costs for future years are discounted at the rate of 10% per annum.

Which mines should be operated each year and how much ore should be extracted from each mine?

This problem is based on a larger one arising in deciding which pits to work in the firm of English China Clays. In that problem (in the 1970s), their goal was to work up to 4 mines out of 20 in each year.

Model Formulation

Let \( T \) be a set of time periods, where \( t_0 \in T \) is the first period and \( t_e \in T \) the last period. Let \( M\) be a set of mines.

For each time period \(t \in T\) and each skill level \(m \in M\) we introduce a continuous non negative variable \(a_{t,m}\) and binary variables \(w_{t,m}\), \(o_{t,m}\). \(a_{t,m}\) describes the amount, which is mined form mine \(m \) in the time period \(t\). \(w_{t,m}\) describes wether the mine \(m \) is working or not in the time period \(t\). \(o_{t,m}\) describes wether the mine \(m \) is open or not in the time period \(t\).

For each time period \(t \in T\) we introduce a continuous non negative variable \(b_{t}\), which describes the amount of processed oil in each year. For each time period \(t \in T\) we have given the BlendedQuality \(d_t\) and the discount \(e_t\). For each mine \(m \in M\) we have given Royalties \(f_m\), the ExtractLimit \(g_m\) and the OreQuality \(h_m\).

The variables can be formulated as:
\begin{equation} a_{t,m} \geq 0 \quad \forall t \in T, \ \forall m \in M \end{equation} \begin{equation} w_{t,m},o_{t,m} \in \{0,1\} \quad \forall t \in T, \ \forall m \in M \end{equation} \begin{equation} b_{t} \geq 0 \quad \forall t \in T, \end{equation}


In each year only three mines can be working.
\begin{equation} \sum_{m \in M} w_{t,m} \leq 3 \quad \forall t \in T \end{equation}


The quality of the ore from the mines multiplied by the amount which is mined must equal the needed quality multiplied by the quantity of blended ore. This ensures that the quality standards are satisfied.
\begin{equation} \sum_{m \in M} h_mo_{t,m} = d_tb_t \quad \forall t \in T \end{equation}


This constraint ensures that the tonnage of blended ore in each year equals the combined tonnage of the constituents:
\begin{equation} \sum_{m \in M} o_{t,m} - b_t = 0 \quad \forall t \in T \end{equation}


This constraint ensures that the mine can extract no more than the Extract limit and also that there is only an output if the mine is working:
\begin{equation} \sum_{m \in M} o_{t,m} \leq g_mw_{t,m} \quad \forall t \in T, \ \forall m \in M \end{equation}


This constraint ensures that when the mine is working, it also needs to be open:
\begin{equation} \sum_{m \in M} w_{t,m} \leq o_{t,m} \quad \forall t \in T, \ \forall m \in M \end{equation}


This constraint forces a mine to be closed in all years subsequent to that in which it is first closed. If the mine is closed, it can't be re-opened later:
\begin{equation} \sum_{m \in M} o_{t+1,m} \leq o_{t,m} \quad \forall t \in T \setminus t_e, \ \forall m \in M \end{equation}


The total profit consists of the income from selling the blended ore minus the royalties payable. This is to be maximized. It can be written:
\begin{equation} \max \sum_{t \in T}\sum_{m \in M} 10 e_t b_t - f_m e_t o_{t,m} \end{equation}


The full model (with the second objective) can be stated as:
\begin{equation} \max \sum_{t \in T}\sum_{m \in M} 10 e_t b_t - f_m e_t o_{t,m} \end{equation} \begin{equation} \sum_{m \in M} o_{t+1,m} \leq o_{t,m} \quad \forall t \in T \setminus t_e, \ \forall m \in M \end{equation} \begin{equation} \sum_{m \in M} w_{t,m} \leq o_{t,m} \quad \forall t \in T, \ \forall m \in M \end{equation} \begin{equation} \sum_{m \in M} o_{t,m} \leq g_mw_{t,m} \quad \forall t \in T, \ \forall m \in M \end{equation} \begin{equation} \sum_{m \in M} o_{t,m} - b_t = 0 \quad \forall t \in T \end{equation} \begin{equation} \sum_{m \in M} h_mo_{t,m} = d_tb_t \quad \forall t \in T \end{equation} \begin{equation} \sum_{m \in M} w_{t,m} \leq 3 \quad \forall t \in T \end{equation} \begin{equation} a_{t,m} \geq 0 \quad \forall t \in T, \ \forall m \in M \end{equation} \begin{equation} w_{t,m},o_{t,m} \in \{0,1\} \quad \forall t \in T, \ \forall m \in M \end{equation} \begin{equation} b_{t} \geq 0 \quad \forall t \in T, \end{equation}

Implementation with comments

First, we import the Gurobi Python Modue and initalize the datastructures.

from gurobipy import *

# tested with Python 3.5.2 & Gurobi 7.0.1

mines = range(3+1)
years = range(4+1)

Royalties = [5e6, 4e6, 4e6, 5e6]
ExtractLimit = [2e6, 2.5e6, 1.3e6, 3e6]
OreQuality  = [1, .7, 1.5, .5]
BlendedQuality = [0.9, 0.8, 1.2, 0.6, 1.0]
discount = [(1/(1+1/10.0)) ** year for year in years]

mines_limit = 3
sell_price = 10

Next, we create a model and the variables. We set the UpdateMode parameter to 1 (which simplifies the code - see the documentation for more details). For each year and each mine we have a variable for the output in millions of tons the mine produces, a decision variable, which tells us if the mine is open and a decision variable, which tells us if the mine is working.

model = Model('Mining')

out = model.addVars(mines, years, name="output")
quan = model.addVars(years, name="quantity")
work = model.addVars(mines, years, vtype=GRB.BINARY, name="working")
open = model.addVars(mines, years, vtype=GRB.BINARY, name="open")

Next, we insert the constraints. In each year only three mines can be working.

# At most three mines open
model.addConstrs((work.sum('*',year) <= mines_limit for year in years), "AtMost3Mines")

The quality of the ore from the mines multiplied by the amount which is mined must equal the needed blend quality multiplied by the quantity of blended ore. This ensures that the quality standards are satisfied.

# Maintain Quality
model.addConstrs(
    (quicksum(OreQuality[mine]*out[mine, year] for mine in mines) == BlendedQuality[year]*quan[year]
     	for year in years), "Quality")

This constraint ensures that the tonnage of blended ore in each year equals the combined tonnage of the constituents.

# Quantity produced equals output
model.addConstrs((out.sum('*',year) == quan[year] for year in years), "OutQty")

This constraint ensures that the mine can extract no more than the Extract limit and also that there is only an output if there the mine is working.

# Restrict ExtractLimit
#Modeled as described in the HP Williams book
model.addConstrs(
    (out[mine, year] <= ExtractLimit[mine]*work[mine, year]    for mine, year in out), "ExtractLimit")

Alternatively, this constraint can also be modeled using Gurobi 7.0's new general constraints
for year in years:
    for mine in mines:
        out[mine, year].ub= ExtractLimit[mine]
        model. addGenConstrIndicator(work[mine, year], 0, out[mine, year] == 0, name="ExtractLimit" ) 

This constraint ensures that when the mine is working, it also needs to be open.

# Mine Working => Mine Open
#Modeled as described in the HP Williams book
model.addConstrs((work[mine, year] <= open[mine, year] for mine, year in open), "WorkingOpen") 

Modeled using General Constraints

#Modeled using Gurobi General Constraints
for year in years:
    for mine in mines:
        model. addGenConstrIndicator(work[mine, year], 1, open[mine, year] == 1, name="ExtractLimit" ) 

This constraint forces a mine to be closed in all years subsequent to that in which it is first closed. If the mine is closed, it can't be re-opened later:

# Mine Open in Year+1 => Mine was also open in Year
#Modeled as described in the HP Williams book
model.addConstrs(
    (open[mine, year+1] <= open[mine, year]
     for mine, year in open if year < years[-1]), "SubsequentOpen")

Using General Constraints

years2 = (year for year in years if year < years[-1])
for year in years2:
    for mine in mines:
        model. addGenConstrIndicator(open[mine, year + 1], 1, open[mine, year] == 1, name="SubsequentOpen" ) 

The total profit consists of the income from selling the blended ore minus the royalties payable. This is to be maximized. It can be written:

# Maximize Profit
obj = quicksum(sell_price*discount[year]*quan[year] for year in years) \
      - quicksum(Royalties[mine] * discount[year] * open[mine, year]
                 for mine, year in open)
model.setObjective(obj, GRB.MAXIMIZE)

Next, this starts the optimization and Gurobi tries to find the optimal solution.

model.optimize()
for v in model.getVars():
    if v.X != 0:
        print("%s %f" % (v.Varname, v.X))

Note: If you want to write your solution to a file, rather than print it to the terminal, you can use the model.write() command. An example implementation is:

 model.write("mining-output.sol")

Gurobi Output and Analysis

Optimize a model with 15 rows, 65 columns and 70 nonzeros
Model has 56 general constraints
Variable types: 25 continuous, 40 integer (40 binary)
Coefficient statistics:
  Matrix range     [5e-01, 2e+00]
  Objective range  [7e+00, 5e+06]
  Bounds range     [1e+00, 3e+06]
  RHS range        [3e+00, 3e+00]
Found heuristic solution: objective -7.50576e+07
Presolve removed 25 rows and 25 columns
Presolve time: 0.00s
Presolved: 66 rows, 116 columns, 151 nonzeros
Presolved model has 56 SOS constraint(s)
Variable types: 20 continuous, 96 integer (60 binary)

Root relaxation: objective 2.543811e+08, 23 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 2.5438e+08    0   16 -7.506e+07 2.5438e+08   439%     -    0s
H    0     0                      -0.0000000 2.5438e+08      -     -    0s
H    0     0                    6.290909e+07 2.5438e+08   304%     -    0s
H    0     0                    6.681632e+07 2.5438e+08   281%     -    0s
     0     2 2.5438e+08    0   16 6.6816e+07 2.5438e+08   281%     -    0s
*   56    48              16    1.297493e+08 2.3688e+08  82.6%   0.2    0s
*  129    80              16    1.306469e+08 2.0971e+08  60.5%   0.2    0s
*  131    74              17    1.372493e+08 2.0971e+08  52.8%   0.2    0s
*  340    98              16    1.397378e+08 1.9736e+08  41.2%   0.1    0s
*  342    90              17    1.463402e+08 1.9736e+08  34.9%   0.1    0s
*  730    48              15    1.468620e+08 1.6573e+08  12.8%   0.1    0s

Explored 799 nodes (140 simplex iterations) in 0.06 seconds
Thread count was 4 (of 4 available processors)

Solution count 10: 1.46862e+08 1.4634e+08 1.39738e+08 ... -7.50576e+07
Pool objective bound 1.46862e+08

Optimal solution found (tolerance 1.00e-04)
Best objective 1.468619743642e+08, best bound 1.468619743642e+08, gap 0.0000%

Analyzing the Output

The optimal solution results in a profit of $146.86 million with the following production plan for each mine in each year:

  Mines
Year 1 1, 3, 4
Year 2 2, 3, 4
Year 3 1, 3
Year 4 1, 2, 4
Year 5 1, 2, 3

The plan indicates that each mine should be kept open each year, with the exception of mine four in year five. In addition, each mine should produce the following quantities of ore (in millions of tons):

  Mine 1 Mine 2 Mine 3 Mine 4
Year 1 2.0 - 1.3 2.45
Year 2 - 2.5 1.3 2.2
Year 3 1.95 - 1.3 -
Year 4 0.12 2.5 - 3.0
Year 5 2.0 2.17 1.3 -

This results in the production of the following quantities of blended ore (in millions of tons) each year:

  Total
Year 1 5.75
Year 2 6.00
Year 3 3.25
Year 4 5.62
Year 5 5.47