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