MineLib in AMPL and amplpy#
Description: A notebook that works with the open-pit mining problems defined in the MineLib library and shows how to read and parse the instances using Python and amplpy
Tags: educational, mining, amplpy, gurobi, minelib
Notebook author: Eduardo Salazar <eduardo@ampl.com>
Model author: Daniel Espinoza (Universidad de Chile), Marcos Goycoolea (Universidad Adolfo Ibañez), Eduardo Moreno (Universidad Adolfo Ibañez), Alexandra Newman (Colorado School of Mines).
References:
Help on Minelib: https://mansci-web.uai.cl/minelib/Help.xhtml
MineLib: A Library of Open Pit Mining Problems. Daniel Espinoza, Marcos Goycoolea, Eduardo Moreno, Alexandra Newman https://mgoycool.github.io/papers/13espinoza_aor.pdf.
# Install dependencies
%pip install -q amplpy requests
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["gurobi"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
1. Introduction to Open Pit Mining #
Open pit mining is a surface mining technique used to extract minerals from the earth by removing layers of rock and soil. The extracted material is discretized into blocks, each containing estimated amounts of ore (valuable mineral) and waste. The key decisions involve determining which blocks to extract and when, subject to various operational and geometric constraints.
Key Concepts#
Block Model: The mineral deposit is represented as a three-dimensional grid of blocks. Each block has associated attributes such as position \((x, y, z)\), tonnage, ore grade, and estimated economic value.
Precedence Constraints: Due to geotechnical requirements for pit wall stability and equipment access, blocks cannot be extracted in arbitrary order. If block \(a\) lies above block \(b\), then \(a\) must be extracted before \(b\). These relationships form a directed acyclic graph where edges represent “must be extracted before” relationships.
Pit Slope: The angle at which the pit walls are maintained to prevent collapse. This angle determines which blocks must be removed to access deeper blocks.
Cutoff Grade: The minimum ore concentration that makes a block economically viable for processing. Blocks below this threshold are typically sent to waste dumps.
Why Optimization?#
Mining operations represent massive capital investments spanning decades. A typical open pit mine may contain millions of blocks, and the extraction sequence significantly impacts profitability due to the time value of money. Early extraction of high-value ore maximizes net present value (NPV). Mathematical optimization helps mine planners make informed decisions about pit design and production scheduling.
2. The MineLib Library#
MineLib was introduced by Espinoza, Goycoolea, Moreno, and Newman (2012) as a standardized library of open pit mining problem instances. Similar to how MIPLIB serves the mixed-integer programming community, MineLib provides benchmark data for evaluating and comparing optimization algorithms.
Available Instances#
Instance |
Blocks |
Precedences |
Periods |
Description |
|---|---|---|---|---|
Newman1 |
1,060 |
3,922 |
6 |
Small academic dataset |
Zuck_small |
9,400 |
145,640 |
20 |
Fictitious mine |
D |
14,153 |
219,778 |
12 |
North American copper deposit |
Zuck_medium |
29,277 |
1,271,207 |
15 |
Fictitious mine |
P4HD |
40,947 |
738,609 |
10 |
North American gold/copper mine |
Marvin |
53,271 |
650,631 |
20 |
Whittle software test mine |
W23 |
74,260 |
764,786 |
12 |
North American gold mine |
Zuck_large |
96,821 |
1,053,105 |
30 |
Fictitious mine |
SM2 |
99,014 |
96,642 |
30 |
Based on Brazilian nickel mine |
McLaughlin_limit |
112,687 |
3,035,483 |
15 |
California gold mine (final pit) |
McLaughlin |
2,140,342 |
73,143,770 |
20 |
California gold mine (full) |
Problem Types#
MineLib supports three problem formulations of increasing complexity:
UPIT (Ultimate Pit Limit): Determines the optimal pit boundary ignoring time.
CPIT (Constrained Pit Limit): Adds time periods and resource constraints with fixed cutoff grades.
PCPSP (Precedence Constrained Production Scheduling): Variable cutoff grades with multiple destinations and general side constraints.
3. Mathematical Formulations#
We now present the mathematical formulations for the three problem types. The following notation is used throughout:
3.1 Notation#
Indices and Sets:
Symbol |
Description |
|---|---|
\(t \in \mathcal{T}\) |
Set of time periods \(t\) in the horizon |
\(b \in \mathcal{B}\) |
Set of blocks \(b\) |
\(b' \in \mathcal{B}_b\) |
Set of blocks \(b'\) that are predecessor blocks for block \(b\) |
\(r \in \mathcal{R}\) |
Set of operational resource types \(r\) |
\(d \in \mathcal{D}\) |
Set of destinations \(d\) |
Parameters:
Symbol |
Description |
|---|---|
\(p_b\) (\(\hat{p}_{bt}\), \(\check{p}_{bd}\), \(\tilde{p}_{bdt}\)) |
Profit obtained from extracting (and processing) block \(b\) (at time period \(t\) and/or sending it to destination \(d\)) ($) |
\(\alpha\) |
Discount rate used in computing the objective function (profit) coefficients |
\(q_{br}\) (\(\hat{q}_{brd}\)) |
Amount of operational resource \(r\) used to extract and, if applicable, process block \(b\) (when sent to destination \(d\)) (tons) |
\(\underline{R}_{rt}\) |
Minimum availability of operational resource \(r\) in time period \(t\) (tons) |
\(\overline{R}_{rt}\) |
Maximum availability of operational resource \(r\) in time period \(t\) (tons) |
\(A\) |
Arbitrary constraint coefficients on general side constraints |
\(\underline{a}, \overline{a}\) |
Arbitrary lower and upper bounds, respectively, on general side constraints |
Variables:
Symbol |
Description |
|---|---|
\(\hat{x}_b\) |
1 if block \(b\) is in the final pit design, 0 otherwise |
\(x_{bt}\) |
1 if we extract block \(b\) in time period \(t\), 0 otherwise |
\(y_{bdt}\) |
The amount of block \(b\) sent to destination \(d\) in time period \(t\) (%) |
3.2 The Ultimate Pit Problem (UPIT)#
The simplest model we consider is known as the ultimate pit limit problem, or the maximum-weight closure problem. The problem entails determining only the envelope of profitable blocks within the orebody and, hence, there is no temporal dimension and there are no operational resource constraints. The constraint set consists merely of precedences between blocks; the corresponding matrix of left-hand-side coefficients is totally unimodular, rendering this problem a network flow problem.
Formulation:
\[
\text{(UPIT)} \quad \max \sum_{b \in \mathcal{B}} p_b \, \hat{x}_b
\]
subject to:
\[
\hat{x}_b \leq \hat{x}_{b'} \quad \forall b \in \mathcal{B}, \; \forall b' \in \mathcal{B}_b \tag{1}
\]
\[
\hat{x}_b \in \{0, 1\} \quad \forall b \in \mathcal{B} \tag{2}
\]
Notes:
The objective maximizes the undiscounted value of all extracted blocks.
Constraints (1) ensure that each block is extracted only if its predecessor blocks are extracted.
Variables need not be restricted to be binary because of the total unimodularity of the constraint matrix.
The Lerchs-Grossmann algorithm (1965) and successors provide efficient specialized solutions.
3.3 The Constrained Pit Limit Problem (CPIT)#
The constrained pit limit problem generalizes the ultimate pit limit problem by introducing a time dimension, and associated constraints, into the model. The underlying assumption is that a block can be mined in its entirety in a single time period. Not only are precedence constraints considered, but per-period operational resource restrictions are present as well.
Formulation:
\[
\text{(CPIT)} \quad \max \sum_{b \in \mathcal{B}} \sum_{t \in \mathcal{T}} \hat{p}_{bt} \, x_{bt}
\]
subject to:
\[
\sum_{\tau \leq t} x_{b\tau} \leq \sum_{\tau \leq t} x_{b'\tau} \quad \forall b \in \mathcal{B}, \; b' \in \mathcal{B}_b, \; t \in \mathcal{T} \tag{3}
\]
\[
\sum_{t \in \mathcal{T}} x_{bt} \leq 1 \quad \forall b \in \mathcal{B} \tag{4}
\]
\[
\underline{R}_{rt} \leq \sum_{b \in \mathcal{B}} q_{br} \, x_{bt} \leq \overline{R}_{rt} \quad \forall t \in \mathcal{T}, \; r \in \mathcal{R} \tag{5}
\]
\[
x_{bt} \in \{0, 1\} \quad \forall b \in \mathcal{B}, \; t \in \mathcal{T} \tag{6}
\]
where \(\hat{p}_{bt} = \dfrac{p_b}{(1+\alpha)^t}\).
Constraint Interpretation:
Precedence (3): If block \(b'\) is an immediate predecessor of block \(b\), then \(b'\) must be extracted in the same time period as or prior to \(b\).
Single extraction (4): Each block can be extracted no more than once.
Resource bounds (5): Minimum and maximum operational resource constraints are satisfied each period.
Complexity: Unlike UPIT, CPIT is strongly NP-hard.
3.4 The Precedence Constrained Production Scheduling Problem (PCPSP)#
A generalization of CPIT determines whether a block, if extracted, is sent to the processing plant or to the waste dump. In addition to variable \(x_{bt}\), we employ a second variable \(y_{bdt}\), which equals the amount of block \(b\) we send to destination \(d\) in time period \(t\). This enables variable cutoff grades where the destination is optimized rather than predetermined.
Formulation:
\[
\text{(PCPSP)} \quad \max \sum_{b \in \mathcal{B}} \sum_{d \in \mathcal{D}} \sum_{t \in \mathcal{T}} \tilde{p}_{bdt} \, y_{bdt}
\]
subject to:
\[
\sum_{\tau \leq t} x_{b\tau} \leq \sum_{\tau \leq t} x_{b'\tau} \quad \forall b \in \mathcal{B}, \; b' \in \mathcal{B}_b, \; t \in \mathcal{T} \tag{7}
\]
\[
x_{bt} = \sum_{d \in \mathcal{D}} y_{bdt} \quad \forall b \in \mathcal{B}, \; t \in \mathcal{T} \tag{8}
\]
\[
\sum_{t \in \mathcal{T}} x_{bt} \leq 1 \quad \forall b \in \mathcal{B} \tag{9}
\]
\[
\underline{R}_{rt} \leq \sum_{b \in \mathcal{B}} \sum_{d \in \mathcal{D}} \hat{q}_{brd} \, y_{bdt} \leq \overline{R}_{rt} \quad \forall r \in \mathcal{R}, \; t \in \mathcal{T} \tag{10}
\]
\[
\underline{a} \leq A y \leq \overline{a} \tag{11}
\]
\[
y_{bdt} \in [0, 1] \quad \forall b \in \mathcal{B}, \; d \in \mathcal{D}, \; t \in \mathcal{T} \tag{12}
\]
\[
x_{bt} \in \{0, 1\} \quad \forall b \in \mathcal{B}, \; t \in \mathcal{T} \tag{13}
\]
where \(\tilde{p}_{bdt} = \dfrac{\check{p}_{bd}}{(1+\alpha)^t}\).
Constraint Interpretation:
Precedence (7): Enforces precedence requirements for all blocks and time periods.
Linking (8): If a block is not extracted, its contents cannot be sent to any destination; if extracted, the entirety must be sent somewhere.
Single extraction (9): A block can be extracted at most once over the horizon.
Resource bounds (10): No more operational resource than available is used for extraction purposes.
General side constraints (11): Can model blending, grade constraints, stockpiling, etc.
Example - Grade Constraint:
Letting \(g_{bm}\) be the amount of mineral \(m\) in block \(b\), and \(\underline{G}_m\), \(\overline{G}_m\) be the minimum/maximum acceptable average grades:
\[
\underline{G}_m \sum_{b \in \mathcal{B}} y_{bdt} \leq \sum_{b \in \mathcal{B}} g_{bm} \, y_{bdt} \leq \overline{G}_m \sum_{b \in \mathcal{B}} y_{bdt} \quad \forall d \in \mathcal{D}, \; m \in \mathcal{M}, \; t \in \mathcal{T} \tag{14}
\]
4. MineLib File Format#
MineLib instances are distributed in three file types:
4.1 Block-Model Descriptor File (.blocks)#
Contains geological data for each block:
<id> <x> <y> <z> [additional fields...]
id: Unique block identifier (0-indexed)x, y, z: Block coordinates (z=0 is the bottom)Additional fields may include tonnage, grade, etc.
4.2 Block-Precedence Descriptor File (.prec)#
Defines precedence relationships:
<block_id> <num_predecessors> <pred_1> <pred_2> ... <pred_n>
Example:
100 5 10 11 12 20 21
Block 100 has 5 predecessors: blocks 10, 11, 12, 20, and 21 must be extracted first.
4.3 Optimization-Model Descriptor File (.upit, .cpit, .pcpsp)#
Contains optimization parameters:
NAME: instance_name
TYPE: UPIT|CPIT|PCPSP
NBLOCKS: n
NPERIODS: T (CPIT/PCPSP only)
NDESTINATIONS: D (PCPSP only)
NRESOURCE_SIDE_CONSTRAINTS: R (CPIT/PCPSP only)
DISCOUNT_RATE: alpha (CPIT/PCPSP only)
OBJECTIVE_FUNCTION:
<block_id> <value_dest_0> [<value_dest_1> ...]
...
RESOURCE_CONSTRAINT_LIMITS:
<resource> <period> <type> <bound1> [<bound2>]
...
RESOURCE_CONSTRAINT_COEFFICIENTS:
<block> [<dest>] <resource> <coefficient>
...
EOF
The constraint type can be:
L: Less than or equal (\(\leq\) bound1)G: Greater than or equal (\(\geq\) bound1)I: Interval (bound1 \(\leq \cdot \leq\) bound2)
In the hidden code cells we define the MineLib parser tools for instances and solutions in Python.
For this notebook, we will be solving the smallest instance which is newman1.
from pathlib import Path
from typing import Dict, List, Optional, Tuple, Union, TextIO
from dataclasses import dataclass, field
import math
import requests
from io import StringIO
from contextlib import contextmanager
# @title
@dataclass
class Precedences:
"""Block precedence relationships."""
num_blocks: int
predecessors: Dict[int, List[int]] = field(default_factory=dict)
def get_predecessors(self, block_id: int) -> List[int]:
return self.predecessors.get(block_id, [])
def num_precedence_arcs(self) -> int:
return sum(len(preds) for preds in self.predecessors.values())
@dataclass
class ResourceLimits:
"""Resource constraint bounds by resource and time period."""
lower_bounds: Dict[int, Dict[int, float]] = field(default_factory=dict)
upper_bounds: Dict[int, Dict[int, float]] = field(default_factory=dict)
@dataclass
class UPITData:
"""Data for Ultimate Pit Limit Problem."""
name: str
num_blocks: int
objective: Dict[int, float] = field(default_factory=dict)
precedences: Optional[Precedences] = None
@dataclass
class CPITData:
"""Data for Constrained Pit Limit Problem."""
name: str
num_blocks: int
num_periods: int
num_resources: int
discount_rate: float
objective: Dict[int, float] = field(default_factory=dict)
resource_limits: ResourceLimits = field(default_factory=ResourceLimits)
resource_coefficients: Dict[int, Dict[int, float]] = field(default_factory=dict)
precedences: Optional[Precedences] = None
@dataclass
class PCPSPData:
"""Data for Precedence Constrained Production Scheduling Problem."""
name: str
num_blocks: int
num_periods: int
num_destinations: int
num_resources: int
num_general_constraints: int
discount_rate: float
objective: Dict[int, Dict[int, float]] = field(default_factory=dict)
resource_limits: ResourceLimits = field(default_factory=ResourceLimits)
resource_coefficients: Dict[int, Dict[int, Dict[int, float]]] = field(
default_factory=dict
)
general_coefficients: Dict[int, Dict[int, Dict[int, Dict[int, float]]]] = field(
default_factory=dict
)
general_limits: Dict[int, Tuple[float, float]] = field(default_factory=dict)
precedences: Optional[Precedences] = None
# @title
def parse_precedences(file_source: Union[str, Path, TextIO]) -> Precedences:
"""
Parse a MineLib precedence file.
Args:
file_source: Either a file path (str/Path) or a file-like object (StringIO)
"""
predecessors: Dict[int, List[int]] = {}
num_blocks = 0
# Handle both file paths and file-like objects
if isinstance(file_source, (str, Path)):
f = open(file_source, "r")
should_close = True
else:
f = file_source
should_close = False
try:
for line in f:
line = line.strip()
if not line or line.startswith("%"):
continue
parts = line.split()
block_id = int(parts[0])
num_preds = int(parts[1])
preds = [int(parts[i]) for i in range(2, 2 + num_preds)]
predecessors[block_id] = preds
num_blocks = max(num_blocks, block_id + 1)
finally:
if should_close:
f.close()
return Precedences(num_blocks=num_blocks, predecessors=predecessors)
def _parse_float(value: str) -> float:
"""Parse float, handling infinity."""
value = value.strip().lower()
if value in ("infinity", "inf", "+infinity", "+inf"):
return math.inf
if value in ("-infinity", "-inf"):
return -math.inf
return float(value)
@contextmanager
def open_source(source: Union[str, Path, TextIO]):
"""Context manager that handles both file paths and file-like objects."""
if isinstance(source, (str, Path)):
f = open(source, "r")
try:
yield f
finally:
f.close()
else:
# Already a file-like object (StringIO), just yield it
yield source
def parse_upit(
opt_source: Union[str, Path, TextIO],
prec_source: Optional[Union[str, Path, TextIO, Precedences]] = None,
) -> UPITData:
"""Parse a UPIT optimization file."""
data = UPITData(name="", num_blocks=0)
# Use context manager instead of Path()
with open_source(opt_source) as f:
for line in f:
line = line.strip()
if not line or line.startswith("%"):
continue
if ":" in line:
key, value = line.split(":", 1)
key, value = key.strip().upper(), value.strip()
if key == "NAME":
data.name = value
elif key == "NBLOCKS":
data.num_blocks = int(value)
elif key == "OBJECTIVE_FUNCTION":
for _ in range(data.num_blocks):
obj_line = f.readline().strip()
while not obj_line or obj_line.startswith("%"):
obj_line = f.readline().strip()
parts = obj_line.split()
block_id = int(parts[0])
data.objective[block_id] = _parse_float(parts[1])
# Handle precedences flexibly
if prec_source is not None:
if isinstance(prec_source, Precedences):
data.precedences = prec_source
else:
data.precedences = parse_precedences(prec_source)
return data
def parse_cpit(
opt_source: Union[str, Path, TextIO],
prec_source: Optional[Union[str, Path, TextIO, Precedences]] = None,
) -> CPITData:
"""Parse a CPIT optimization file."""
data = CPITData(
name="", num_blocks=0, num_periods=0, num_resources=0, discount_rate=0.0
)
with open_source(opt_source) as f:
for line in f:
line = line.strip()
if not line or line.startswith("%"):
continue
if ":" in line:
key, value = line.split(":", 1)
key, value = key.strip().upper().replace(" ", "_"), value.strip()
if key == "NAME":
data.name = value
elif key == "NBLOCKS":
data.num_blocks = int(value)
elif key == "NPERIODS":
data.num_periods = int(value)
elif key in (
"NRESOURCE_SIDE_CONSTRAINTS",
"NRESOURCE_SIDE_CONSTRAINTS",
):
data.num_resources = int(value)
elif key == "DISCOUNT_RATE":
data.discount_rate = _parse_float(value)
elif key == "OBJECTIVE_FUNCTION":
for _ in range(data.num_blocks):
obj_line = f.readline().strip()
while not obj_line or obj_line.startswith("%"):
obj_line = f.readline().strip()
parts = obj_line.split()
data.objective[int(parts[0])] = _parse_float(parts[1])
elif key == "RESOURCE_CONSTRAINT_LIMITS":
data.resource_limits.lower_bounds = {
r: {} for r in range(data.num_resources)
}
data.resource_limits.upper_bounds = {
r: {} for r in range(data.num_resources)
}
for _ in range(data.num_resources * data.num_periods):
lim_line = f.readline().strip()
while not lim_line or lim_line.startswith("%"):
lim_line = f.readline().strip()
parts = lim_line.split()
r, t, ctype = int(parts[0]), int(parts[1]), parts[2].upper()
if ctype == "L":
data.resource_limits.lower_bounds[r][t] = -math.inf
data.resource_limits.upper_bounds[r][t] = _parse_float(
parts[3]
)
elif ctype == "G":
data.resource_limits.lower_bounds[r][t] = _parse_float(
parts[3]
)
data.resource_limits.upper_bounds[r][t] = math.inf
elif ctype == "I":
data.resource_limits.lower_bounds[r][t] = _parse_float(
parts[3]
)
data.resource_limits.upper_bounds[r][t] = _parse_float(
parts[4]
)
elif key == "RESOURCE_CONSTRAINT_COEFFICIENTS":
while True:
coef_line = f.readline()
if not coef_line:
break
coef_line = coef_line.strip()
if not coef_line or coef_line.startswith("%"):
continue
if coef_line.upper() == "EOF":
break
parts = coef_line.split()
if len(parts) < 3:
break
block_id, resource_id = int(parts[0]), int(parts[1])
coef_value = _parse_float(parts[2])
if block_id not in data.resource_coefficients:
data.resource_coefficients[block_id] = {}
data.resource_coefficients[block_id][resource_id] = coef_value
# Handle precedences flexibly
if prec_source is not None:
if isinstance(prec_source, Precedences):
data.precedences = prec_source
else:
data.precedences = parse_precedences(prec_source)
return data
def parse_pcpsp(
opt_source: Union[str, Path, TextIO],
prec_source: Optional[Union[str, Path, TextIO, Precedences]] = None,
) -> PCPSPData:
"""Parse a PCPSP optimization file."""
data = PCPSPData(
name="",
num_blocks=0,
num_periods=0,
num_destinations=0,
num_resources=0,
num_general_constraints=0,
discount_rate=0.0,
)
with open_source(opt_source) as f:
for line in f:
line = line.strip()
if not line or line.startswith("%"):
continue
if ":" in line:
key, value = line.split(":", 1)
key, value = key.strip().upper().replace(" ", "_"), value.strip()
if key == "NAME":
data.name = value
elif key == "NBLOCKS":
data.num_blocks = int(value)
elif key == "NPERIODS":
data.num_periods = int(value)
elif key == "NDESTINATIONS":
data.num_destinations = int(value)
elif key in ("NRESOURCE_SIDE_CONSTRAINTS",):
data.num_resources = int(value)
elif key in ("NGENERAL_SIDE_CONSTRAINTS",):
data.num_general_constraints = int(value)
elif key == "DISCOUNT_RATE":
data.discount_rate = _parse_float(value)
elif key == "OBJECTIVE_FUNCTION":
for _ in range(data.num_blocks):
obj_line = f.readline().strip()
while not obj_line or obj_line.startswith("%"):
obj_line = f.readline().strip()
parts = obj_line.split()
block_id = int(parts[0])
data.objective[block_id] = {}
for d in range(data.num_destinations):
data.objective[block_id][d] = _parse_float(parts[1 + d])
elif key == "RESOURCE_CONSTRAINT_LIMITS":
data.resource_limits.lower_bounds = {
r: {} for r in range(data.num_resources)
}
data.resource_limits.upper_bounds = {
r: {} for r in range(data.num_resources)
}
for _ in range(data.num_resources * data.num_periods):
lim_line = f.readline().strip()
while not lim_line or lim_line.startswith("%"):
lim_line = f.readline().strip()
parts = lim_line.split()
r, t, ctype = int(parts[0]), int(parts[1]), parts[2].upper()
if ctype == "L":
data.resource_limits.lower_bounds[r][t] = -math.inf
data.resource_limits.upper_bounds[r][t] = _parse_float(
parts[3]
)
elif ctype == "G":
data.resource_limits.lower_bounds[r][t] = _parse_float(
parts[3]
)
data.resource_limits.upper_bounds[r][t] = math.inf
elif ctype == "I":
data.resource_limits.lower_bounds[r][t] = _parse_float(
parts[3]
)
data.resource_limits.upper_bounds[r][t] = _parse_float(
parts[4]
)
elif key == "RESOURCE_CONSTRAINT_COEFFICIENTS":
while True:
coef_line = f.readline()
if not coef_line:
break
coef_line = coef_line.strip()
if not coef_line or coef_line.startswith("%"):
continue
if coef_line.upper() == "EOF":
break
parts = coef_line.split()
if len(parts) < 4:
break
block_id, dest_id = int(parts[0]), int(parts[1])
resource_id, coef_value = int(parts[2]), _parse_float(parts[3])
if block_id not in data.resource_coefficients:
data.resource_coefficients[block_id] = {}
if dest_id not in data.resource_coefficients[block_id]:
data.resource_coefficients[block_id][dest_id] = {}
data.resource_coefficients[block_id][dest_id][
resource_id
] = coef_value
# Handle precedences flexibly
if prec_source is not None:
if isinstance(prec_source, Precedences):
data.precedences = prec_source
else:
data.precedences = parse_precedences(prec_source)
return data
UPIT_MODEL = """
# Ultimate Pit Limit Problem (UPIT)
# Maximum-weight closure problem - polynomial time solvable
param NBLOCKS;
set setBL = 0..NBLOCKS-1;
param OBJECTIVE_FUNCTION {setBL};
set PREC{setBL};
var X{setBL} binary;
maximize Profit: sum{b in setBL} OBJECTIVE_FUNCTION[b] * X[b];
subject to Precedence {b in setBL, i in PREC[b]}:
X[b] <= X[i];
"""
CPIT_MODEL = """
# Constrained Pit Limit Problem (CPIT)
# NP-hard - maximizes NPV with resource constraints
param NBLOCKS;
param NPERIODS;
param NRESOURCE_SIDE_CONSTRAINTS;
param DISCOUNT_RATE;
set setBL = 0..NBLOCKS-1;
set setT = 0..NPERIODS-1;
set setSIDE = 0..NRESOURCE_SIDE_CONSTRAINTS-1;
param RESOURCE_CONSTRAINT_UB_LIMITS {setSIDE,setT};
param RESOURCE_CONSTRAINT_LB_LIMITS {setSIDE,setT};
param RESOURCE_CONSTRAINT_COEFFICIENTS {setBL, setSIDE} default 0;
param OBJECTIVE_FUNCTION {setBL};
set PREC{setBL};
var X{setBL,setT} binary;
maximize Profit: sum{t in setT} sum{b in setBL} (1.0/(1.0 + DISCOUNT_RATE))^t * OBJECTIVE_FUNCTION[b] * X[b,t];
subject to Resource { r in setSIDE, t in setT}:
RESOURCE_CONSTRAINT_LB_LIMITS[r,t] <= sum{b in setBL} RESOURCE_CONSTRAINT_COEFFICIENTS[b,r] * X[b,t] <= RESOURCE_CONSTRAINT_UB_LIMITS[r,t];
subject to CliqueBlock {b in setBL}:
sum{t in setT} X[b,t] <= 1;
subject to Precedence {b in setBL, i in PREC[b], t in setT}:
sum{s in 0..t} X[b,s] <= sum{s in 0..t} X[i,s];
"""
PCPSP_MODEL = """
# Precedence Constrained Production Scheduling Problem (PCPSP)
# Variable cutoff grade with multiple destinations
param NBLOCKS;
param NPERIODS;
param NDESTINATIONS;
param NRESOURCE_SIDE_CONSTRAINTS;
param DISCOUNT_RATE;
set setBL = 0..NBLOCKS-1;
set setT = 0..NPERIODS-1;
set setSIDE = 0..NRESOURCE_SIDE_CONSTRAINTS-1;
set setDEST = 0..NDESTINATIONS-1;
param RESOURCE_CONSTRAINT_UB_LIMITS {setSIDE,setT};
param RESOURCE_CONSTRAINT_LB_LIMITS {setSIDE,setT};
param RESOURCE_CONSTRAINT_COEFFICIENTS {setBL, setSIDE, setDEST} default 0;
param OBJECTIVE_FUNCTION {setBL,setDEST};
set PREC{setBL};
var X{setBL,setT} binary;
var Y{setBL,setDEST,setT} >=0, <=1;
maximize Profit: sum{t in setT} sum {d in setDEST} sum{b in setBL} (1.0/(1.0 + DISCOUNT_RATE))^t * OBJECTIVE_FUNCTION[b,d] * Y[b,d,t];
subject to Resource { r in setSIDE, t in setT}:
RESOURCE_CONSTRAINT_LB_LIMITS[r,t] <= sum{b in setBL} sum {d in setDEST} RESOURCE_CONSTRAINT_COEFFICIENTS[b,r,d] * Y[b,d,t] <= RESOURCE_CONSTRAINT_UB_LIMITS[r,t];
subject to CliqueBlock {b in setBL}:
sum{t in setT} X[b,t] <= 1;
subject to SumDest {b in setBL, t in setT}:
sum {d in setDEST} Y[b,d,t] = X[b,t];
subject to Precedence {b in setBL, i in PREC[b], t in setT}:
sum{s in 0..t} X[b,s] <= sum{s in 0..t} X[i,s];
"""
# @title
@dataclass
class SolutionUPIT:
"""Solution for UPIT problem."""
objective: float
solve_time: float
status: str
extracted: Dict[int, bool] # block_id -> in pit?
def num_extracted(self) -> int:
return sum(1 for v in self.extracted.values() if v)
@dataclass
class SolutionCPIT:
"""Solution for CPIT problem."""
objective: float
solve_time: float
status: str
extraction_period: Dict[
int, Optional[int]
] # block_id -> period (None if not extracted)
def num_extracted(self) -> int:
return sum(1 for v in self.extraction_period.values() if v is not None)
def blocks_in_period(self, t: int) -> List[int]:
return [b for b, period in self.extraction_period.items() if period == t]
@dataclass
class SolutionPCPSP:
"""Solution for PCPSP problem."""
objective: float
solve_time: float
status: str
extraction_period: Dict[int, Optional[int]]
destination_allocation: Dict[
int, Dict[int, float]
] # block_id -> {dest_id: fraction}
def num_extracted(self) -> int:
return sum(1 for v in self.extraction_period.values() if v is not None)
# @title
def _format_bound(value: float) -> float:
"""Convert infinity to large number for AMPL."""
if value == math.inf:
return 1e20
elif value == -math.inf:
return -1e20
return value
def solve_upit(
data: UPITData,
solver: str = "gurobi",
solver_options: str = "",
verbose: bool = True,
) -> SolutionUPIT:
"""Solve a UPIT problem."""
if data.precedences is None:
raise ValueError("Precedences must be loaded")
ampl = AMPL()
ampl.eval(UPIT_MODEL)
# Load parameters
ampl.param["NBLOCKS"] = data.num_blocks
for b in range(data.num_blocks):
ampl.param["OBJECTIVE_FUNCTION"][b] = data.objective.get(b, 0.0)
# Load precedences as indexed sets
for b in range(data.num_blocks):
ampl.set["PREC"][b] = data.precedences.get_predecessors(b)
# Configure and solve
ampl.option["solver"] = solver
if solver_options:
ampl.option[f"{solver}_options"] = solver_options
if verbose:
print(f"Solving UPIT: {data.name}")
print(
f" Blocks: {data.num_blocks}, Arcs: {data.precedences.num_precedence_arcs()}"
)
ampl.solve()
# Extract solution
status = ampl.get_value("solve_result")
solve_time = ampl.get_value("_solve_time")
objective = ampl.obj["Profit"].value()
extracted = {}
X = ampl.var["X"]
for b in range(data.num_blocks):
extracted[b] = X[b].value() > 0.5
if verbose:
print(
f" Status: {status}, Objective: ${objective:,.2f}, Time: {solve_time:.2f}s"
)
return SolutionUPIT(
objective=objective, solve_time=solve_time, status=status, extracted=extracted
)
def solve_cpit(
data: CPITData,
solver: str = "gurobi",
solver_options: str = "",
verbose: bool = True,
) -> SolutionCPIT:
"""Solve a CPIT problem."""
if data.precedences is None:
raise ValueError("Precedences must be loaded")
ampl = AMPL()
ampl.eval(CPIT_MODEL)
# Load scalar parameters
ampl.param["NBLOCKS"] = data.num_blocks
ampl.param["NPERIODS"] = data.num_periods
ampl.param["NRESOURCE_SIDE_CONSTRAINTS"] = data.num_resources
ampl.param["DISCOUNT_RATE"] = data.discount_rate
# Load profits
for b in range(data.num_blocks):
ampl.param["OBJECTIVE_FUNCTION"][b] = data.objective.get(b, 0.0)
# Load resource limits
for r in range(data.num_resources):
for t in range(data.num_periods):
lb = data.resource_limits.lower_bounds.get(r, {}).get(t, -math.inf)
ub = data.resource_limits.upper_bounds.get(r, {}).get(t, math.inf)
ampl.param["RESOURCE_CONSTRAINT_LB_LIMITS"][r, t] = _format_bound(lb)
ampl.param["RESOURCE_CONSTRAINT_UB_LIMITS"][r, t] = _format_bound(ub)
# Load resource coefficients
for b, res_dict in data.resource_coefficients.items():
for r, coef in res_dict.items():
ampl.param["RESOURCE_CONSTRAINT_COEFFICIENTS"][b, r] = coef
# Load precedences
for b in range(data.num_blocks):
ampl.set["PREC"][b] = data.precedences.get_predecessors(b)
# Configure and solve
ampl.option["solver"] = solver
if solver_options:
ampl.option[f"{solver}_options"] = solver_options
if verbose:
print(f"Solving CPIT: {data.name}")
print(
f" Blocks: {data.num_blocks}, Periods: {data.num_periods}, Resources: {data.num_resources}"
)
ampl.solve()
# Extract solution
status = ampl.get_value("solve_result")
solve_time = ampl.get_value("_solve_time")
objective = ampl.obj["Profit"].value()
extraction_period = {}
X = ampl.var["X"]
for b in range(data.num_blocks):
extraction_period[b] = None
for t in range(data.num_periods):
if X[b, t].value() > 0.5:
extraction_period[b] = t
break
if verbose:
print(f" Status: {status}, NPV: ${objective:,.2f}, Time: {solve_time:.2f}s")
return SolutionCPIT(
objective=objective,
solve_time=solve_time,
status=status,
extraction_period=extraction_period,
)
def solve_pcpsp(
data: PCPSPData,
solver: str = "gurobi",
solver_options: str = "",
verbose: bool = True,
) -> SolutionPCPSP:
"""Solve a PCPSP problem."""
if data.precedences is None:
raise ValueError("Precedences must be loaded")
ampl = AMPL()
ampl.eval(PCPSP_MODEL)
# Load scalar parameters
ampl.param["NBLOCKS"] = data.num_blocks
ampl.param["NPERIODS"] = data.num_periods
ampl.param["NDESTINATIONS"] = data.num_destinations
ampl.param["NRESOURCE_SIDE_CONSTRAINTS"] = data.num_resources
ampl.param["DISCOUNT_RATE"] = data.discount_rate
# Load profits (per block and destination)
for b in range(data.num_blocks):
for d in range(data.num_destinations):
val = data.objective.get(b, {}).get(d, 0.0)
ampl.param["OBJECTIVE_FUNCTION"][b, d] = val
# Load resource limits
for r in range(data.num_resources):
for t in range(data.num_periods):
lb = data.resource_limits.lower_bounds.get(r, {}).get(t, -math.inf)
ub = data.resource_limits.upper_bounds.get(r, {}).get(t, math.inf)
ampl.param["RESOURCE_CONSTRAINT_LB_LIMITS"][r, t] = _format_bound(lb)
ampl.param["RESOURCE_CONSTRAINT_UB_LIMITS"][r, t] = _format_bound(ub)
# Load resource coefficients (per block, resource, destination)
for b, dest_dict in data.resource_coefficients.items():
for d, res_dict in dest_dict.items():
for r, coef in res_dict.items():
ampl.param["RESOURCE_CONSTRAINT_COEFFICIENTS"][b, r, d] = coef
# Load precedences
for b in range(data.num_blocks):
ampl.set["PREC"][b] = data.precedences.get_predecessors(b)
# Configure and solve
ampl.option["solver"] = solver
if solver_options:
ampl.option[f"{solver}_options"] = solver_options
if verbose:
print(f"Solving PCPSP: {data.name}")
print(
f" Blocks: {data.num_blocks}, Periods: {data.num_periods}, Destinations: {data.num_destinations}"
)
ampl.solve()
# Extract solution
status = ampl.get_value("solve_result")
solve_time = ampl.get_value("_solve_time")
objective = ampl.obj["Profit"].value()
extraction_period = {}
destination_allocation = {}
X = ampl.var["X"]
Y = ampl.var["Y"]
for b in range(data.num_blocks):
extraction_period[b] = None
for t in range(data.num_periods):
if X[b, t].value() > 0.5:
extraction_period[b] = t
destination_allocation[b] = {}
for d in range(data.num_destinations):
y_val = Y[b, d, t].value()
if y_val > 1e-6:
destination_allocation[b][d] = y_val
break
if verbose:
print(f" Status: {status}, NPV: ${objective:,.2f}, Time: {solve_time:.2f}s")
return SolutionPCPSP(
objective=objective,
solve_time=solve_time,
status=status,
extraction_period=extraction_period,
destination_allocation=destination_allocation,
)
# Base URL for raw GitHub content
BASE_URL = "https://raw.githubusercontent.com/ampl/colab.ampl.com/master/authors/eduardosalaz/minelib/data"
def fetch_minelib_file(instance: str, extension: str) -> StringIO:
"""
Fetch a MineLib file from GitHub and return as file-like object.
Args:
instance: Instance name (e.g., 'newman1', 'zuck_small')
extension: File extension (e.g., 'prec', 'upit', 'cpit', 'pcpsp', 'blocks')
Returns:
StringIO object that can be used like a file
"""
url = f"{BASE_URL}/{instance}/{instance}.{extension}"
response = requests.get(url)
response.raise_for_status() # Raise exception if download fails
return StringIO(response.text)
# Example usage
INSTANCE = "newman1"
# Fetch files as file-like objects
prec_file = fetch_minelib_file(INSTANCE, "prec")
upit_file = fetch_minelib_file(INSTANCE, "upit")
cpit_file = fetch_minelib_file(INSTANCE, "cpit")
pcpsp_file = fetch_minelib_file(INSTANCE, "pcpsp")
precedences = parse_precedences(prec_file)
# Parse each problem type and assign precedences
upit_data = parse_upit(upit_file)
upit_data.precedences = precedences
cpit_data = parse_cpit(cpit_file)
cpit_data.precedences = precedences
pcpsp_data = parse_pcpsp(pcpsp_file)
pcpsp_data.precedences = precedences
print(f"Instance: {upit_data.name}")
print(f"Blocks: {upit_data.num_blocks}")
print(f"Precedence arcs: {upit_data.precedences.num_precedence_arcs()}")
# Count ore vs waste blocks
positive = sum(1 for v in upit_data.objective.values() if v > 0)
negative = sum(1 for v in upit_data.objective.values() if v < 0)
print(f"Ore blocks (positive value): {positive}")
print(f"Waste blocks (negative value): {negative}")
# %%
# Solve UPIT
upit_solution = solve_upit(upit_data, solver="gurobi", verbose=True)
print(
f"\nBlocks in optimal pit: {upit_solution.num_extracted()} / {upit_data.num_blocks}"
)
Instance: Newman1
Blocks: 1060
Precedence arcs: 3922
Ore blocks (positive value): 545
Waste blocks (negative value): 515
Solving UPIT: Newman1
Blocks: 1060, Arcs: 3922
Gurobi 13.0.0:Gurobi 13.0.0: optimal solution; objective 26086899.03
0 simplex iterations
Status: solved, Objective: $26,086,899.03, Time: 0.07s
Blocks in optimal pit: 1059 / 1060
print(f"Instance: {cpit_data.name}")
print(f"Blocks: {cpit_data.num_blocks}, Periods: {cpit_data.num_periods}")
print(f"Resources: {cpit_data.num_resources}, Discount rate: {cpit_data.discount_rate}")
# %%
# Solve CPIT (may take longer due to time-indexed variables)
cpit_solution = solve_cpit(
cpit_data,
solver="gurobi",
solver_options="timelimit=300 outlev=1 MIPGap=1e-9",
verbose=True,
)
print(f"\nBlocks extracted: {cpit_solution.num_extracted()} / {cpit_data.num_blocks}")
# %%
# Analyze extraction schedule
print("\nExtraction Schedule:")
print(f"{'Period':<10} {'Blocks':<15} {'Cumulative':<15}")
print("-" * 40)
cumulative = 0
for t in range(cpit_data.num_periods):
count = len(cpit_solution.blocks_in_period(t))
cumulative += count
if count > 0:
print(f"{t:<10} {count:<15} {cumulative:<15}")
Instance: Newman1
Blocks: 1060, Periods: 6
Resources: 2, Discount rate: 0.08
Solving CPIT: Newman1
Blocks: 1060, Periods: 6, Resources: 2
Gurobi 13.0.0: lim:time = 300
Set parameter LogToConsole to value 1
tech:outlev = 1
Set parameter MIPGap to value 1e-09
mip:gap = 1.0000000000000001e-09
AMPL MP initial flat model has 6360 variables (0 integer, 6360 binary);
Objectives: 1 linear;
Constraints: 24604 linear;
AMPL MP final model has 6360 variables (0 integer, 6360 binary);
Objectives: 1 linear;
Constraints: 24604 linear;
Set parameter InfUnbdInfo to value 1
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (linux64 - "Ubuntu 22.04.5 LTS")
CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Non-default parameters:
TimeLimit 300
MIPGap 1e-09
InfUnbdInfo 1
Optimize a model with 24604 rows, 6360 columns and 180876 nonzeros (Max)
Model fingerprint: 0xc6f724c7
Model has 6360 linear objective coefficients
Variable types: 0 continuous, 6360 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 6e+03]
Objective range [5e+00, 2e+05]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+06]
Found heuristic solution: objective 989307.55494
Found heuristic solution: objective 1822786.3759
Presolve removed 605 rows and 147 columns
Presolve time: 0.87s
Presolved: 23999 rows, 6213 columns, 173942 nonzeros
Variable types: 0 continuous, 6213 integer (6213 binary)
Found heuristic solution: objective 2014319.6047
Root relaxation: objective 2.441235e+07, 7877 iterations, 2.47 seconds (1.98 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 2.4412e+07 0 2270 2014319.60 2.4412e+07 1112% - 4s
0 0 2.4401e+07 0 2252 2014319.60 2.4401e+07 1111% - 4s
0 0 2.4396e+07 0 2199 2014319.60 2.4396e+07 1111% - 5s
H 0 0 5009329.3557 2.4396e+07 387% - 5s
H 0 0 5083504.7013 2.4396e+07 380% - 5s
H 0 0 1.016018e+07 2.4396e+07 140% - 5s
H 0 0 1.223241e+07 2.4396e+07 99.4% - 5s
H 0 0 1.223841e+07 2.4396e+07 99.3% - 5s
H 0 0 1.259800e+07 2.4396e+07 93.7% - 5s
H 0 0 1.262926e+07 2.4396e+07 93.2% - 5s
H 0 0 1.830650e+07 2.4396e+07 33.3% - 5s
H 0 0 2.061441e+07 2.4396e+07 18.3% - 5s
H 0 0 2.070077e+07 2.4396e+07 17.9% - 5s
H 0 0 2.073573e+07 2.4396e+07 17.7% - 5s
H 0 0 2.073648e+07 2.4396e+07 17.6% - 5s
H 0 0 2.219034e+07 2.4396e+07 9.94% - 5s
H 0 0 2.257133e+07 2.4396e+07 8.08% - 5s
0 0 2.4378e+07 0 2162 2.2571e+07 2.4378e+07 8.00% - 8s
H 0 0 2.270092e+07 2.4378e+07 7.39% - 8s
H 0 0 2.279696e+07 2.4378e+07 6.93% - 8s
H 0 0 2.280138e+07 2.4378e+07 6.91% - 8s
0 0 2.4378e+07 0 2162 2.2801e+07 2.4378e+07 6.91% - 11s
H 0 0 2.336037e+07 2.4378e+07 4.35% - 12s
H 0 0 2.402383e+07 2.4378e+07 1.47% - 12s
H 0 0 2.402575e+07 2.4378e+07 1.46% - 12s
H 0 0 2.402845e+07 2.4378e+07 1.45% - 12s
H 0 0 2.402886e+07 2.4378e+07 1.45% - 12s
H 0 0 2.403116e+07 2.4378e+07 1.44% - 13s
H 0 0 2.403198e+07 2.4378e+07 1.44% - 13s
0 2 2.4378e+07 0 2162 2.4032e+07 2.4378e+07 1.44% - 13s
H 3 3 2.408751e+07 2.4378e+07 1.20% 331 15s
H 3 3 2.409004e+07 2.4378e+07 1.19% 331 15s
H 3 3 2.409669e+07 2.4378e+07 1.17% 331 15s
H 3 3 2.409861e+07 2.4378e+07 1.16% 331 15s
H 3 3 2.410128e+07 2.4378e+07 1.15% 331 15s
H 3 3 2.411263e+07 2.4378e+07 1.10% 331 15s
H 4 2 2.411452e+07 2.4378e+07 1.09% 436 15s
H 4 2 2.413434e+07 2.4378e+07 1.01% 436 15s
H 4 2 2.413533e+07 2.4378e+07 1.00% 436 15s
H 4 2 2.413660e+07 2.4378e+07 1.00% 436 15s
H 8 4 2.414115e+07 2.4378e+07 0.98% 304 16s
H 8 4 2.414283e+07 2.4378e+07 0.97% 304 16s
H 8 4 2.414450e+07 2.4378e+07 0.97% 304 16s
H 10 6 2.414599e+07 2.4378e+07 0.96% 319 16s
H 10 6 2.415008e+07 2.4378e+07 0.94% 319 16s
H 10 6 2.416976e+07 2.4378e+07 0.86% 319 17s
H 10 6 2.417021e+07 2.4378e+07 0.86% 319 17s
H 10 6 2.417391e+07 2.4378e+07 0.84% 319 17s
H 10 6 2.417654e+07 2.4378e+07 0.83% 319 17s
37 9 2.4274e+07 5 1016 2.4177e+07 2.4286e+07 0.45% 196 20s
H 84 26 2.417686e+07 2.4276e+07 0.41% 133 23s
125 29 cutoff 8 2.4177e+07 2.4208e+07 0.13% 119 25s
H 209 16 2.417686e+07 2.4178e+07 0.00% 84.9 26s
Cutting planes:
Gomory: 1
Cover: 7
Explored 323 nodes (30025 simplex iterations) in 27.72 seconds (19.88 work units)
Thread count was 2 (of 2 available processors)
Solution count 10: 2.41769e+07 2.41769e+07 2.41765e+07 ... 2.41428e+07
Optimal solution found (tolerance 1.00e-09)
Best objective 2.417686482482e+07, best bound 2.417686482482e+07, gap 0.0000%
Gurobi 13.0.0: optimal solution; objective 24176864.82
30025 simplex iterations
323 branching nodes
Status: solved, NPV: $24,176,864.82, Time: 27.64s
Blocks extracted: 1059 / 1060
Extraction Schedule:
Period Blocks Cumulative
----------------------------------------
0 407 407
1 365 772
2 287 1059
print(f"Instance: {pcpsp_data.name}")
print(f"Blocks: {pcpsp_data.num_blocks}, Destinations: {pcpsp_data.num_destinations}")
# %%
# Solve PCPSP
pcpsp_solution = solve_pcpsp(
pcpsp_data,
solver="gurobi",
solver_options="timelimit=600 outlev=1 mipgap=1e-9",
verbose=True,
)
print(f"\nBlocks extracted: {pcpsp_solution.num_extracted()}")
# Destination analysis
dest_names = {0: "Processing", 1: "Waste"}
dest_counts = {d: 0 for d in range(pcpsp_data.num_destinations)}
for b, alloc in pcpsp_solution.destination_allocation.items():
for d, frac in alloc.items():
if frac > 0.5:
dest_counts[d] += 1
print("\nDestination breakdown:")
for d, count in dest_counts.items():
name = dest_names.get(d, f"Dest {d}")
print(f" {name}: {count} blocks")
Instance: Newman1
Blocks: 1060, Destinations: 2
Solving PCPSP: Newman1
Blocks: 1060, Periods: 6, Destinations: 2
Gurobi 13.0.0: lim:time = 600
Set parameter LogToConsole to value 1
tech:outlev = 1
Set parameter MIPGap to value 1e-09
mip:gap = 1.0000000000000001e-09
AMPL MP initial flat model has 19080 variables (0 integer, 6360 binary);
Objectives: 1 linear;
Constraints: 30964 linear;
AMPL MP final model has 19080 variables (0 integer, 6360 binary);
Objectives: 1 linear;
Constraints: 30964 linear;
Set parameter InfUnbdInfo to value 1
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (linux64 - "Ubuntu 22.04.5 LTS")
CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Non-default parameters:
TimeLimit 600
MIPGap 1e-09
InfUnbdInfo 1
Optimize a model with 30964 rows, 19080 columns and 209244 nonzeros (Max)
Model fingerprint: 0x7f65152b
Model has 12720 linear objective coefficients
Variable types: 12720 continuous, 6360 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 6e+03]
Objective range [5e+00, 6e+06]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+06]
Found heuristic solution: objective -0.0000000
Presolve removed 3583 rows and 6222 columns
Presolve time: 1.78s
Presolved: 27381 rows, 12858 columns, 187916 nonzeros
Variable types: 6630 continuous, 6228 integer (6228 binary)
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
10783 5.5903458e+07 2.137874e+06 0.000000e+00 5s
20286 2.4538405e+07 1.058044e+05 0.000000e+00 10s
20939 2.4418598e+07 0.000000e+00 0.000000e+00 11s
20939 2.4418598e+07 0.000000e+00 0.000000e+00 11s
Root relaxation: objective 2.441860e+07, 20939 iterations, 8.71 seconds (6.91 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 2.4419e+07 0 2297 -0.00000 2.4419e+07 - - 11s
H 0 0 1.876798e+07 2.4419e+07 30.1% - 13s
H 0 0 1.877627e+07 2.4419e+07 30.1% - 13s
H 0 0 1.948454e+07 2.4419e+07 25.3% - 13s
H 0 0 1.948538e+07 2.4419e+07 25.3% - 13s
0 0 2.4410e+07 0 2271 1.9485e+07 2.4410e+07 25.3% - 16s
0 0 2.4392e+07 0 2236 1.9485e+07 2.4392e+07 25.2% - 19s
H 0 0 2.032787e+07 2.4392e+07 20.0% - 19s
H 0 0 2.034655e+07 2.4392e+07 19.9% - 19s
H 0 0 2.034673e+07 2.4392e+07 19.9% - 19s
H 0 0 2.034837e+07 2.4392e+07 19.9% - 19s
H 0 0 2.035169e+07 2.4392e+07 19.9% - 19s
H 0 0 2.035422e+07 2.4392e+07 19.8% - 19s
0 0 2.4392e+07 0 2236 2.0354e+07 2.4392e+07 19.8% - 21s
0 0 2.4392e+07 0 2236 2.0354e+07 2.4392e+07 19.8% - 21s
H 0 0 2.199142e+07 2.4392e+07 10.9% - 23s
H 0 0 2.199838e+07 2.4392e+07 10.9% - 23s
H 0 0 2.219597e+07 2.4392e+07 9.89% - 24s
H 0 0 2.324899e+07 2.4392e+07 4.92% - 24s
H 0 0 2.325183e+07 2.4392e+07 4.90% - 24s
H 0 0 2.366827e+07 2.4392e+07 3.06% - 25s
H 0 0 2.367461e+07 2.4392e+07 3.03% - 25s
H 0 0 2.367470e+07 2.4392e+07 3.03% - 25s
H 0 0 2.367506e+07 2.4392e+07 3.03% - 25s
H 0 0 2.386873e+07 2.4392e+07 2.19% - 25s
H 0 0 2.387615e+07 2.4392e+07 2.16% - 25s
H 0 0 2.387619e+07 2.4392e+07 2.16% - 25s
H 0 0 2.388455e+07 2.4392e+07 2.12% - 26s
H 0 0 2.388482e+07 2.4392e+07 2.12% - 26s
H 0 0 2.391956e+07 2.4392e+07 1.97% - 26s
0 2 2.4392e+07 0 2236 2.3920e+07 2.4392e+07 1.97% - 27s
H 2 2 2.392245e+07 2.4391e+07 1.96% 723 29s
3 5 2.4336e+07 2 1709 2.3922e+07 2.4390e+07 1.95% 658 30s
H 5 5 2.393256e+07 2.4336e+07 1.69% 477 32s
H 5 5 2.394418e+07 2.4336e+07 1.64% 477 32s
H 5 5 2.394777e+07 2.4336e+07 1.62% 477 32s
H 5 5 2.395983e+07 2.4336e+07 1.57% 477 32s
H 8 8 2.397002e+07 2.4299e+07 1.37% 460 33s
H 8 8 2.397093e+07 2.4299e+07 1.37% 460 33s
H 8 8 2.397093e+07 2.4299e+07 1.37% 460 33s
H 8 8 2.398968e+07 2.4299e+07 1.29% 460 33s
H 8 8 2.399236e+07 2.4299e+07 1.28% 460 33s
H 8 8 2.399341e+07 2.4299e+07 1.27% 460 33s
H 8 8 2.399426e+07 2.4299e+07 1.27% 460 33s
H 11 11 2.401220e+07 2.4299e+07 1.19% 433 34s
H 11 11 2.401416e+07 2.4299e+07 1.19% 433 34s
H 11 11 2.401974e+07 2.4299e+07 1.16% 433 34s
H 11 11 2.402042e+07 2.4299e+07 1.16% 433 34s
12 14 2.4096e+07 6 554 2.4020e+07 2.4299e+07 1.16% 445 35s
H 14 14 2.410346e+07 2.4299e+07 0.81% 407 35s
H 14 14 2.410476e+07 2.4299e+07 0.81% 407 35s
H 14 14 2.411542e+07 2.4299e+07 0.76% 407 35s
H 21 17 2.412290e+07 2.4299e+07 0.73% 328 36s
H 21 17 2.413239e+07 2.4299e+07 0.69% 328 36s
H 21 17 2.413241e+07 2.4299e+07 0.69% 328 37s
H 65 42 2.413263e+07 2.4299e+07 0.69% 132 40s
H 65 26 2.413419e+07 2.4299e+07 0.68% 132 40s
H 66 27 2.413443e+07 2.4297e+07 0.68% 144 40s
H 66 27 2.413454e+07 2.4297e+07 0.68% 144 40s
H 66 14 2.414473e+07 2.4297e+07 0.63% 144 40s
H 146 74 2.414691e+07 2.4219e+07 0.30% 98.1 43s
H 146 74 2.415052e+07 2.4219e+07 0.28% 98.1 43s
H 146 74 2.415054e+07 2.4219e+07 0.28% 98.1 43s
H 146 62 2.415406e+07 2.4219e+07 0.27% 98.1 43s
166 70 2.4178e+07 13 14 2.4154e+07 2.4200e+07 0.19% 95.5 45s
H 189 36 2.416159e+07 2.4200e+07 0.16% 87.3 45s
H 189 36 2.416533e+07 2.4200e+07 0.14% 87.3 45s
H 189 36 2.417362e+07 2.4200e+07 0.11% 87.3 45s
H 189 36 2.417563e+07 2.4200e+07 0.10% 87.3 45s
* 202 40 40 2.417570e+07 2.4200e+07 0.10% 84.0 45s
* 247 49 27 2.417576e+07 2.4178e+07 0.01% 73.7 46s
* 333 72 32 2.417644e+07 2.4178e+07 0.01% 61.8 48s
* 367 65 33 2.417686e+07 2.4178e+07 0.00% 57.8 48s
417 72 2.4177e+07 25 14 2.4177e+07 2.4178e+07 0.00% 54.3 50s
* 610 102 28 2.417686e+07 2.4177e+07 0.00% 45.7 54s
657 100 2.4177e+07 25 14 2.4177e+07 2.4177e+07 0.00% 45.9 55s
852 32 2.4177e+07 17 6 2.4177e+07 2.4177e+07 0.00% 43.7 60s
Cutting planes:
Gomory: 2
Explored 920 nodes (59702 simplex iterations) in 61.57 seconds (45.56 work units)
Thread count was 2 (of 2 available processors)
Solution count 10: 2.41769e+07 2.41769e+07 2.41764e+07 ... 2.41541e+07
Optimal solution found (tolerance 1.00e-09)
Best objective 2.417686482482e+07, best bound 2.417686482482e+07, gap 0.0000%
Gurobi 13.0.0: optimal solution; objective 24176864.82
59702 simplex iterations
920 branching nodes
Status: solved, NPV: $24,176,864.82, Time: 61.28s
Blocks extracted: 1059
Destination breakdown:
Processing: 487 blocks
Waste: 572 blocks
print("=" * 65)
print("Comparison of Formulations")
print("=" * 65)
print(f"{'Model':<10} {'Objective ($)':<20} {'Blocks':<12} {'Time (s)':<10}")
print("-" * 65)
print(
f"{'UPIT':<10} {upit_solution.objective:>18,.0f} {upit_solution.num_extracted():>10} {upit_solution.solve_time:>9.2f}"
)
print(
f"{'CPIT':<10} {cpit_solution.objective:>18,.0f} {cpit_solution.num_extracted():>10} {cpit_solution.solve_time:>9.2f}"
)
print(
f"{'PCPSP':<10} {pcpsp_solution.objective:>18,.0f} {pcpsp_solution.num_extracted():>10} {pcpsp_solution.solve_time:>9.2f}"
)
print("-" * 65)
# NPV reduction due to discounting and constraints
reduction = (
(upit_solution.objective - cpit_solution.objective) / upit_solution.objective * 100
)
print(f"\nNPV reduction (UPIT -> CPIT): {reduction:.1f}%")
=================================================================
Comparison of Formulations
=================================================================
Model Objective ($) Blocks Time (s)
-----------------------------------------------------------------
UPIT 26,086,899 1059 0.07
CPIT 24,176,865 1059 27.64
PCPSP 24,176,865 1059 61.28
-----------------------------------------------------------------
NPV reduction (UPIT -> CPIT): 7.3%
7. Conclusions#
This notebook demonstrated solving MineLib open pit mining problems with AMPL:
UPIT: Polynomial-time pit boundary optimization (network flow)
CPIT: NP-hard scheduling with time periods and resource constraints
PCPSP: Variable cutoff grades with destination optimization
Key observations:
Strong precedence formulation improves LP bounds
Discounting and constraints reduce achievable NPV vs. UPIT
Variable cutoff (PCPSP) can recover some value vs. fixed cutoff (CPIT)