None

This notebook contains material from CBE60499; content is available on Github.

< 2.10 Pyomo Homework 1 | Contents | Tag Index | 2.12 Pyomo Homework 3 >

2.11 Pyomo Homework 2

Due Date: 2/23/2021

# This code cell installs packages on Colab

import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py"
    import helper
    helper.install_idaes()
    helper.install_ipopt()
    helper.install_glpk()

## IMPORT LIBRARIES
from pyomo.environ import *
import pandas as pd

Special thanks to the Pyomo team for create these excercises as part of their excellent PyomoFest workshop.

2.11.1 Some Advanced Pyomo Tricks

2.11.1.1 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 = 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 = Var(W, C, bounds=(0,1))
model.y = Var(W, within=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

# 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.

# YOUR SOLUTION HERE
def num_warehouses(m):
    return sum(m.y[w] for w in W) <= P

SolverFactory('glpk').solve(model)

model.y.pprint()
model.x.pprint()

2.11.1.2 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 = ConcreteModel()
model.x = Var( A, within=Binary )
# YOUR SOLUTION HERE

2.11.1.3 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 = ConcreteModel()
model.x = Var( A, within=Binary )

def obj_rule(m):
    return sum( b[i]*m.x[i] for i in A )
model.obj = 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)

opt = SolverFactory('glpk')

# YOUR SOLUTION HERE

2.11.1.4 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 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 = 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}

# YOUR SOLUTION HERE

# solve the problem
solver = 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]))) 

2.11.2 Nonlinear programs: initial and problem formulation are very important!

2.11.2.1 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 variables. 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 = Objective(rule=rosenbrock, sense=minimize)


solver = SolverFactory('ipopt')

print('x_init, y_init, x_soln, y_soln')

# 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!).

2.11.2.2 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}$

1. 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?

# YOUR SOLUTION HERE

Question Answers

Fill in here

1. 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.

# YOUR SOLUTION HERE

Discussion

Fill in here

1. Think about other solutions for this problem and attempt to implement one of these solutions. Hint: $x \geq 1.001$.

model = ConcreteModel()

model.x = Var(initialize=5.0, bounds=(1.001,None))
model.y = Var(initialize=5.0)

def obj_rule(m):
    return (m.x-1.01)**2 + m.y**2
model.obj = Objective(rule=obj_rule)

def con_rule(m):
    return m.y == sqrt(m.x - 1.0)
model.con = Constraint(rule=con_rule)

solver = SolverFactory('ipopt')
solver.options['halt_on_ampl_error'] = 'yes'
solver.solve(model, tee=True)

print( value(model.x) )
print( value(model.y) )

2.11.2.3 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?

1. $\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $\frac{x-1}{y}=1$

1. $\min_{x,y} f(x,y)=(x-1.01)^2+y^2$

s.t. $\frac{x}{y+1}=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.

2.11.2.3.1 Formulation 1

# YOUR SOLUTION HERE

2.11.2.3.2 Formulation 2

# YOUR SOLUTION HERE

2.11.2.3.3 Formulation 3

# 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.

# YOUR SOLUTION HERE

Discussion

Fill in here

2.11.2.4 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.

# YOUR SOLUTION HERE

5 Var Declarations
    ca : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 3874.2588672317133 :  None : False : False : PositiveReals
    cb : Size=1, Index=None
        Key  : Lower : Value             : Upper : Fixed : Stale : Domain
        None :     0 : 1072.437200108632 :  None : False : False : PositiveReals
    cc : Size=1, Index=None
        Key  : Lower : Value             : Upper : Fixed : Stale : Domain
        None :     0 : 1330.093533408881 :  None : False : False : PositiveReals
    cd : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 1861.6051996253875 :  None : False : False : PositiveReals
    sv : Size=1, Index=None
        Key  : Lower : Value             : Upper : Fixed : Stale : Domain
        None :     0 : 1.343811761067278 :  None : False : False : PositiveReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize :         cb

4 Constraint Declarations
    ca_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                                       : Upper : Active
        None :   0.0 : 10000.0*sv - sv*ca - 0.8333333333333334*ca - 0.0003333333333333333*ca**2.0 :   0.0 :   True
    cb_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                    : Upper : Active
        None :   0.0 : - sv*cb + 0.8333333333333334*ca - 1.6666666666666667*cb :   0.0 :   True
    cc_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                            : Upper : Active
        None :   0.0 : - sv*cc + 1.6666666666666667*cb :   0.0 :   True
    cd_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                                     : Upper : Active
        None :   0.0 : - sv*cd + 0.00016666666666666666*ca**2.0 :   0.0 :   True

10 Declarations: sv ca cb cc cd obj ca_bal cb_bal cc_bal cd_bal

< 2.10 Pyomo Homework 1 | Contents | Tag Index | 2.12 Pyomo Homework 3 >