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()

# add a model objective
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: 
.ix is deprecated. Please use
.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:
    model.constraints.add(0 <= model.y[p] <= max_production[p])
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:
    model.constraints.add(0 <= model.y[p] <= max_production[p])
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()

# for access to dual solution for constraints
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