# Pyomo Examples¶

## Example: Nonlinear Optimization of Series Reaction in a Continuous Stirred Tank Reactor¶

In [14]:
from pyomo.environ import *

V = 40     # liters
kA = 0.5   # 1/min
kB = 0.1   # l/min
CAf = 2.0  # moles/liter

# create a model instance
model = ConcreteModel()

# create x and y variables in the model
model.q = Var()

model.objective = Objective(expr = model.q*V*kA*CAf/(model.q + V*kB)/(model.q + V*kA), sense=maximize)

# compute a solution using ipopt for nonlinear optimization
results = SolverFactory('ipopt').solve(model)
model.pprint()

# print solutions
qmax = model.q()
CBmax = model.objective()
print('\nFlowrate at maximum CB = ', qmax, 'liters per minute.')
print('\nMaximum CB =', CBmax, 'moles per liter.')
print('\nProductivity = ', qmax*CBmax, 'moles per minute.')

1 Var Declarations
q : Size=1, Index=None
Key  : Lower : Value             : Upper : Fixed : Stale : Domain
None :  None : 8.944271909985442 :  None : False : False :  Reals

1 Objective Declarations
objective : Size=1, Index=None, Active=True
Key  : Active : Sense    : Expression
None :   True : maximize : 40.0 * q / ( ( 4.0 + q ) * ( 20.0 + q ) )

2 Declarations: q objective

Flowrate at maximum CB =  8.944271909985442 liters per minute.

Maximum CB = 0.954915028125263 moles per liter.

Productivity =  8.541019662483748 moles per minute.


## Example 19.3: Linear Programming Refinery¶

In [51]:
import pandas as pd

PRODUCTS = ['Gasoline', 'Kerosine', 'Fuel Oil', 'Residual']
FEEDS = ['Crude 1', 'Crude 2']

products = pd.DataFrame(index=PRODUCTS)
products['Price'] = [72, 48, 42, 20]
products['Max Production'] = [24000, 2000, 6000, 100000]

crudes = pd.DataFrame(index=FEEDS)
crudes['Processing Cost'] = [1.00, 2.00]
crudes['Feed Costs'] = [48, 30]

yields = pd.DataFrame(index=PRODUCTS)
yields['Crude 1'] = [0.80, 0.05, 0.10, 0.05]
yields['Crude 2'] = [0.44, 0.10, 0.36, 0.10]

print('\n', products)
print('\n', crudes)
print('\n', yields)

           Price  Max Production
Gasoline     72           24000
Kerosine     48            2000
Fuel Oil     42            6000
Residual     20          100000

Processing Cost  Feed Costs
Crude 1              1.0          48
Crude 2              2.0          30

Crude 1  Crude 2
Gasoline     0.80     0.44
Kerosine     0.05     0.10
Fuel Oil     0.10     0.36
Residual     0.05     0.10

In [64]:
price = {}
price['Gasoline'] = 72
price['Kerosine'] = 48
price

Out[64]:
{'Gasoline': 72, 'Kerosine': 48}
In [56]:
products.loc['Gasoline','Price']

Out[56]:
72
In [54]:
# model formulation
model = ConcreteModel()

# variables
model.x = Var(FEEDS,    domain=NonNegativeReals)
model.y = Var(PRODUCTS, domain=NonNegativeReals)

# objective
income = sum(products.ix[p, 'Price'] * model.y[p] for p in PRODUCTS)
raw_materials_cost = sum(crudes.ix[f,'Feed Costs'] * model.x[f] for f in FEEDS)

/Users/jeff/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:10: DeprecationWarning:
.loc for label based indexing or
.iloc for positional indexing

See the documentation here:
http://pandas.pydata.org/pandas-docs/stable/indexing.html#ix-indexer-is-deprecated
# Remove the CWD from sys.path while we load stuff.

In [49]:
processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS)
profit = income - raw_materials_cost - processing_cost
model.objective = Objective(expr = profit, sense=maximize)

# constraints
model.constraints = ConstraintList()
for p in PRODUCTS:
for p in PRODUCTS:
model.constraints.add(model.y[p] == sum(model.x[f] * product_yield[(f,p)] for f in FEEDS))

solver = SolverFactory('glpk')
solver.solve(model)
model.pprint()

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/anaconda3/lib/python3.6/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
2441             try:
-> 2442                 return self._engine.get_loc(key)
2443             except KeyError:

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: ('Gasoline', 'Price')

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
<ipython-input-49-472329535d57> in <module>()
7
8 # objective
----> 9 income = sum(products[p, 'Price'] * model.y[p] for p in PRODUCTS)
10 raw_materials_cost = sum(feed_costs[f] * model.x[f] for f in FEEDS)
11 processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS)

<ipython-input-49-472329535d57> in <genexpr>(.0)
7
8 # objective
----> 9 income = sum(products[p, 'Price'] * model.y[p] for p in PRODUCTS)
10 raw_materials_cost = sum(feed_costs[f] * model.x[f] for f in FEEDS)
11 processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS)

~/anaconda3/lib/python3.6/site-packages/pandas/core/frame.py in __getitem__(self, key)
1962             return self._getitem_multilevel(key)
1963         else:
-> 1964             return self._getitem_column(key)
1965
1966     def _getitem_column(self, key):

~/anaconda3/lib/python3.6/site-packages/pandas/core/frame.py in _getitem_column(self, key)
1969         # get column
1970         if self.columns.is_unique:
-> 1971             return self._get_item_cache(key)
1972
1973         # duplicate columns & possible reduce dimensionality

~/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py in _get_item_cache(self, item)
1643         res = cache.get(item)
1644         if res is None:
-> 1645             values = self._data.get(item)
1646             res = self._box_item_values(item, values)
1647             cache[item] = res

~/anaconda3/lib/python3.6/site-packages/pandas/core/internals.py in get(self, item, fastpath)
3588
3589             if not isnull(item):
-> 3590                 loc = self.items.get_loc(item)
3591             else:
3592                 indexer = np.arange(len(self.items))[isnull(self.items)]

~/anaconda3/lib/python3.6/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
2442                 return self._engine.get_loc(key)
2443             except KeyError:
-> 2444                 return self._engine.get_loc(self._maybe_cast_indexer(key))
2445
2446         indexer = self.get_indexer([key], method=method, tolerance=tolerance)

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: ('Gasoline', 'Price')
In [5]:
from pyomo.environ import *
import numpy as np

# problem data
FEEDS = ['Crude #1', 'Crude #2']
PRODUCTS = ['Gasoline', 'Kerosine', 'Fuel Oil', 'Residual']

# feed costs
feed_costs = {'Crude #1': 48,
'Crude #2': 30}

# processing costs
processing_costs = {'Crude #1': 1.00,
'Crude #2': 2.00}

# yield data
product_yield =
product_yield = {('Crude #1', 'Gasoline'): 0.80,
('Crude #1', 'Kerosine'): 0.05,
('Crude #1', 'Fuel Oil'): 0.10,
('Crude #1', 'Residual'): 0.05,
('Crude #2', 'Gasoline'): 0.44,
('Crude #2', 'Kerosine'): 0.10,
('Crude #2', 'Fuel Oil'): 0.36,
('Crude #2', 'Residual'): 0.10}

# product sales prices
sales_price = {'Gasoline': 72,
'Kerosine': 48,
'Fuel Oil': 42,
'Residual': 20}

# production limits
max_production = {'Gasoline': 24000,
'Kerosine': 2000,
'Fuel Oil': 6000,
'Residual': 100000}

# model formulation
model = ConcreteModel()

# variables
model.x = Var(FEEDS, domain=NonNegativeReals)
model.y = Var(PRODUCTS, domain=NonNegativeReals)

# objective
income = sum(sales_price[p] * model.y[p] for p in PRODUCTS)
raw_materials_cost = sum(feed_costs[f] * model.x[f] for f in FEEDS)
processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS)

profit = income - raw_materials_cost - processing_cost
model.objective = Objective(expr = profit, sense=maximize)

# constraints
model.constraints = ConstraintList()
for p in PRODUCTS:
for p in PRODUCTS:
model.constraints.add(model.y[p] == sum(model.x[f] * product_yield[(f,p)] for f in FEEDS))

solver = SolverFactory('glpk')
solver.solve(model)
model.pprint()

3 Set Declarations
constraints_index : Dim=0, Dimen=1, Size=8, Domain=None, Ordered=False, Bounds=None
[1, 2, 3, 4, 5, 6, 7, 8]
x_index : Dim=0, Dimen=1, Size=2, Domain=None, Ordered=False, Bounds=None
['Crude #1', 'Crude #2']
y_index : Dim=0, Dimen=1, Size=4, Domain=None, Ordered=False, Bounds=None
['Fuel Oil', 'Gasoline', 'Kerosine', 'Residual']

