This notebook contains material from cbe30338-2021; content is available on Github.
A Pyomo implementation of this blending model is shown in the next cell. The model is contained within a Python function so that it can be more easily reused for additional calculations, or eventually for use by the process operator.
Note that the pyomo library has been imported with the prefix pyomo. This is good programming practive to avoid namespace collisions with problem data.
import pyomo.environ as pyomo
vol = 100
abv = 0.040
def brew_blend(vol, abv, data):
C = data.keys()
model = pyomo.ConcreteModel()
model.x = pyomo.Var(C, domain=pyomo.NonNegativeReals)
model.cost = pyomo.Objective(expr = sum(model.x[c]*data[c]['cost'] for c in C))
model.vol = pyomo.Constraint(expr = vol == sum(model.x[c] for c in C))
model.abv = pyomo.Constraint(expr = 0 == sum(model.x[c]*(data[c]['abv'] - abv) for c in C))
solver = pyomo.SolverFactory('glpk')
solver.solve(model)
print('Optimal Blend')
for c in data.keys():
print(' ', c, ':', model.x[c](), 'gallons')
print()
print('Volume = ', model.vol(), 'gallons')
print('Cost = $', model.cost())
brew_blend(vol, abv, data)
WARNING: Could not locate the 'glpsol' executable, which is required for
solver 'glpk'
--------------------------------------------------------------------------- ApplicationError Traceback (most recent call last) <ipython-input-222-f1ef07ecaf09> in <module> 27 print('Cost = $', model.cost()) 28 ---> 29 brew_blend(vol, abv, data) <ipython-input-222-f1ef07ecaf09> in brew_blend(vol, abv, data) 18 19 solver = pyomo.SolverFactory('glpk') ---> 20 solver.solve(model) 21 22 print('Optimal Blend') ~/opt/anaconda3/lib/python3.8/site-packages/pyomo/opt/base/solvers.py in solve(self, *args, **kwds) 516 """ Solve the problem """ 517 --> 518 self.available(exception_flag=True) 519 # 520 # If the inputs are models, then validate that they have been ~/opt/anaconda3/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py in available(self, exception_flag) 116 if exception_flag: 117 msg = "No executable found for solver '%s'" --> 118 raise ApplicationError(msg % self.name) 119 return False 120 return True ApplicationError: No executable found for solver 'glpk'
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/(q + 4.0)/(q + 20.0)
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.
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
price = {}
price['Gasoline'] = 72
price['Kerosine'] = 48
price
{'Gasoline': 72, 'Kerosine': 48}
products.loc['Gasoline','Price']
72
# 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/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: .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/user_guide/indexing.html#ix-indexer-is-deprecated if __name__ == '__main__': /Users/jeff/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: .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/user_guide/indexing.html#ix-indexer-is-deprecated # Remove the CWD from sys.path while we load stuff.
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()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-6-54aa9b28c7ab> in <module> ----> 1 processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS) 2 profit = income - raw_materials_cost - processing_cost 3 model.objective = Objective(expr = profit, sense=maximize) 4 5 # constraints <ipython-input-6-54aa9b28c7ab> in <genexpr>(.0) ----> 1 processing_cost = sum(processing_costs[f] * model.x[f] for f in FEEDS) 2 profit = income - raw_materials_cost - processing_cost 3 model.objective = Objective(expr = profit, sense=maximize) 4 5 # constraints NameError: name 'processing_costs' is not defined
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()
File "<ipython-input-7-9d0bbdfd689e>", line 17 product_yield = ^ SyntaxError: invalid syntax
profit()
573517.2413793121
income()
2078344.8275862068
raw_materials_cost()
1464827.5862068948
processing_cost()
39999.99999999996
One of the important modeling features of Pyomo is the ability to index variables and constraints. The
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)
{'dime': 0, 'nickel': 1, 'penny': 4, 'quarter': 41}
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
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