# Mining

## Objective and Prerequisites

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

### What You Will Learn
In this example, we will model and solve a multi-period production planning problem. In this case, the aim of 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, operational, 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 five years to maximize the total profit, while considering production capacity, mine restrictions, and mine costs.

More information on this type of model can be found in the fifth edition of Model Building in Mathematical Programming, by H. Paul Williams.

This modeling example is at the intermediate level, where we assume that you know Python and are familiar with the Gurobi Python API. In addition, you should have some knowledge about building mathematical optimization models.

**Note:** You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). In order to run this Jupyter Notebook properly, you must have a Gurobi license. If you do not have one, you can request an [evaluation license](https://www.gurobi.com/downloads/request-an-evaluation-license/?utm_source=Github&utm_medium=website_JupyterME&utm_campaign=CommercialDataScience) as a *commercial user*, or download a [free license](https://www.gurobi.com/academia/academic-program-and-licenses/?utm_source=Github&utm_medium=website_JupyterME&utm_campaign=AcademicDataScience) as an *academic user*.

---
## 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 still must pay royalties on that mine if there is an expectation that it will operate again in the future. Otherwise, it can be permanently closed and no more royalties need to be paid.

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

| <i></i> | 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:

| <i></i> | Max Production |
| --- | --- |
| Mine 1 | $2.0\times10^6$ Tons |
| Mine 2 | $2.5\times10^6$ Tons |
| Mine 3 | $1.3\times10^6$ Tons |
| Mine 4 | $3.0\times10^6$ 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. The grade for each mine’s ore is as follows:

| <i></i> | 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:

| <i></i> | 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.

The key question for the mining company is: 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 faced by the firm of English China Clays, which had to in decide which mines to work. In that problem (in the 1970s), the goal was to work up to four mines out of 20 in each year.

---
## Model Formulation

### Sets and Indices

$t \in \text{Years}=\{1,2,\dots,5\}$: Set of years.

$m \in \text{Mines}=\{1,2,\dots,4\}$: Set of mines.

### Parameters

$\text{price} \in \mathbb{R}^+$: Selling price (in USD) of one ton of blended ore.

$\text{max_mines} \in \mathbb{N}$: Maximum number of mines that can operate in any given year.

$\text{royalties}_m \in \mathbb{R}^+$: Yearly royalties (in USD) for having mine $m$ open.

$\text{capacity}_m \in \mathbb{R}^+$: Maximum tons of ore that can be extracted from mine $m$ in any given year.

$\text{quality}_m \in \mathbb{R}^+$: Quality of the ore extracted from mine $m$.

$\text{target} \in \mathbb{R}^+$: Quality target of the blended ore in year $t$.

$\text{time_discount}_t \in [0,1] \subset \mathbb{R}^+$: Time discount for revenue and cost in year $t$.

### Decision Variables

$\text{blend}_t \in \mathbb{R}^+$: Tons of blended ore in year $t$.

$\text{extract}_{t,m} \in \mathbb{R}^+$: Tons of ore extracted from mine $m$ in year $t$.

$\text{working}_{t,m} \in \{0,1\}$: 1 if mine $m$ is working in year $t$, 0 otherwise.

$\text{available}_{t,m} \in \{0,1\}$: 1 if mine $m$ is open in year $t$, 0 otherwise.

### Objective Function

- **Profit**: Maximize the total profit (in USD) of the planning horizon.

\begin{equation}
\text{Maximize} \quad Z = \sum_{t \in \text{Years}}\sum_{m \in \text{Mines}}{\text{time_discount}_t*(\text{price}*\text{blend}_t-\text{royalties}_m*\text{extract}_{t,m})}
\tag{0}
\end{equation}

### Constraints

- **Operating Mines**: The total number of operating mines in year $t$ cannot exceed the limit.

\begin{equation}
\sum_{m \in \text{Mines}}{\text{working}_{t,m}} \leq \text{max_mines} \quad \forall t \in \text{Years}
\tag{1}
\end{equation}

- **Quality**: The final quality of the ore blended in year $t$ must meet the target.

\begin{equation}
\sum_{m \in \text{Mines}}{\text{quality}_m*\text{extract}_{t,m}} = \text{target}_t*\text{blended}_t \quad \forall t \in \text{Years}
\tag{2}
\end{equation}

- **Mass Conservation**: Total tons of ore extracted in year $t$ should be equal to the Tons of the ore blended in that year.

\begin{equation}
\sum_{m \in \text{Mines}}{\text{extract}_{t,m}} = \text{blend}_t \quad \forall t \in \text{Years}
\tag{3}
\end{equation}

- **Mine Capacity**: Total tons of ore extracted from mine $m$ in year $t$ cannot exceed the yearly capacity of that mine.

\begin{equation}
\sum_{m \in \text{Mines}}{\text{extract}_{t,m}} \leq \text{capacity}_m*\text{working}_{t,m} \quad \forall t \in \text{Years}
\tag{4}
\end{equation}

- **Open to Operate**: Mine $m$ can be operated in year $t$ only if it is open in that year. 

\begin{equation}
\text{working}_{t,m} \leq \text{available}_{t,m} \quad \forall (t,m) \in \text{Years} \times \text{Mines}
\tag{5}
\end{equation}

- **Shut Down**: If mine $m$ is closed in year $t$, it cannot be opened again in the future.

\begin{equation}
\text{available}_{t+1,m} \leq \text{available}_{t,m} \quad \forall (t < 5,m) \in \text{Years} \times \text{Mines}
\tag{6}
\end{equation}

---
## Python Implementation

We import the Gurobi Python Module and other Python libraries.


```python
import numpy as np
import pandas as pd

import gurobipy as gp
from gurobipy import GRB

# tested with Python 3.7.0 & Gurobi 9.0
```

## Input Data
We define all the input data of the model.


```python
# Parameters

years = [1, 2, 3, 4, 5]
mines = [1, 2, 3, 4]

royalties = {1: 5e6, 2: 4e6, 3: 4e6, 4: 5e6}
capacity = {1: 2e6, 2: 2.5e6, 3: 1.3e6, 4: 3e6}
quality  = {1: 1.0, 2: 0.7, 3: 1.5, 4: 0.5}
target = {1: 0.9, 2: 0.8, 3: 1.2, 4: 0.6, 5: 1.0}
time_discount = {year: (1/(1+1/10.0)) ** (year-1) for year in years}

max_mines = 3
price = 10
```

**Note:** The present value of a future amount at period $n$ is given by the formula: $\text{present_value} = \frac{1}{(1+\text{interest_rate})^n}*\text{future_value} $

## Model Deployment
We create a model and the variables. For each year and each mine, we have (i) a variable that captures production, in millions of tons, (ii) a decision variable which tells us if the mine is open, and (iii) a decision variable which tells us if the mine is operational.


```python
mining = gp.Model('Mining')

blend = mining.addVars(years, name="Blend")
extract = mining.addVars(years, mines, name="Extract")
working = mining.addVars(years, mines, vtype=GRB.BINARY, name="Working")
available = mining.addVars(years, mines, vtype=GRB.BINARY, name="Available")
```

    Using license file c:\gurobi\gurobi.lic
    Set parameter TokenServer to value SANTOS-SURFACE-


Next, we insert the constraints.

In each year only three mines can be operational.


```python
#1. Operating Mines

OperatingMines = mining.addConstrs((working.sum(year, '*') <= max_mines for year in years), "Operating_mines")
```

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


```python
#2. Quality

Quality = mining.addConstrs((gp.quicksum(quality[mine]*extract[year, mine] for mine in mines)
                   == target[year]*blend[year] for year in years), "Quality")
```

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


```python
#3. Mass Conservation

MassConservation = mining.addConstrs((extract.sum(year, '*') == blend[year] for year in years), "Mass Conservation")
```

The following 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 operational.


```python
#4. Mine Capacity

MineCapacity = mining.addConstrs((extract[year, mine] <= capacity[mine]*working[year, mine] for year, mine in extract), "Capacity")
```

The following constraint ensures that when the mine is operational, it also needs to be open.


```python
# Open to operate
OpenToOperate = mining.addConstrs((working[year, mine] <= available[year, mine] for year, mine in available), "Open to Operate")
```

The following 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 cannot be re-opened later:


```python
# Shutdown Mine
ShutdownMine = mining.addConstrs((available[year+1, mine] <= available[year, mine]
                   for year, mine in available if year < years[-1]), "Shut down")
```

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:


```python
#0. Objective function
obj = gp.quicksum(price*time_discount[year]*blend[year] for year in years) \
- gp.quicksum(royalties[mine] * time_discount[year] * available[year, mine] for year, mine in available)
mining.setObjective(obj, GRB.MAXIMIZE)
```

Next, the optimization process starts and Gurobi finds the optimal solution.


```python
mining.optimize()
```

    Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (win64)
    Optimize a model with 71 rows, 65 columns and 182 nonzeros
    Model fingerprint: 0xb3698537
    Variable types: 25 continuous, 40 integer (40 binary)
    Coefficient statistics:
      Matrix range     [5e-01, 3e+06]
      Objective range  [7e+00, 5e+06]
      Bounds range     [1e+00, 1e+00]
      RHS range        [3e+00, 3e+00]
    Found heuristic solution: objective -0.0000000
    Presolve removed 13 rows and 13 columns
    Presolve time: 0.00s
    Presolved: 58 rows, 52 columns, 135 nonzeros
    Variable types: 16 continuous, 36 integer (36 binary)
    
    Root relaxation: objective 1.577309e+08, 40 iterations, 0.00 seconds
    
        Nodes    |    Current Node    |     Objective Bounds      |     Work
     Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
    
         0     0 1.5773e+08    0    4   -0.00000 1.5773e+08      -     -    0s
    H    0     0                    1.189074e+08 1.5773e+08  32.7%     -    0s
    H    0     0                    1.302711e+08 1.5773e+08  21.1%     -    0s
    H    0     0                    1.393620e+08 1.5524e+08  11.4%     -    0s
         0     0 1.5462e+08    0    4 1.3936e+08 1.5462e+08  10.9%     -    0s
    H    0     0                    1.468620e+08 1.5462e+08  5.28%     -    0s
    
    Cutting planes:
      Implied bound: 4
      MIR: 3
      Flow cover: 1
      RLT: 3
      Relax-and-lift: 3
    
    Explored 1 nodes (48 simplex iterations) in 0.06 seconds
    Thread count was 8 (of 8 available processors)
    
    Solution count 5: 1.46862e+08 1.39362e+08 1.30271e+08 ... -0
    
    Optimal solution found (tolerance 1.00e-04)
    Best objective 1.468619743642e+08, best bound 1.468619743642e+08, gap 0.0000%


---
## Analysis

The optimal solution results in a profit of $\$146.8620$ million with the following production plan for each mine in each year (quantities expressed in millions of tons):

### Extraction Plan
This plan determines the millions of tons of ore extracted from each mine (columns)  in each year (rows) of the planning horizon. For example, 1.3 millions of tons of ore will be extracted from mine 3 during the year 2.


```python
rows = years.copy()
columns = mines.copy()
extraction = pd.DataFrame(columns=columns, index=rows, data=0.0)

for year, mine in extract.keys():
    if (abs(extract[year, mine].x) > 1e-6):
        extraction.loc[year, mine] = np.round(extract[year, mine].x / 1e6, 2)
extraction
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1</th>
      <td>2.00</td>
      <td>0.00</td>
      <td>1.3</td>
      <td>2.45</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.00</td>
      <td>2.50</td>
      <td>1.3</td>
      <td>2.20</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1.95</td>
      <td>0.00</td>
      <td>1.3</td>
      <td>0.00</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.12</td>
      <td>2.50</td>
      <td>0.0</td>
      <td>3.00</td>
    </tr>
    <tr>
      <th>5</th>
      <td>2.00</td>
      <td>2.17</td>
      <td>1.3</td>
      <td>0.00</td>
    </tr>
  </tbody>
</table>
</div>



### Sales Plan
This plan defines the millions of tons of blended ore to sell during each year of the planning horizon. For example, we plan to sell 5.62 million tons of blended ore during year 4.


```python
rows = years.copy()
sales = pd.DataFrame(columns=['Sales'], index=rows, data=0.0)

for year in blend.keys():
    if (abs(blend[year].x) > 1e-6):
        sales.loc[year, 'Sales'] = np.round(blend[year].x / 1e6, 2)
sales
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Sales</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1</th>
      <td>5.75</td>
    </tr>
    <tr>
      <th>2</th>
      <td>6.00</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3.25</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5.62</td>
    </tr>
    <tr>
      <th>5</th>
      <td>5.47</td>
    </tr>
  </tbody>
</table>
</div>



**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:

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

---
## References

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

Copyright © 2020 Gurobi Optimization, LLC


```python

```

f8e16027f36a2e340f777b08fcbe4893e2fb6fe2
ipynb
Jupyter Notebook
hp_williams/mining.ipynb
Zhong-HY/modeling-examples
66355938dde764a320f0a6a70a9e5a69c696a18f
Apache-2.0
hp_williams/mining.ipynb
Zhong-HY/modeling-examples
66355938dde764a320f0a6a70a9e5a69c696a18f
Apache-2.0
hp_williams/mining.ipynb
Zhong-HY/modeling-examples
66355938dde764a320f0a6a70a9e5a69c696a18f
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES