Economic Planning¶
Reference¶
SAS/OR example: http://go.documentation.sas.com/?docsetId=ormpex&docsetTarget=ormpex_ex9_toc.htm&docsetVersion=15.1&locale=en
SAS/OR code for example: http://support.sas.com/documentation/onlinedoc/or/ex_code/151/mpex09.html
Model¶
import sasoptpy as so
import pandas as pd
def test(cas_conn):
m = so.Model(name='economic_planning', session=cas_conn)
industry_data = pd.DataFrame([
['coal', 150, 300, 60],
['steel', 80, 350, 60],
['transport', 100, 280, 30]
], columns=['industry', 'init_stocks', 'init_productive_capacity',
'demand']).set_index(['industry'])
production_data = pd.DataFrame([
['coal', 0.1, 0.5, 0.4],
['steel', 0.1, 0.1, 0.2],
['transport', 0.2, 0.1, 0.2],
['manpower', 0.6, 0.3, 0.2],
], columns=['input', 'coal',
'steel', 'transport']).set_index(['input'])
productive_capacity_data = pd.DataFrame([
['coal', 0.0, 0.7, 0.9],
['steel', 0.1, 0.1, 0.2],
['transport', 0.2, 0.1, 0.2],
['manpower', 0.4, 0.2, 0.1],
], columns=['input', 'coal',
'steel', 'transport']).set_index(['input'])
manpower_capacity = 470
num_years = 5
YEARS = list(range(1, num_years+1))
YEARS0 = [0] + list(YEARS)
INDUSTRIES = industry_data.index.tolist()
init_stocks = industry_data['init_stocks']
init_productive_capacity = industry_data['init_productive_capacity']
demand = industry_data['demand']
production_coeff = so.flatten_frame(production_data)
productive_capacity_coeff = so.flatten_frame(productive_capacity_data)
static_production = m.add_variables(INDUSTRIES, lb=0,
name='static_production')
m.set_objective(0, sense=so.MIN, name='Zero')
m.add_constraints((static_production[i] == demand[i] +
so.expr_sum(
production_coeff[i, j] * static_production[j]
for j in INDUSTRIES) for i in INDUSTRIES),
name='static_con')
m.solve()
print(so.get_value_table(static_production))
final_demand = so.get_value_table(
static_production)['static_production']
production = m.add_variables(INDUSTRIES, range(0, num_years+2), lb=0,
name='production')
stock = m.add_variables(INDUSTRIES, range(0, num_years+2), lb=0,
name='stock')
extra_capacity = m.add_variables(INDUSTRIES, range(2, num_years+3), lb=0,
name='extra_capacity')
productive_capacity = so.ImplicitVar(
(init_productive_capacity[i] +
so.expr_sum(extra_capacity[i, y] for y in range(2, year+1))
for i in INDUSTRIES for year in range(1, num_years+2)),
name='productive_capacity'
)
for i in INDUSTRIES:
production[i, 0].set_bounds(ub=0)
stock[i, 0].set_bounds(lb=init_stocks[i], ub=init_stocks[i])
total_productive_capacity = sum(productive_capacity[i, num_years]
for i in INDUSTRIES)
total_production = so.expr_sum(production[i, year] for i in INDUSTRIES
for year in [4, 5])
total_manpower = so.expr_sum(production_coeff['manpower', i] *
production[i, year+1] +
productive_capacity_coeff['manpower', i] *
extra_capacity[i, year+2]
for i in INDUSTRIES for year in YEARS)
continuity_con = m.add_constraints((
stock[i, year] + production[i, year] ==
(demand[i] if year in YEARS else 0) +
so.expr_sum(production_coeff[i, j] * production[j, year+1] +
productive_capacity_coeff[i, j] *
extra_capacity[j, year+2] for j in INDUSTRIES) +
stock[i, year+1]
for i in INDUSTRIES for year in YEARS0), name='continuity_con')
manpower_con = m.add_constraints((
so.expr_sum(production_coeff['manpower', j] * production[j, year] +
productive_capacity_coeff['manpower', j] *
extra_capacity[j, year+1]
for j in INDUSTRIES)
<= manpower_capacity for year in range(1, num_years+2)),
name='manpower_con')
capacity_con = m.add_constraints((production[i, year] <=
productive_capacity[i, year]
for i in INDUSTRIES
for year in range(1, num_years+2)),
name='capacity_con')
for i in INDUSTRIES:
production[i, num_years+1].set_bounds(lb=final_demand[i])
for i in INDUSTRIES:
for year in [num_years+1, num_years+2]:
extra_capacity[i, year].set_bounds(ub=0)
problem1 = so.Model(name='Problem1', session=cas_conn)
problem1.include(
production, stock, extra_capacity, continuity_con, manpower_con,
capacity_con, productive_capacity)
problem1.set_objective(total_productive_capacity, sense=so.MAX,
name='total_productive_capacity')
problem1.solve()
so.pd.display_dense()
print(so.get_value_table(production, stock, extra_capacity,
productive_capacity).sort_index())
print(so.get_value_table(manpower_con.get_expressions()))
# Problem 2
problem2 = so.Model(name='Problem2', session=cas_conn)
problem2.include(problem1)
problem2.set_objective(total_production, name='total_production',
sense=so.MAX)
for i in INDUSTRIES:
for year in YEARS:
continuity_con[i, year].set_rhs(0)
problem2.solve()
print(so.get_value_table(production, stock, extra_capacity,
productive_capacity).sort_index())
print(so.get_value_table(manpower_con.get_expressions()))
# Problem 3
problem3 = so.Model(name='Problem3', session=cas_conn)
problem3.include(production, stock, extra_capacity, continuity_con,
capacity_con)
problem3.set_objective(total_manpower, sense=so.MAX, name='total_manpower')
for i in INDUSTRIES:
for year in YEARS:
continuity_con[i, year].set_rhs(demand[i])
problem3.solve()
print(so.get_value_table(production, stock, extra_capacity,
productive_capacity).sort_index())
print(so.get_value_table(manpower_con.get_expressions()))
return problem3.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.economic_planning import test
In [8]: test(cas_conn)
NOTE: Initialized model economic_planning.
NOTE: Added action set 'optimization'.
NOTE: Converting model economic_planning to OPTMODEL.
NOTE: Submitting OPTMODEL code to CAS server.
NOTE: Problem generation will use 8 threads.
NOTE: The problem has 3 variables (0 free, 0 fixed).
NOTE: The problem has 3 linear constraints (0 LE, 3 EQ, 0 GE, 0 range).
NOTE: The problem has 9 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 LP presolver value AUTOMATIC is applied.
NOTE: The LP presolver time is 0.00 seconds.
NOTE: The LP presolver removed all variables and constraints.
NOTE: Optimal.
NOTE: Objective = 0.
NOTE: The output table 'SOLUTION' in caslib 'CASUSER(casuser)' has 3 rows and 6 columns.
NOTE: The output table 'DUAL' in caslib 'CASUSER(casuser)' has 3 rows and 4 columns.
NOTE: The CAS table 'solutionSummary' in caslib 'CASUSER(casuser)' has 13 rows and 4 columns.
NOTE: The CAS table 'problemSummary' in caslib 'CASUSER(casuser)' has 18 rows and 4 columns.
static_production
coal 166.396761
steel 105.668016
transport 92.307692
NOTE: Initialized model Problem1.
NOTE: Added action set 'optimization'.
NOTE: Converting model Problem1 to OPTMODEL.
NOTE: Submitting OPTMODEL code to CAS server.
NOTE: Problem generation will use 8 threads.
NOTE: The problem has 60 variables (0 free, 12 fixed).
NOTE: The problem has 42 linear constraints (24 LE, 18 EQ, 0 GE, 0 range).
NOTE: The problem has 255 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 LP presolver value AUTOMATIC is applied.
NOTE: The LP presolver time is 0.00 seconds.
NOTE: The LP presolver removed 15 variables and 3 constraints.
NOTE: The LP presolver removed 37 constraint coefficients.
NOTE: The presolved problem has 45 variables, 39 constraints, and 218 constraint coefficients.
NOTE: The LP solver is called.
NOTE: The Dual Simplex algorithm is used.
Objective
Phase Iteration Value Time
D 2 1 2.683246E+04 0
P 2 42 2.141875E+03 0
NOTE: Optimal.
NOTE: Objective = 2141.8751967.
NOTE: The Dual Simplex solve time is 0.02 seconds.
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 42 rows and 4 columns.
NOTE: The CAS table 'solutionSummary' in caslib 'CASUSER(casuser)' has 13 rows and 4 columns.
NOTE: The CAS table 'problemSummary' in caslib 'CASUSER(casuser)' has 18 rows and 4 columns.
production stock extra_capacity productive_capacity
coal 0 0.000000 1.500000e+02 NaN NaN
coal 1 260.402615 0.000000e+00 NaN 300.000000
coal 2 293.406208 0.000000e+00 0.000000 300.000000
coal 3 300.000000 0.000000e+00 0.000000 300.000000
coal 4 17.948718 1.484480e+02 189.203132 489.203132
coal 5 166.396761 0.000000e+00 1022.672065 1511.875197
coal 6 166.396761 1.421085e-14 0.000000 1511.875197
coal 7 NaN NaN 0.000000 NaN
steel 0 0.000000 8.000000e+01 NaN NaN
steel 1 135.341540 1.228110e+01 NaN 350.000000
steel 2 181.659854 0.000000e+00 0.000000 350.000000
steel 3 193.090418 0.000000e+00 0.000000 350.000000
steel 4 105.668016 0.000000e+00 0.000000 350.000000
steel 5 105.668016 0.000000e+00 0.000000 350.000000
steel 6 105.668016 -1.456613e-13 0.000000 350.000000
steel 7 NaN NaN 0.000000 NaN
transport 0 0.000000 1.000000e+02 NaN NaN
transport 1 140.722422 6.240839e+00 NaN 280.000000
transport 2 200.580168 0.000000e+00 0.000000 280.000000
transport 3 267.152497 0.000000e+00 0.000000 280.000000
transport 4 92.307692 0.000000e+00 0.000000 280.000000
transport 5 92.307692 0.000000e+00 0.000000 280.000000
transport 6 92.307692 -3.907985e-14 0.000000 280.000000
transport 7 NaN NaN 0.000000 NaN
manpower_con
1 224.988515
2 270.657715
3 367.038878
4 470.000000
5 150.000000
6 150.000000
NOTE: Initialized model Problem2.
NOTE: Added action set 'optimization'.
NOTE: Converting model Problem2 to OPTMODEL.
NOTE: Submitting OPTMODEL code to CAS server.
NOTE: Problem generation will use 8 threads.
NOTE: The problem has 60 variables (0 free, 12 fixed).
NOTE: The problem has 42 linear constraints (24 LE, 18 EQ, 0 GE, 0 range).
NOTE: The problem has 255 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 LP presolver value AUTOMATIC is applied.
NOTE: The LP presolver time is 0.00 seconds.
NOTE: The LP presolver removed 15 variables and 3 constraints.
NOTE: The LP presolver removed 37 constraint coefficients.
NOTE: The presolved problem has 45 variables, 39 constraints, and 218 constraint coefficients.
NOTE: The LP solver is called.
NOTE: The Dual Simplex algorithm is used.
Objective
Phase Iteration Value Time
D 2 1 1.504360E+04 0
P 2 46 2.618579E+03 0
NOTE: Optimal.
NOTE: Objective = 2618.5791147.
NOTE: The Dual Simplex solve time is 0.01 seconds.
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 42 rows and 4 columns.
NOTE: The CAS table 'solutionSummary' in caslib 'CASUSER(casuser)' has 13 rows and 4 columns.
NOTE: The CAS table 'problemSummary' in caslib 'CASUSER(casuser)' has 18 rows and 4 columns.
production stock extra_capacity productive_capacity
coal 0 0.000000 150.000000 NaN NaN
coal 1 184.818327 31.628509 NaN 300.000000
coal 2 430.504654 16.372454 1.305047e+02 430.504654
coal 3 430.504654 0.000000 0.000000e+00 430.504654
coal 4 430.504654 0.000000 -1.207321e-12 430.504654
coal 5 430.504654 0.000000 0.000000e+00 430.504654
coal 6 166.396761 324.107893 0.000000e+00 430.504654
coal 7 NaN NaN 0.000000e+00 NaN
steel 0 0.000000 80.000000 NaN NaN
steel 1 86.729504 11.532298 NaN 350.000000
steel 2 155.337478 0.000000 0.000000e+00 350.000000
steel 3 182.867219 0.000000 0.000000e+00 350.000000
steel 4 359.402270 0.000000 9.402270e+00 359.402270
steel 5 359.402270 176.535051 0.000000e+00 359.402270
steel 6 105.668016 490.269305 0.000000e+00 359.402270
steel 7 NaN NaN 0.000000e+00 NaN
transport 0 0.000000 100.000000 NaN NaN
transport 1 141.312267 0.000000 NaN 280.000000
transport 2 198.387943 0.000000 0.000000e+00 280.000000
transport 3 225.917684 0.000000 0.000000e+00 280.000000
transport 4 519.382633 0.000000 2.393826e+02 519.382633
transport 5 519.382633 293.464949 0.000000e+00 519.382633
transport 6 92.307692 750.539890 0.000000e+00 519.382633
transport 7 NaN NaN 0.000000e+00 NaN
manpower_con
1 217.374162
2 344.581624
3 384.165212
4 470.000000
5 470.000000
6 150.000000
NOTE: Initialized model Problem3.
NOTE: Added action set 'optimization'.
NOTE: Converting model Problem3 to OPTMODEL.
NOTE: Submitting OPTMODEL code to CAS server.
NOTE: Problem generation will use 8 threads.
NOTE: The problem has 60 variables (0 free, 12 fixed).
NOTE: The problem has 36 linear constraints (18 LE, 18 EQ, 0 GE, 0 range).
NOTE: The problem has 219 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 LP presolver value AUTOMATIC is applied.
NOTE: The LP presolver time is 0.00 seconds.
NOTE: The LP presolver removed 15 variables and 3 constraints.
NOTE: The LP presolver removed 31 constraint coefficients.
NOTE: The presolved problem has 45 variables, 33 constraints, and 188 constraint coefficients.
NOTE: The LP solver is called.
NOTE: The Dual Simplex algorithm is used.
Objective
Phase Iteration Value Time
D 2 1 1.464016E+05 0
D 2 54 2.450706E+03 0
P 2 56 2.450027E+03 0
NOTE: Optimal.
NOTE: Objective = 2450.0266228.
NOTE: The Dual Simplex solve time is 0.01 seconds.
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 36 rows and 4 columns.
NOTE: The CAS table 'solutionSummary' in caslib 'CASUSER(casuser)' has 13 rows and 4 columns.
NOTE: The CAS table 'problemSummary' in caslib 'CASUSER(casuser)' has 18 rows and 4 columns.
production stock extra_capacity productive_capacity
coal 0 0.000000 150.000000 NaN NaN
coal 1 251.792754 0.000000 NaN 300.000000
coal 2 316.015222 0.000000 16.015222 316.015222
coal 3 319.832020 0.000000 3.816798 319.832020
coal 4 366.349753 0.000000 46.517734 366.349753
coal 5 859.359606 0.000000 493.009853 859.359606
coal 6 859.359606 460.207993 0.000000 859.359606
coal 7 NaN NaN 0.000000 NaN
steel 0 0.000000 80.000000 NaN NaN
steel 1 134.794583 11.028028 NaN 350.000000
steel 2 175.041379 0.000000 0.000000 350.000000
steel 3 224.064039 0.000000 0.000000 350.000000
steel 4 223.136289 0.000000 0.000000 350.000000
steel 5 220.043787 0.000000 0.000000 350.000000
steel 6 350.000000 0.000000 0.000000 350.000000
steel 7 NaN NaN 0.000000 NaN
transport 0 0.000000 100.000000 NaN NaN
transport 1 143.558583 4.247230 NaN 280.000000
transport 2 181.676355 0.000000 0.000000 280.000000
transport 3 280.000000 0.000000 0.000000 280.000000
transport 4 279.072249 0.000000 0.000000 280.000000
transport 5 275.979748 0.000000 0.000000 280.000000
transport 6 195.539132 0.000000 0.000000 280.000000
transport 7 NaN NaN 0.000000 NaN
manpower_con
1 226.631832
2 279.983537
3 333.725517
4 539.769130
5 636.824849
6 659.723590
Out[8]: 2450.026622821294