Pyomo Homework 2#
# This code cell installs packages on Colab
import sys
if "google.colab" in sys.modules:
!wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
import helper
helper.easy_install()
helper.install_glpk()
else:
sys.path.insert(0, '../')
import helper
helper.set_plotting_style()
## IMPORT LIBRARIES
import pyomo.environ as pyo
import pandas as pd
Special thanks to the Pyomo team for create these excercises as part of their excellent PyomoFest workshop.
Some Advanced Pyomo Tricks#
Using the decorator notation for rules#
Alternative notation for declaring and defining Pyomo components using decorators exists. Starting with the warehouse location problem code below, change the model to use the decorator notation.
# warehouse_location.py: Warehouse location determination problem
model = pyo.ConcreteModel(name="(WL)")
W = ['Harlingen', 'Memphis', 'Ashland']
C = ['NYC', 'LA', 'Chicago', 'Houston']
d = {('Harlingen', 'NYC'): 1956, \
('Harlingen', 'LA'): 1606, \
('Harlingen', 'Chicago'): 1410, \
('Harlingen', 'Houston'): 330, \
('Memphis', 'NYC'): 1096, \
('Memphis', 'LA'): 1792, \
('Memphis', 'Chicago'): 531, \
('Memphis', 'Houston'): 567, \
('Ashland', 'NYC'): 485, \
('Ashland', 'LA'): 2322, \
('Ashland', 'Chicago'): 324, \
('Ashland', 'Houston'): 1236 }
P = 2
model.x = pyo.Var(W, C, bounds=(0,1))
model.y = pyo.Var(W, within=pyo.Binary)
@model.Objective()
def obj(m):
return sum(d[w,c]*m.x[w,c] for w in W for c in C)
@model.Constraint(C)
def one_per_cust(m, c):
return sum(m.x[w,c] for w in W) == 1
# Add your solution here
def warehouse_active(m, w, c):
return m.x[w,c] <= m.y[w]
# Note: This is only split across cells because of a bug in nbpages (notebook/website software).
# There is no other reason to split your code across cells.
# Add your solution here
def num_warehouses(m):
return sum(m.y[w] for w in W) <= P
pyo.SolverFactory('glpk').solve(model)
model.y.pprint()
model.x.pprint()
Changing parameter values#
A parameter can be specified to be mutable. This tells Pyomo that the value of the parameter may change in the future, and allows the user to change the parameter value and resolve the problem without the need to rebuild the entire model each time. We will use this functionality to find a better solution to the knapsack problem. We would like to find when the wrench becomes valuable enough to be a part of the optimal solution. Create a Pyomo Parameter for the value of the items, make it mutable, and then write a loop that prints the solution for different wrench values.
A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14
model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )
# Add your solution here
def obj_rule(m):
return sum( m.item_benefit[i]*m.x[i] for i in A )
model.obj = pyo.Objective(rule=obj_rule, sense = pyo.maximize )
def weight_rule(m):
return sum( w[i]*m.x[i] for i in A ) <= W_max
model.weight = pyo.Constraint(rule=weight_rule)
# You may instead use 'cbc' as the solver
opt = pyo.SolverFactory('glpk')
for wrench_benefit in range(1,11):
model.item_benefit['wrench'] = wrench_benefit
result_obj = opt.solve(model)
# Add your solution here
Integer cuts#
Often, it can be important to find not only the “best” solution, but a number of solutions that are equally optimal, or close to optimal. For discrete optimization problems, this can be done using something known as an integer cut. Consider again the knapsack problem where the choice of which items to select is a discrete variable \(x_i \forall i \in A\). Let \(x_i^*\) be a particular set of \(x\) values we want to remove from the feasible solution space. We define an integer cut using two sets. The first set \(S_0\) contains the indices for those variables whose current solution is 0, and the second set \(S_1\) consists of indices for those variables whose current solution is 1. Given these two sets, an integer cut constraint that would prevent such a solution from appearing again is defined by,
\(\sum_{i \in S_0}x[i] + \sum_{i \in S_1}(1-x_i) \geq 1\)
Write a loop that solves the problem 5 times, adding an integer cut to remove the previous solution, and printing the value of the objective function and the solution at each iteration of the loop.
from pyomo.environ import *
A = ['hammer', 'wrench', 'screwdriver', 'towel']
b = {'hammer':8, 'wrench':3, 'screwdriver':6, 'towel':11}
w = {'hammer':5, 'wrench':7, 'screwdriver':4, 'towel':3}
W_max = 14
model = pyo.ConcreteModel()
model.x = pyo.Var( A, within=pyo.Binary )
def obj_rule(m):
return sum( b[i]*m.x[i] for i in A )
model.obj = pyo.Objective(rule=obj_rule, sense = maximize )
def weight_con_rule(m):
return sum( w[i]*m.x[i] for i in A ) <= W_max
model.weight_con = Constraint(rule=weight_con_rule)
# You may instead use 'cbc' as the solver
opt = pyo.SolverFactory('glpk')
# create the ConstraintList to hold the integer cuts
model.int_cuts = ConstraintList()
# Add your solution here
Putting it all together: Lot sizing example (Hart et al., 2017)#
We will now write a complete model from scratch using a well-known multi-period optimization problem for optimal lot-sizing adapted from Hagen et al. (2001) shown below.
\(\min \sum_{t \in T}c_ty_t+h_t^+I_t^+ +h_t^-I_t^-\)
s.t. \(I_t=I_{t-1}+X_t-d_t, \forall t \in T\)
\(I_t=I_t^+-I_t^-, \forall t \in T\)
\(X_t \leq Py_t, \forall t \in T\)
\(X_t, I_t^+, I_t^- \geq 0, \forall t \in T\)
\(y_t \in \{0,1\}, \forall t \in T\)
Our goal is to finnd the optimal production \(X_t\) given known demands \(d_t\), fixed cost \(c_t\) associated with active production in a particular time period, an inventory holding cost \(h_t^+\) and a shortage cost \(h_t^-\) (cost of keeping a backlog) of orders. The variable \(y_t\) (binary) determines if we produce in time \(t\) or not, and \(I_t^+\) represents inventory that we are storing across time period \(t\), while \(I_t^-\) represents the magnitude of the backlog. Note that \(X_t \leq Py_t\) is a constraint that only allows production in time period \(t\) if the indicator variable \(y_t\)=1.
Write a Pyomo model for this problem and solve it using glpk (or cbc) using the data provided below.
Parameter |
Description |
Value |
|---|---|---|
\(c\) |
fixed cost of production |
4.6 |
\(I_0^+\) |
initial value of positive inventory |
5.0 |
\(I_0^-\) |
initial value of backlogged orders |
0.0 |
\(h^+\) |
cost (per unit) of holding inventory |
0.7 |
\(h^-\) |
shortage cost (per unit) |
1.2 |
\(P\) |
maximum production amoung (big-M value) |
5 |
\(d\) |
demand |
[5,7,6.2,3.1,1,7] |
Reference: Hart, W. E., Laird, C. D., Watson, J. P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., and Siirola, J. D. Pyomo: Optimization Modeling in Python (Second Edition), Vol (67), Springer Verlag, 2017.
model = pyo.ConcreteModel()
model.T = RangeSet(5) # time periods
i0 = 5.0 # initial inventory
c = 4.6 # setup cost
h_pos = 0.7 # inventory holding cost
h_neg = 1.2 # shortage cost
P = 5.0 # maximum production amount
# demand during period t
d = {1: 5.0, 2:7.0, 3:6.2, 4:3.1, 5:1.7}
# Add your solution here
# solve the problem
# You may instead use 'cbc' as the solver
solver = pyo.SolverFactory('glpk')
solver.solve(model)
# print the results
for t in model.T:
print('Period: {0}, Prod. Amount: {1}'.format(t, value(model.x[t])))
Nonlinear programs: initial and problem formulation are very important!#
Alternative initialization#
Effective initialization can be critical for solving nonlinear problems, since they can have several local solutions and numerical diffculties. Solve the Rosenbrock problem using different initial values for the x (and optionally y) variable(s). Write a loop that varies the initial value from 2.0 to 6.0, solves the problem, and prints the solution for each iteration of the loop.
model = ConcreteModel()
model.x = Var()
model.y = Var()
def rosenbrock(m):
return (1.0-m.x)**2 + 10000.0*(m.y - m.x**2)**2
model.obj = pyo.Objective(rule=rosenbrock, sense=minimize)
solver = pyo.SolverFactory('ipopt')
print('x_init, y_init, x_soln, y_soln')
# Add your solution here
As elaborated here, the Rosenbrock problem is a classic “hard” test case for optimization algorithms. Your results may surprise you (and show the effectiveness of Pyomo and Ipopt!).
Evaluation errors#
Consider the following problem with initial values \(x\)=5, \(y\)=5.
\(\min_{x,y} f(x,y)=(x-1.01)^2+y^2\)
s.t. \(y=\sqrt{x-1.0}\)
Formulate this Pyomo model and solve using IPOPT. You should get a list of errors from the solver. Add the IPOPT solver option solver.options[‘halt_on_ampl_error’]=’yes’ to find the problem. Hint: the error output might be ordered strangely, look this up in the console output. What did you discover? How might you fix this?
# Add your solution here
Question Answers
Fill in here
Add bounds \(x \geq 1\) to fix this problem and resolve. Comment on the number of iterations and the quality of solution. Note, the problem still occurs because \(x \geq 1\) is not enforced exactly, and small numerical values still cause the error.
# Add your solution here
Discussion
Fill in here
Think about other solutions for this problem and attempt to implement one of these solutions. Hint: \(x \geq 1.001\).
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0, bounds=(1.001,None))
model.y = pyo.Var(initialize=5.0)
def obj_rule(m):
return (m.x-1.01)**2 + m.y**2
model.obj = pyo.Objective(rule=obj_rule)
def con_rule(m):
return m.y == sqrt(m.x - 1.0)
model.con = pyo.Constraint(rule=con_rule)
solver = pyo.SolverFactory('ipopt')
solver.options['halt_on_ampl_error'] = 'yes'
solver.solve(model, tee=True)
print( pyo.value(model.x) )
print( pyo.value(model.y) )
Alternative formulations#
Consider the following problem with initial values \(x\)=5, \(y\)=5.
\(\min_{x,y} f(x,y)=(x-1.01)^2+y^2\)
s.t. \(\frac{x-1}{y}=1\)
Note, the solution to this problem is \(x\)=1.005 and \(y\)=0.005. There are several ways that the problem above can be reformulated. Some examples are shown below. Which ones do you expect to be better? Why?
\(\min_{x,y} f(x,y)=(x-1.01)^2+y^2\)
s.t. \(\frac{x-1}{y}=1\)
\(\min_{x,y} f(x,y)=(x-1.01)^2+y^2\)
s.t. \(\frac{x}{y+1}=1\)
\(\min_{x,y} f(x,y)=(x-1.01)^2+y^2\)
s.t. \(y=x-1\)
Implement Pyomo models for each formulation and solve with IPOPT.
Formulation 1#
# Add your solution here
Formulation 2#
# Add your solution here
Formulation 3#
# Add your solution here
Note the nuber of iterations and quality of solutions. What can you learn about the problem formulation from these examples?
Discussion
Fill in here
Bounds and initialization can be very helpful when solving nonlinear optimization problems. Resolve the original problem below, but add bounds, \(y \geq 0\). Note the number of iterations and quality of solution, and compare with what you found for Formulation 1.
# Add your solution here
Discussion
Fill in here
Reactor design problem (Hart et al., 2017; Bequette, 2003)#
Here we will consider a chemical reactor designed to produce product B from reactant A using a reaction scheme known as the Van de Vusse reaction:
\(A \overset{k_1}{\rightarrow} B \overset{k_2}{\rightarrow} C\)
\(2A \overset{k_3}{\rightarrow} D\)
Under appropriate assumptions, \(F\) is the volumetric flowrate through the tank. The concentation of component A in the feed is \(c_{Af}\), and the concentrations in the reactor are equivalent to the concentrations of each component flowing out of the reactor, given by \(c_A\), \(c_B\), \(c_C\), and \(c_D\).
If the reactor is too small, we will not produce sufficient quantity of B, and if the reactor is too large, much of B will be further reacted to form the undesired product C. Therefore, our goal is to solve for the reactor volume that maximizes the outlet concentration for product B.
The steady-state mole balances for each of the four components are given by
\(0=\frac{F}{V}c_{Af}-\frac{F}{V}c_A-k_1C_A-2k_3c^2_A\)
\(0=-\frac{F}{V}c_{B}+k_1C_A-k_2c_B\)
\(0=-\frac{F}{V}c_{C}+k_2c_B\)
\(0=-\frac{F}{V}c_{D}+k_3c^2_A\)
The known parameters for the system are:
\(c_{Af}=10\frac{\mathrm{gmol}}{\mathrm{m}^3}\)
\(k_1=\frac{5}{6}\mathrm{min}^{-1}\)
\(k_2=\frac{5}{3}\mathrm{min}^{-1}\)
\(k_3=\frac{1}{6000}\frac{\mathrm{m}^3}{\mathrm{gmol}~\mathrm{min}^{-1}}\)
Formulate and solve this optimization problem using Pyomo. Since the volumetric flowrate \(F\) always appears as the numerator over the reactor volume \(V\), it is common to consider this ratio as a single variable, called the space-velocity \(SV\).
References: Hart, W. E., Laird, C. D., Watson, J. P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., and Siirola, J. D. Pyomo: Optimization Modeling in Python (Second Edition), Vol (67), Springer Verlag, 2017.
B.W. Bequette. Process control: modeling, design, and simulation. Prentice Hall, 2003.
# Add your solution here