Pyomo Nonlinear Problems#
This notebook contains material from the Pyomo Summer Workshop 2018 by the Pyomo Developement Team and the notebook is developed by David Bernal (dbernaln@purdue.edu), Zedong Peng (peng372@purdue.edu), and Albert Lee (lee4382@purdue.edu); the content is available on Github. The original sources of this material are available on Pyomo.
import pyomo.environ as pyo
Example: Rosenbrock function#
The Rosenbrock function is defined as below:
\[
f(x,y) = (1-x)^2 + 100(y-x^2)^2
\]
Minimize the Rosenbrock function using Pyomo and IPOPT. Initialize \(x\) and \(y\) to 1.5 and 1.5, respectively.
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=1.5)
model.y = pyo.Var(initialize=1.5)
def rosenbrock(model):
return (1.0-model.x)**2 \
+ 100.0*(model.y - model.x**2)**2
model.obj = pyo.Objective(rule=rosenbrock, sense=pyo.minimize)
pyo.SolverFactory('ipopt').solve(model)
model.pprint()
2 Var Declarations
x : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : 1.0000000000008233 : None : False : False : Reals
y : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : None : 1.0000000000016314 : None : False : False : Reals
1 Objective Declarations
obj : Size=1, Index=None, Active=True
Key : Active : Sense : Expression
None : True : minimize : (1.0 - x)**2 + 100.0*(y - x**2)**2
3 Declarations: x y obj
Exercise 1: Nonlinear Problems#
1.1. Alternative Initialization#
Effective initialization can be critical for solving nonlinear problems, since they can have several local solutions and numerical difficulties.
Solve the Rosenbrock example 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.
(A soution for this problem can be found in rosenbrock.init.sol.py.)
model = pyo.ConcreteModel()
model.x = pyo.Var()
model.y = pyo.Var()
def rosenbrock(m):
return (1.0-m.x)**2 + 100.0*(m.y - m.x**2)**2
model.obj = pyo.Objective(rule=rosenbrock, sense=pyo.minimize)
solver = pyo.SolverFactory('ipopt')
print('x_init, y_init, x_soln, y_soln')
y_init = 5.0
for x_init in range(2, 6):
model.x = x_init
model.y = 5.0
solver.solve(model)
print("{0:6.2f} {1:6.2f} {2:6.2f} {3:6.2f}".format(x_init, \
y_init, pyo.value(model.x), pyo.value(model.y)))
x_init, y_init, x_soln, y_soln
2.00 5.00 1.00 1.00
3.00 5.00 1.00 1.00
4.00 5.00 1.00 1.00
5.00 5.00 1.00 1.00
1.2. Evaluation errors#
Consider the following problem with initial values \(x=5\), \(y=5\).
\[\begin{split}
\begin{aligned}
\min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\
\text{s.t.} &~~ y = \sqrt{x - 1.0}
\end{aligned}
\end{split}\]
Starting with
evalueation_error_incomplete.py, formulate this Pyomo model and solve using IPOPT. You should get list of errors from the solver. Add the IPOPT solver optionssolver.options[halt_on_ampl_error]=yes`` to find the problem. (Hint: error output might be ordered strangely, look up in the console output.) What did you discover? How might you fix this? (Solution for this can be found inevaluation_error_bounds_soln.py)
# evaluation_error_incomplete.py
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
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')
# TODO: ADD SOLVER OPTIONS HERE
solver.solve(model, tee=True)
print( pyo.value(model.x) )
print( pyo.value(model.y) )
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
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 == pyo.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) )
Ipopt 3.14.12: halt_on_ampl_error=yes
Error evaluating constraint 1: can't evaluate sqrt(-0.752432).
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 2
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 3.00e+00 9.86e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
ERROR: Solver (ipopt) returned non-zero return code (1)
ERROR: See the solver log above for diagnostic information.
---------------------------------------------------------------------------
ApplicationError Traceback (most recent call last)
Cell In[10], line 16
14 solver = pyo.SolverFactory('ipopt')
15 solver.options['halt_on_ampl_error'] = 'yes'
---> 16 solver.solve(model, tee=True)
18 print( pyo.value(model.x) )
19 print( pyo.value(model.y) )
File ~/Github/pyomo/pyomo/opt/base/solvers.py:627, in OptSolver.solve(self, *args, **kwds)
625 elif hasattr(_status, 'log') and _status.log:
626 logger.error("Solver log:\n" + str(_status.log))
--> 627 raise ApplicationError("Solver (%s) did not exit normally" % self.name)
628 solve_completion_time = time.time()
629 if self._report_timing:
ApplicationError: Solver (ipopt) did not exit normally
Add bounds \(x \geq 1\) to fix this problem. Resolve the problem. Comment on the number of iterations and the quality of the solution. (Note: The problem still occurs because \(x \geq 1\) is not enforced exactly, and small numerical values still cause the error.) (A soution for this can be found in
evaluation_error_bounds_soln.py.)
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0, bounds=(1,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 == pyo.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) )
Ipopt 3.14.12: halt_on_ampl_error=yes
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 2
Total number of variables............................: 2
variables with only lower bounds: 1
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 3.00e+00 8.92e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.3964766e+00 9.81e-01 4.57e+00 -1.0 4.20e+00 - 1.00e+00 9.43e-01f 1
2 1.3265288e+00 4.58e-01 9.56e+00 -1.0 3.54e-01 2.0 1.31e-01 1.00e+00f 1
3 3.2528786e-01 1.23e-02 1.25e+00 -1.0 5.70e-01 - 1.00e+00 1.00e+00f 1
4 4.9701812e-03 7.82e-02 1.36e-01 -1.0 3.78e-01 - 1.00e+00 1.00e+00F 1
5 6.0628408e-03 4.53e-02 1.66e+00 -2.5 2.08e-02 - 1.00e+00 1.00e+00h 1
6 6.7063949e-03 4.53e-02 2.24e+00 -2.5 1.05e+00 - 1.53e-02 3.91e-03h 9
7 7.6973110e-03 1.67e-02 5.10e-01 -2.5 6.16e-03 1.5 1.00e+00 1.00e+00h 1
8 3.5855462e-03 1.09e-03 2.51e-01 -2.5 2.81e-02 - 1.00e+00 1.00e+00h 1
9 3.0456986e-03 1.24e-04 1.10e-03 -2.5 4.78e-03 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 9.2960387e-04 1.58e-02 1.21e+00 -3.8 2.59e-02 - 1.00e+00 1.00e+00h 1
11 2.2722219e-04 1.18e-04 1.42e-01 -3.8 1.75e-02 - 1.00e+00 1.00e+00h 1
12 2.4804811e-04 4.21e-05 4.38e-03 -3.8 8.99e-04 - 1.00e+00 1.00e+00h 1
13 2.5041988e-04 8.02e-07 1.04e-04 -3.8 9.89e-05 - 1.00e+00 1.00e+00h 1
14 1.3928326e-04 4.89e-03 3.54e+00 -5.7 6.12e-03 - 1.00e+00 1.00e+00h 1
15 1.6721345e-04 3.66e-03 4.16e+00 -5.7 1.95e-03 1.0 1.01e-01 1.00e+00h 1
16 1.7756299e-04 1.12e-03 2.28e+00 -5.7 6.51e-04 1.5 1.00e+00 1.00e+00h 1
17 1.5169810e-04 1.68e-05 1.58e-01 -5.7 1.61e-03 1.0 1.00e+00 1.00e+00h 1
Error evaluating constraint 1: can't evaluate sqrt(-9.90312e-09).
ERROR: Solver (ipopt) returned non-zero return code (1)
ERROR: See the solver log above for diagnostic information.
---------------------------------------------------------------------------
ApplicationError Traceback (most recent call last)
Cell In[9], line 16
14 solver = pyo.SolverFactory('ipopt')
15 solver.options['halt_on_ampl_error'] = 'yes'
---> 16 solver.solve(model, tee=True)
18 print( pyo.value(model.x) )
19 print( pyo.value(model.y) )
File ~/Github/pyomo/pyomo/opt/base/solvers.py:627, in OptSolver.solve(self, *args, **kwds)
625 elif hasattr(_status, 'log') and _status.log:
626 logger.error("Solver log:\n" + str(_status.log))
--> 627 raise ApplicationError("Solver (%s) did not exit normally" % self.name)
628 solve_completion_time = time.time()
629 if self._report_timing:
ApplicationError: Solver (ipopt) did not exit normally
Think about other solutions for this problem. (e.g., \( x \geq 1.001\)). (A solution for this can be found in
evaluation_error_bounds2_soln.py.)
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 == pyo.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) )
Ipopt 3.14.12: halt_on_ampl_error=yes
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 2
Total number of variables............................: 2
variables with only lower bounds: 1
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 3.00e+00 8.92e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.3985928e+00 9.80e-01 4.51e+00 -1.0 4.20e+00 - 1.00e+00 9.43e-01f 1
2 1.3200328e+00 4.52e-01 1.00e+01 -1.0 3.56e-01 2.0 1.31e-01 1.00e+00f 1
3 3.3527706e-01 1.13e-02 1.27e+00 -1.0 5.60e-01 - 1.00e+00 1.00e+00f 1
4 4.9066403e-03 4.74e-02 2.52e-02 -1.0 3.83e-01 - 1.00e+00 1.00e+00F 1
5 4.4928608e-03 2.36e-02 7.26e-01 -2.5 1.19e-02 - 1.00e+00 1.00e+00h 1
6 5.9528153e-03 7.95e-03 5.53e-01 -2.5 1.04e-02 - 1.00e+00 1.00e+00h 1
7 3.2290130e-03 1.43e-03 4.99e-03 -2.5 2.06e-02 - 1.00e+00 1.00e+00h 1
8 1.5403307e-03 3.85e-03 1.04e-01 -3.8 1.82e-02 - 1.00e+00 1.00e+00h 1
9 1.2315611e-03 2.80e-06 1.34e-03 -3.8 4.28e-03 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0880355e-03 7.80e-05 2.46e-03 -5.7 2.23e-03 - 1.00e+00 1.00e+00h 1
11 1.0828351e-03 2.74e-10 3.15e-07 -5.7 8.21e-05 - 1.00e+00 1.00e+00h 1
12 1.0809936e-03 1.39e-08 4.39e-07 -8.6 2.96e-05 - 1.00e+00 1.00e+00h 1
13 1.0809927e-03 1.34e-15 6.88e-14 -8.6 1.40e-08 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 1.0809926760836025e-03 1.0809926760836025e-03
Dual infeasibility......: 6.8833827526759706e-14 6.8833827526759706e-14
Constraint violation....: 1.3392065234540951e-15 1.3392065234540951e-15
Variable bound violation: 7.4581631981374130e-09 7.4581631981374130e-09
Complementarity.........: 2.5059036424968433e-09 2.5059036424968433e-09
Overall NLP error.......: 2.5059036424968433e-09 2.5059036424968433e-09
Number of objective function evaluations = 15
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 15
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 13
Total seconds in IPOPT = 0.010
EXIT: Optimal Solution Found.
1.0009999925418367
0.0316226586775465
1.3. Alternative Formulations#
Consider the following problem with initial values \(x=5\) and \(y=5\).
\[\begin{split}
\begin{aligned}
\min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\
\text{s.t.} &~~ \frac{x-1}{y} = 1
\end{aligned}
\end{split}\]
Note that 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?
Starting with formulation_incomplete.py finish the Pyomo model for each of the iterations and quality of solutions.
What can you learn about problem formulation from this examples?
# formulation_incomplete.py
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
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):
# TODO: SPECIFY CONSTRAINT HERE
model.con = pyo.Constraint(rule=con_rule)
solver = pyo.SolverFactory('ipopt')
solver.solve(model, tee=True)
print(pyo.value(model.x))
print(pyo.value(model.y))
(A solution can be found in
formulation_1_soln.py.) $\( \begin{aligned} \min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\ \text{s.t.} &~~ \frac{x-1}{y} = 1 \end{aligned} \)$
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
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.x - 1.0) / m.y == 1.0
model.con = pyo.Constraint(rule=con_rule)
solver = pyo.SolverFactory('ipopt')
solver.solve(model, tee=True)
print(pyo.value(model.x))
print(pyo.value(model.y))
Ipopt 3.14.12:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 2.00e-01 9.99e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.9762286e-01 2.27e+00 1.52e+01 -1.0 5.48e+00 - 1.00e+00 1.00e+00f 1
2 1.9787064e-01 2.97e+00 2.40e+01 -1.0 8.02e+00 - 1.00e+00 1.25e-01f 4
3 6.5241867e+01 2.82e+00 2.98e+01 -1.0 6.68e+00 - 1.00e+00 1.00e+00h 1
4 9.5583161e+01 1.54e+00 2.31e+01 -1.0 4.70e+00 - 1.00e+00 1.00e+00h 1
5 1.8959813e+02 4.39e-01 2.77e+01 -1.0 1.14e+01 - 1.00e+00 1.00e+00h 1
6 2.2666991e+01 2.02e+00 1.24e+01 -1.0 1.54e+01 - 1.00e+00 1.00e+00f 1
7 3.7976779e+01 9.23e-01 1.21e+01 -1.0 3.88e+00 - 1.00e+00 1.00e+00h 1
8 2.0620942e+01 6.84e-02 1.64e+01 -1.0 5.65e+00 - 1.00e+00 5.00e-01f 2
9 9.8180623e-02 8.02e-01 1.24e+02 -1.0 3.63e+00 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.3622224e+03 7.97e-01 1.42e+02 -1.0 4.73e+01 - 1.00e+00 1.00e+00h 1
11 3.6759509e+03 1.64e-01 2.37e+02 -1.0 3.63e+01 - 1.00e+00 1.00e+00h 1
12 7.0378669e+01 6.52e-01 5.20e+02 -1.0 4.32e+01 - 1.00e+00 1.00e+00f 1
13 2.5142436e+04 6.17e-01 5.94e+02 -1.0 1.40e+02 - 1.00e+00 1.00e+00h 1
14 2.6505724e+04 3.06e-01 1.02e+03 -1.7 7.26e+01 - 1.00e+00 1.00e+00h 1
15 3.7335004e+03 2.99e-01 5.43e+02 -1.7 9.42e+01 - 1.00e+00 1.00e+00f 1
16 3.5127061e+03 7.10e-02 3.16e+02 -1.7 9.60e+00 - 1.00e+00 1.00e+00f 1
17 3.2134630e+01 4.96e-01 1.50e+03 -1.7 4.08e+01 - 1.00e+00 1.00e+00f 1
18 1.0915648e+05 4.88e-01 1.76e+03 -1.7 2.89e+02 - 1.00e+00 1.00e+00h 1
19 1.1667899e+05 1.53e-01 1.76e+03 -1.7 1.07e+02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 3.7449172e+03 4.78e-01 3.37e+03 -1.7 2.30e+02 - 1.00e+00 1.00e+00f 1
21 4.9463570e+05 4.36e-01 3.65e+03 -1.7 5.58e+02 - 1.00e+00 1.00e+00h 1
22 3.6988313e+05 2.75e-01 5.13e+03 -1.7 2.37e+02 - 1.00e+00 1.00e+00f 1
23 8.9309627e+04 1.73e-01 1.92e+03 -1.7 2.88e+02 - 1.00e+00 1.00e+00f 1
24 4.0590995e+04 3.57e-02 4.55e+02 -1.7 1.71e+02 - 1.00e+00 5.00e-01f 2
25 1.8907970e+01 5.44e+00 1.24e+05 -1.7 1.46e+02 - 1.00e+00 1.00e+00f 1
26 1.4670968e+01 5.44e+00 1.24e+05 -1.7 1.63e+04 - 1.00e+00 4.88e-04f 12
27 2.3540992e+01 5.44e+00 1.24e+05 -1.7 1.63e+04 - 1.00e+00 6.10e-05h 15
28 2.8758563e+01 5.44e+00 1.24e+05 -1.7 1.63e+04 - 1.00e+00 3.05e-05h 16
29 3.1562971e+01 5.44e+00 1.24e+05 -1.7 1.63e+04 - 1.00e+00 1.53e-05h 17
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 2.8032208e+08 5.44e+00 1.32e+05 -1.7 1.63e+04 - 1.00e+00 1.00e+00s 22
31 2.9975467e+08 2.93e+00 7.23e+04 -1.7 4.29e+03 - 1.00e+00 1.00e+00s 22
32 5.0495436e+08 1.67e+00 5.58e+04 -1.7 1.07e+04 - 1.00e+00 1.00e+00s 22
33r 5.0495436e+08 1.67e+00 9.99e+02 0.2 0.00e+00 - 0.00e+00 0.00e+00R 1
34r 1.7350836e+09 9.31e-01 7.02e+02 0.2 6.83e+06 - 1.00e+00 3.35e-03f 1
35r 9.2739578e+08 4.28e-01 8.63e+00 0.2 4.56e+04 - 1.00e+00 3.31e-01f 1
36 2.5595903e+08 9.68e-02 4.45e+04 -1.7 3.13e+04 - 1.00e+00 5.00e-01f 2
37 9.5896907e+05 9.69e-01 2.09e+05 -2.5 1.18e+04 - 1.00e+00 1.00e+00f 1
38 1.0215080e+10 9.60e-01 2.15e+05 -2.5 1.00e+05 - 1.00e+00 1.00e+00h 1
39 1.9190691e+10 5.97e-02 4.31e+05 -2.5 9.67e+04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 3.6729262e+07 8.83e-01 3.34e+06 -2.5 1.00e+05 - 1.00e+00 1.00e+00f 1
41 2.1363532e+12 8.79e-01 3.67e+06 -2.5 1.45e+06 - 1.00e+00 1.00e+00h 1
42 3.7010171e+12 1.13e-01 6.64e+06 -2.5 1.26e+06 - 1.00e+00 1.00e+00h 1
43 2.8512383e+10 7.70e-01 2.44e+07 -2.5 1.39e+06 - 1.00e+00 1.00e+00f 1
44 8.3324853e+13 7.56e-01 2.81e+07 -2.5 8.70e+06 - 1.00e+00 1.00e+00h 1
45 1.1867571e+14 2.10e-01 4.78e+07 -2.5 6.23e+06 - 1.00e+00 1.00e+00h 1
46 4.3825285e+12 5.52e-01 6.94e+07 -2.5 7.54e+06 - 1.00e+00 1.00e+00f 1
47 2.9193897e+14 4.83e-01 7.29e+07 -2.5 1.33e+07 - 1.00e+00 1.00e+00h 1
48 1.9270328e+14 4.45e-01 1.55e+08 -2.5 7.28e+06 - 1.00e+00 1.00e+00f 1
49 1.5322309e+14 3.60e-02 4.23e+07 -2.5 2.51e+06 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 4.0213041e+13 1.47e-03 1.01e+07 -2.5 8.84e+06 - 1.00e+00 5.00e-01f 2
51 2.1909501e+07 1.78e+00 1.36e+09 -5.7 4.49e+06 - 1.00e+00 1.00e+00f 1
52r 2.1909501e+07 1.78e+00 9.99e+02 0.2 0.00e+00 - 0.00e+00 4.77e-07R 22
53r 6.2264681e+07 9.52e-01 6.60e+02 0.2 1.43e+06 - 1.00e+00 2.94e-03f 1
54 2.9605877e+07 4.56e-02 2.39e+04 -5.7 8.24e+03 - 1.00e+00 5.00e-01f 2
55 7.5041437e+04 6.23e-01 2.12e+05 -5.7 3.84e+03 - 1.00e+00 1.00e+00f 1
56 3.8044524e+09 6.20e-01 2.54e+05 -5.7 5.74e+04 - 1.00e+00 1.00e+00h 1
57 4.8409848e+09 1.71e-01 3.17e+05 -5.7 3.10e+04 - 1.00e+00 1.00e+00h 1
58 1.4096798e+08 5.45e-01 5.83e+05 -5.7 4.80e+04 - 1.00e+00 1.00e+00f 1
59 2.0273680e+10 4.99e-01 6.42e+05 -5.7 1.16e+05 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 1.6357818e+10 3.26e-01 1.05e+06 -5.7 5.03e+04 - 1.00e+00 1.00e+00f 1
61 4.5791481e+09 1.59e-01 2.86e+05 -5.7 5.85e+04 - 1.00e+00 1.00e+00f 1
62 1.5986897e+09 1.44e-02 6.17e+04 -5.7 4.66e+04 - 1.00e+00 5.00e-01f 2
63 3.9919794e+08 1.14e-04 2.86e+04 -5.7 2.87e+04 - 1.00e+00 5.00e-01f 2
64 1.3061165e+00 1.96e+00 6.01e+06 -5.7 1.41e+04 - 1.00e+00 1.00e+00f 1
65r 1.3061165e+00 1.96e+00 9.99e+02 0.3 0.00e+00 - 0.00e+00 4.77e-07R 22
66r 2.6019454e+00 9.77e-01 4.87e+02 0.3 1.00e+03 - 1.00e+00 1.94e-03f 1
67 1.2554928e+00 2.21e-02 4.85e+00 -5.7 1.65e+00 - 1.00e+00 5.00e-01f 2
68 9.7999887e-04 5.40e-01 7.80e+01 -5.7 7.92e-01 - 1.00e+00 1.00e+00f 1
69 3.6572559e+02 5.39e-01 9.38e+01 -5.7 1.73e+01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 4.3470356e+02 1.43e-01 9.63e+01 -5.7 7.69e+00 - 1.00e+00 1.00e+00h 1
71 1.4059476e+02 2.44e-02 2.93e+01 -5.7 1.44e+01 - 1.00e+00 5.00e-01f 2
72 3.5786237e+01 5.20e-04 8.68e+00 -5.7 8.51e+00 - 1.00e+00 5.00e-01f 2
73 7.4391281e-05 3.59e-01 1.54e+02 -5.7 4.23e+00 - 1.00e+00 1.00e+00f 1
74 5.3687560e+02 3.59e-01 1.78e+02 -5.7 1.95e+01 - 1.00e+00 1.00e+00h 1
75 5.7064266e+02 7.02e-02 1.09e+02 -5.7 4.95e+00 - 1.00e+00 1.00e+00h 1
76 1.6408381e+02 6.79e-03 2.64e+01 -5.7 1.67e+01 - 1.00e+00 5.00e-01f 2
77 4.1148273e+01 3.71e-05 9.08e+00 -5.7 9.10e+00 - 1.00e+00 5.00e-01f 2
78 5.1696817e-05 3.31e-02 1.12e+01 -5.7 4.54e+00 - 1.00e+00 1.00e+00f 1
79 1.7827388e-02 3.14e-02 1.08e+01 -5.7 9.57e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 1.4428077e-02 3.91e-03 1.69e+00 -5.7 1.11e-02 - 1.00e+00 1.00e+00h 1
81 3.7713449e-03 2.69e-04 1.93e-01 -5.7 8.37e-02 - 1.00e+00 5.00e-01h 2
82 5.0116310e-05 2.32e-03 9.17e-01 -5.7 4.31e-02 - 1.00e+00 1.00e+00h 1
83 5.0578118e-05 2.18e-04 8.63e-02 -5.7 5.30e-04 - 1.00e+00 1.00e+00h 1
84 4.9998894e-05 2.30e-05 9.01e-03 -5.7 5.28e-04 - 1.00e+00 1.00e+00h 1
85 5.0000001e-05 2.08e-08 8.15e-06 -5.7 4.64e-06 - 1.00e+00 1.00e+00h 1
86 5.0000000e-05 2.07e-13 8.34e-11 -8.6 5.17e-08 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 86
(scaled) (unscaled)
Objective...............: 4.9999999999989723e-05 4.9999999999989723e-05
Dual infeasibility......: 8.3368130107674965e-11 8.3368130107674965e-11
Constraint violation....: 2.0738966099997924e-13 2.0738966099997924e-13
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 8.3368130107674965e-11 8.3368130107674965e-11
Number of objective function evaluations = 240
Number of objective gradient evaluations = 86
Number of equality constraint evaluations = 240
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 90
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 86
Total seconds in IPOPT = 0.049
EXIT: Optimal Solution Found.
1.0050000000000423
0.0050000000000412674
(A solution can be found in
formulation_2_soln.py.)
\[\begin{split}
\begin{aligned}
\min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\
\text{s.t.} &~~ \frac{x}{y+1} = 1
\end{aligned}
\end{split}\]
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
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.x / (m.y + 1.0) == 1.0
model.con = pyo.Constraint(rule=con_rule)
solver = pyo.SolverFactory('ipopt')
solver.solve(model, tee=True)
print(pyo.value(model.x))
print(pyo.value(model.y))
Ipopt 3.14.12:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 1.67e-01 9.83e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 4.0023722e-01 1.49e+00 2.82e+01 -1.0 5.39e+00 - 1.00e+00 1.00e+00f 1
2 1.0494629e+01 1.04e+00 2.24e+01 -1.0 2.58e+00 - 1.00e+00 1.00e+00h 1
3 5.4394796e+00 1.78e-01 8.48e+00 -1.0 1.23e+00 - 1.00e+00 1.00e+00f 1
4 1.0538196e-01 1.52e-01 3.90e+00 -1.0 1.75e+00 - 1.00e+00 1.00e+00f 1
5 4.6918851e-02 2.17e-02 1.20e+00 -1.0 1.63e-01 - 1.00e+00 1.00e+00h 1
6 1.3071752e-04 2.88e-03 1.19e-01 -1.0 1.62e-01 - 1.00e+00 1.00e+00h 1
7 4.9833353e-05 1.79e-05 1.02e-03 -2.5 6.25e-03 - 1.00e+00 1.00e+00h 1
8 5.0000013e-05 1.28e-09 6.30e-08 -5.7 9.00e-05 - 1.00e+00 1.00e+00h 1
9 5.0000000e-05 0.00e+00 3.30e-16 -8.6 4.58e-09 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 4.9999999999999169e-05 4.9999999999999169e-05
Dual infeasibility......: 3.2959746043559335e-16 3.2959746043559335e-16
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 3.2959746043559335e-16 3.2959746043559335e-16
Number of objective function evaluations = 10
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 10
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 10
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 9
Total seconds in IPOPT = 0.007
EXIT: Optimal Solution Found.
1.0050000000000001
0.0050000000000000235
(A solution can be found in
formulation_3_soln.py.)
\[\begin{split}
\begin{aligned}
\min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\
\text{s.t.} &~~ y = x-1 \\
\end{aligned}
\end{split}\]
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
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 == m.x - 1.0
model.con = pyo.Constraint(rule=con_rule)
solver = pyo.SolverFactory('ipopt')
solver.solve(model, tee=True)
print(pyo.value(model.x))
print(pyo.value(model.y))
Ipopt 3.14.12:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 2
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 1.00e+00 8.99e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.0000000e-05 0.00e+00 4.44e-16 -1.0 5.00e+00 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 5.0000000000000090e-05 5.0000000000000090e-05
Dual infeasibility......: 4.4408920985006262e-16 4.4408920985006262e-16
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 4.4408920985006262e-16 4.4408920985006262e-16
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 0.003
EXIT: Optimal Solution Found.
1.005
0.004999999999999893
Bounds and initialization can be very helpful when solving nonlinear optimization problems. Starting with
formulation_incomplete.py, resolver the original problem, but add bounds, \(y \geq 0\). Note the number of iterations and quality of solution, and compare with what you found in 1.2.2. (A solution can be found informulation_4_soln.py.)
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=5.0)
model.y = pyo.Var(initialize=5.0, bounds=(0,None))
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.x - 1.0) / m.y == 1.0
model.con = pyo.Constraint(rule=con_rule)
solver = pyo.SolverFactory('ipopt')
solver.solve(model, tee=True)
print(pyo.value(model.x))
print(pyo.value(model.y))
Ipopt 3.14.12:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.12, running with linear solver MUMPS 5.2.1.
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 1
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 4.0920100e+01 2.00e-01 9.38e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 9.7049410e-01 1.89e+01 2.29e+03 -1.0 5.19e+00 - 1.00e+00 9.54e-01f 1
2 2.8303563e+03 1.85e+01 2.25e+03 -1.0 5.21e+01 - 1.88e-02 1.00e+00h 1
3 2.9704435e+03 9.12e+00 1.14e+03 -1.0 2.64e+00 - 1.00e+00 1.00e+00h 1
4 2.5293245e+03 4.09e+00 5.51e+02 -1.0 4.89e+00 - 8.94e-01 1.00e+00f 1
5 1.7547126e+03 1.52e+00 2.52e+02 -1.0 1.04e+01 - 1.00e+00 1.00e+00f 1
6 6.8047428e+02 1.49e-01 8.49e+01 -1.0 1.93e+01 - 1.00e+00 1.00e+00f 1
7 7.4721920e+00 8.01e-01 2.32e+02 -1.0 1.91e+01 - 1.00e+00 1.00e+00f 1
8 8.1969347e+03 7.77e-01 2.62e+02 -1.0 8.57e+01 - 3.37e-02 1.00e+00h 1
9 1.1605265e+04 2.36e-01 4.89e+02 -1.0 6.40e+01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 5.5411085e+02 5.19e-01 5.83e+02 -1.0 7.36e+01 - 1.00e+00 1.00e+00f 1
11 1.7377623e+04 4.22e-01 5.65e+02 -1.0 9.29e+01 - 3.00e-01 1.00e+00h 1
12 1.1573362e+04 1.25e-01 1.63e+02 -1.0 6.64e+01 - 1.00e+00 5.00e-01f 2
13 3.3128549e+03 8.55e-03 9.20e+01 -1.0 8.52e+01 - 1.00e+00 4.70e-01f 2
14 7.1443870e-01 8.34e-01 1.80e+03 -1.0 4.11e+01 - 1.00e+00 9.85e-01f 1
15 3.2462198e-05 4.75e-01 1.05e+03 -1.0 1.73e+02 - 1.00e+00 4.30e-03f 1
16 9.8963008e-05 2.74e-01 8.52e+02 -1.0 1.41e+00 - 1.00e+00 4.23e-03h 1
17 9.2105755e-05 2.40e-01 5.57e+02 -1.0 3.56e-04 - 1.00e+00 1.00e+00f 1
18 7.9772743e-05 1.58e-01 3.60e+02 -1.0 7.13e-04 - 1.00e+00 1.00e+00h 1
19 6.9741224e-05 6.86e-02 1.47e+02 -1.0 7.40e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 4.9888137e-05 3.97e-02 6.88e+01 -1.0 2.39e-03 - 1.00e+00 1.00e+00h 1
21 3.3950490e-04 3.02e-02 4.89e+01 -1.0 1.32e-02 - 1.00e+00 1.00e+00f 1
22 3.7536170e-03 1.95e-02 3.11e+01 -1.0 3.12e-02 - 1.00e+00 1.00e+00h 1
23 1.2981580e-02 8.58e-03 1.36e+01 -1.0 3.74e-02 - 1.00e+00 1.00e+00h 1
24 3.4681398e-02 3.23e-03 5.01e+00 -1.0 5.13e-02 - 1.00e+00 1.00e+00h 1
25 4.7179466e-02 4.50e-04 7.03e-01 -1.0 2.21e-02 - 1.00e+00 1.00e+00h 1
26 1.7082443e-02 2.83e-04 3.55e-01 -1.7 6.13e-02 - 1.00e+00 1.00e+00f 1
27 1.0216180e-02 7.79e-05 1.09e-01 -1.7 2.10e-02 - 1.00e+00 1.00e+00h 1
28 3.1289054e-03 5.64e-05 3.27e-02 -2.5 3.21e-02 - 1.00e+00 1.00e+00h 1
29 1.4987215e-03 2.18e-05 2.21e-02 -2.5 1.23e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 3.8554063e-04 1.69e-05 4.50e-03 -3.8 1.40e-02 - 1.00e+00 1.00e+00h 1
31 1.3832051e-04 9.18e-06 4.39e-03 -3.8 6.31e-03 - 1.00e+00 1.00e+00h 1
32 8.9968997e-05 2.11e-06 9.99e-06 -3.8 2.18e-03 - 1.00e+00 1.00e+00h 1
33 5.4301800e-05 9.79e-07 2.80e-03 -5.7 3.00e-03 - 1.00e+00 1.00e+00h 1
34 5.0216846e-05 2.09e-07 1.11e-04 -5.7 1.14e-03 - 1.00e+00 1.00e+00h 1
35 5.0020268e-05 9.37e-09 2.55e-05 -5.7 2.29e-04 - 1.00e+00 1.00e+00h 1
36 5.0016425e-05 1.85e-11 1.06e-08 -5.7 1.00e-05 - 1.00e+00 1.00e+00h 1
37 5.0000006e-05 3.39e-13 3.16e-06 -8.6 8.89e-05 - 1.00e+00 1.00e+00h 1
38 5.0000000e-05 5.55e-15 4.02e-13 -8.6 1.58e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 38
(scaled) (unscaled)
Objective...............: 5.0000000031647193e-05 5.0000000031647193e-05
Dual infeasibility......: 4.0225288377992996e-13 4.0225288377992996e-13
Constraint violation....: 5.5511151231257827e-15 5.5511151231257827e-15
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5158913589264014e-09 2.5158913589264014e-09
Overall NLP error.......: 2.5158913589264014e-09 2.5158913589264014e-09
Number of objective function evaluations = 44
Number of objective gradient evaluations = 39
Number of equality constraint evaluations = 44
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 39
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 38
Total seconds in IPOPT = 0.020
EXIT: Optimal Solution Found.
1.0050001257911454
0.005000125791145421
1.4. Reactor design problem (Hart et al., 2017; Bequette, 2003)#
In this example, 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:
\[\begin{split}
A \stackrel{k_1}{\rightarrow} B \stackrel{k_2}{\rightarrow} C \\
2A \stackrel{k_3}{\rightarrow} D
\end{split}\]
Under the appropriate assumptions, \(F\) is the volumetric flowrate through the tank. The concentration 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,
\[\begin{split}
\begin{aligned}
0 =& \frac{F}{V} c_{Af} - \frac{F}{V} c_{A} - k_1 c_A - 2 k_3 c_A^2 \\
0 =& - \frac{F}{V} c_{B} + k_1 c_A - k_2 c_B \\
0 =& - \frac{F}{V} c_{C} + k_2 c_B \\
0 =& - \frac{F}{V} c_{D} + k_3 c_A^2
\end{aligned}
\end{split}\]
The known parameters for the system are,
\[
c_{Af} = 10 \frac{\text{gmol}}{\text{m}^3} \quad k_1 = \frac{5}{6} \text{min}^{-1} \quad k_2 = \frac{5}{6} \text{min}^{-1} \quad k_3 = \frac{1}{6000} \frac{\text{m}^3}{\text{mol} \cdot \text{min}}.
\]
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\).
(A solution to this problem can be found in reactor_design_soln.py.)
# create the concrete model
model = pyo.ConcreteModel()
# set the data (native python data)
k1 = 5.0/6.0 # min^-1
k2 = 5.0/3.0 # min^-1
k3 = 1.0/6000.0 # m^3/(gmol min)
caf = 10000.0 # gmol/m^3
# create the variables
model.sv = pyo.Var(initialize = 1.0, within=pyo.PositiveReals)
model.ca = pyo.Var(initialize = 5000.0, within=pyo.PositiveReals)
model.cb = pyo.Var(initialize = 2000.0, within=pyo.PositiveReals)
model.cc = pyo.Var(initialize = 2000.0, within=pyo.PositiveReals)
model.cd = pyo.Var(initialize = 1000.0, within=pyo.PositiveReals)
# create the objective
model.obj = pyo.Objective(expr = model.cb, sense=pyo.maximize)
# create the constraints
model.ca_bal = pyo.Constraint(expr = (0 == model.sv * caf \
- model.sv * model.ca - k1 * model.ca \
- 2.0 * k3 * model.ca ** 2.0))
model.cb_bal = pyo.Constraint(expr=(0 == -model.sv * model.cb \
+ k1 * model.ca - k2 * model.cb))
model.cc_bal = pyo.Constraint(expr=(0 == -model.sv * model.cc \
+ k2 * model.cb))
model.cd_bal = pyo.Constraint(expr=(0 == -model.sv * model.cd \
+ k3 * model.ca ** 2.0))
pyo.SolverFactory('ipopt').solve(model)
model.pprint()
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.6051996253873 : None : False : False : PositiveReals
sv : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : 0 : 1.3438117610672777 : 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