2 Var Declarations
x : Size=2, Index=x_index
Key      : Lower : Value            : Upper : Fixed : Stale : Domain
Crude #1 :     0 : 26206.8965517241 :  None : False : False : NonNegativeReals
Crude #2 :     0 : 6896.55172413793 :  None : False : False : NonNegativeReals
y : Size=4, Index=y_index
Key      : Lower : Value            : Upper : Fixed : Stale : Domain
Fuel Oil :     0 : 5103.44827586207 :  None : False : False : NonNegativeReals
Gasoline :     0 :          24000.0 :  None : False : False : NonNegativeReals
Kerosine :     0 :           2000.0 :  None : False : False : NonNegativeReals
Residual :     0 :           2000.0 :  None : False : False : NonNegativeReals

1 Objective Declarations
objective : Size=1, Index=None, Active=True
Key  : Active : Sense    : Expression
None :   True : maximize : 72*y[Gasoline] + 48*y[Kerosine] + 42*y[Fuel Oil] + 20*y[Residual] - 48*x[Crude #1] - 30*x[Crude #2] - x[Crude #1] - 2.0*x[Crude #2]

1 Constraint Declarations
constraints : Size=8, Index=constraints_index, Active=True
Key : Lower : Body                                             : Upper    : Active
1 :   0.0 :                                      y[Gasoline] :  24000.0 :   True
2 :   0.0 :                                      y[Kerosine] :   2000.0 :   True
3 :   0.0 :                                      y[Fuel Oil] :   6000.0 :   True
4 :   0.0 :                                      y[Residual] : 100000.0 :   True
5 :   0.0 : y[Gasoline] - 0.8*x[Crude #1] - 0.44*x[Crude #2] :      0.0 :   True
6 :   0.0 : y[Kerosine] - 0.05*x[Crude #1] - 0.1*x[Crude #2] :      0.0 :   True
7 :   0.0 : y[Fuel Oil] - 0.1*x[Crude #1] - 0.36*x[Crude #2] :      0.0 :   True
8 :   0.0 : y[Residual] - 0.05*x[Crude #1] - 0.1*x[Crude #2] :      0.0 :   True

7 Declarations: x_index x y_index y objective constraints_index constraints

In [6]:
profit()

Out[6]:
573517.2413793121
In [7]:
income()

Out[7]:
2078344.8275862068
In [8]:
raw_materials_cost()

Out[8]:
1464827.5862068948
In [9]:
processing_cost()

Out[9]:
39999.99999999996

## Example: Making Change¶

One of the important modeling features of Pyomo is the ability to index variables and constraints. The

In [8]:
from pyomo.environ import *

def make_change(amount, coins):
model = ConcreteModel()
model.x = Var(coins.keys(), domain=NonNegativeIntegers)
model.total = Objective(expr = sum(model.x[c] for c in coins.keys()), sense=minimize)
model.amount = Constraint(expr = sum(model.x[c]*coins[c] for c in coins.keys()) == amount)
SolverFactory('glpk').solve(model)
return {c: int(model.x[c]()) for c in coins.keys()}

# problem data
coins = {
'penny': 1,
'nickel': 5,
'dime': 10,
'quarter': 25
}

make_change(1034, coins)

Out[8]:
{'dime': 0, 'nickel': 1, 'penny': 4, 'quarter': 41}

## Example: Linear Production Model with Constraints with Duals¶

In [3]:
from pyomo.environ import *
model = ConcreteModel()

model.dual = Suffix(direction=Suffix.IMPORT)

# declare decision variables
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)

# declare objective
model.profit = Objective(expr = 40*model.x + 30*model.y, sense = maximize)

# declare constraints
model.demand = Constraint(expr = model.x <= 40)
model.laborA = Constraint(expr = model.x + model.y <= 80)
model.laborB = Constraint(expr = 2*model.x + model.y <= 100)

# solve
SolverFactory('glpk').solve(model).write()

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 2600.0
Upper bound: 2600.0
Number of objectives: 1
Number of constraints: 4
Number of variables: 3
Number of nonzeros: 6
Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Error rc: 0
Time: 0.01401972770690918
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0

In [4]:
str = "   {0:7.2f} {1:7.2f} {2:7.2f} {3:7.2f}"

print("Constraint  value  lslack  uslack    dual")
for c in [model.demand, model.laborA, model.laborB]:
print(c, str.format(c(), c.lslack(), c.uslack(), model.dual[c]))

Constraint  value  lslack  uslack    dual
demand      20.00    -inf   20.00    0.00
laborA      80.00    -inf    0.00   20.00
laborB     100.00    -inf    0.00   10.00