Mining Optimization¶
Reference¶
SAS/OR example: http://go.documentation.sas.com/?docsetId=ormpex&docsetTarget=ormpex_ex7_toc.htm&docsetVersion=15.1&locale=en
SAS/OR code for example: http://support.sas.com/documentation/onlinedoc/or/ex_code/151/mpex07.html
Model¶
import sasoptpy as so
import pandas as pd
def test(cas_conn):
m = so.Model(name='mining_optimization', session=cas_conn)
mine_data = pd.DataFrame([
['mine1', 5, 2, 1.0],
['mine2', 4, 2.5, 0.7],
['mine3', 4, 1.3, 1.5],
['mine4', 5, 3, 0.5],
], columns=['mine', 'cost', 'extract_ub', 'quality']).\
set_index(['mine'])
year_data = pd.DataFrame([
[1, 0.9],
[2, 0.8],
[3, 1.2],
[4, 0.6],
[5, 1.0],
], columns=['year', 'quality_required']).set_index(['year'])
max_num_worked_per_year = 3
revenue_per_ton = 10
discount_rate = 0.10
MINES = mine_data.index.tolist()
cost = mine_data['cost']
extract_ub = mine_data['extract_ub']
quality = mine_data['quality']
YEARS = year_data.index.tolist()
quality_required = year_data['quality_required']
isOpen = m.add_variables(MINES, YEARS, vartype=so.BIN, name='isOpen')
isWorked = m.add_variables(MINES, YEARS, vartype=so.BIN, name='isWorked')
extract = m.add_variables(MINES, YEARS, lb=0, name='extract')
[extract[i, j].set_bounds(ub=extract_ub[i]) for i in MINES for j in YEARS]
extractedPerYear = {j: extract.sum('*', j) for j in YEARS}
discount = {j: 1 / (1+discount_rate) ** (j-1) for j in YEARS}
totalRevenue = revenue_per_ton *\
so.expr_sum(discount[j] * extractedPerYear[j] for j in YEARS)
totalCost = so.expr_sum(discount[j] * cost[i] * isOpen[i, j]
for i in MINES for j in YEARS)
m.set_objective(totalRevenue-totalCost, sense=so.MAX, name='totalProfit')
m.add_constraints((extract[i, j] <= extract[i, j]._ub * isWorked[i, j]
for i in MINES for j in YEARS), name='link')
m.add_constraints((isWorked.sum('*', j) <= max_num_worked_per_year
for j in YEARS), name='cardinality')
m.add_constraints((isWorked[i, j] <= isOpen[i, j] for i in MINES
for j in YEARS), name='worked_implies_open')
m.add_constraints((isOpen[i, j] <= isOpen[i, j-1] for i in MINES
for j in YEARS if j != 1), name='continuity')
m.add_constraints((so.expr_sum(quality[i] * extract[i, j] for i in MINES)
== quality_required[j] * extractedPerYear[j]
for j in YEARS), name='quality_con')
res = m.solve()
if res is not None:
print(so.get_solution_table(isOpen, isWorked, extract))
quality_sol = {j: so.expr_sum(quality[i] * extract[i, j].get_value()
for i in MINES)
/ extractedPerYear[j].get_value() for j in YEARS}
qs = so.dict_to_frame(quality_sol, ['quality_sol'])
epy = so.dict_to_frame(extractedPerYear, ['extracted_per_year'])
print(so.get_solution_table(epy, qs, quality_required))
return m.get_objective_value()
Output¶
In [1]: import os
In [2]: hostname = os.getenv('CASHOST')
In [3]: port = os.getenv('CASPORT')
In [4]: from swat import CAS
In [5]: cas_conn = CAS(hostname, port)
In [6]: import sasoptpy
In [7]: from examples.client_side.mining_optimization import test
In [8]: test(cas_conn)
NOTE: Initialized model mining_optimization.
NOTE: Added action set 'optimization'.
NOTE: Converting model mining_optimization to OPTMODEL.
NOTE: Submitting OPTMODEL code to CAS server.
NOTE: Problem generation will use 8 threads.
NOTE: The problem has 60 variables (0 free, 0 fixed).
NOTE: The problem has 40 binary and 0 integer variables.
NOTE: The problem has 66 linear constraints (61 LE, 5 EQ, 0 GE, 0 range).
NOTE: The problem has 151 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The OPTMODEL presolver is disabled for linear problems.
NOTE: The initial MILP heuristics are applied.
NOTE: The MILP presolver value AUTOMATIC is applied.
NOTE: The MILP presolver removed 8 variables and 8 constraints.
NOTE: The MILP presolver removed 16 constraint coefficients.
NOTE: The MILP presolver modified 8 constraint coefficients.
NOTE: The presolved problem has 52 variables, 58 constraints, and 135 constraint coefficients.
NOTE: The MILP solver is called.
NOTE: The parallel Branch and Cut algorithm is used.
NOTE: The Branch and Cut algorithm is using up to 8 threads.
Node Active Sols BestInteger BestBound Gap Time
0 1 14 96.5802313 157.7309278 38.77% 0
0 1 14 96.5802313 150.9548680 36.02% 0
0 1 14 96.5802313 147.3693449 34.46% 0
0 1 15 146.8619974 146.8619974 0.00% 0
0 0 15 146.8619974 146.8619974 0.00% 0
NOTE: The MILP solver added 7 cuts with 34 cut coefficients at the root.
NOTE: Optimal.
NOTE: Objective = 146.86199735.
NOTE: The output table 'SOLUTION' in caslib 'CASUSER(casuser)' has 60 rows and 6 columns.
NOTE: The output table 'DUAL' in caslib 'CASUSER(casuser)' has 66 rows and 4 columns.
NOTE: The CAS table 'solutionSummary' in caslib 'CASUSER(casuser)' has 18 rows and 4 columns.
NOTE: The CAS table 'problemSummary' in caslib 'CASUSER(casuser)' has 20 rows and 4 columns.
isOpen isWorked extract
(mine1, 1) 1.000000 1.000000 2.000000e+00
(mine1, 2) 0.999994 0.000006 1.252000e-05
(mine1, 3) 0.999994 0.999994 1.950000e+00
(mine1, 4) 0.999994 0.999994 1.250007e-01
(mine1, 5) 0.999994 0.999994 1.999987e+00
(mine2, 1) 1.000000 0.000000 0.000000e+00
(mine2, 2) 1.000000 0.999995 2.499988e+00
(mine2, 3) 0.999998 -0.000000 0.000000e+00
(mine2, 4) 0.999998 0.999998 2.499994e+00
(mine2, 5) 0.999998 0.999998 2.166667e+00
(mine3, 1) 1.000000 1.000000 1.300000e+00
(mine3, 2) 1.000000 0.999998 1.299997e+00
(mine3, 3) 1.000000 1.000000 1.300000e+00
(mine3, 4) 1.000000 0.000001 3.477765e-07
(mine3, 5) 1.000000 1.000000 1.300000e+00
(mine4, 1) 1.000000 1.000000 2.450000e+00
(mine4, 2) 1.000000 1.000000 2.200005e+00
(mine4, 3) 1.000000 -0.000000 0.000000e+00
(mine4, 4) 1.000000 1.000000 3.000000e+00
(mine4, 5) -0.000000 -0.000000 0.000000e+00
extracted_per_year quality_sol quality_required
1 5.750000 0.9 0.9
2 6.000003 0.8 0.8
3 3.250000 1.2 1.2
4 5.624995 0.6 0.6
5 5.466654 1.0 1.0
Out[8]: 146.861997350257