Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Pyomo Nonlinear Problems

Open this notebook in Colab.

This notebook is part of the SECQUOIA Research Group Pyomo Tutorial. It uses IPOPT for nonlinear solves; see Setup and Solvers for local solver notes. In Colab, the first code cell installs the IDAES extension bundle when IPOPT is missing.

Learning objectives

By the end of this chapter, you should be able to:

  • formulate nonlinear Pyomo objectives and constraints;

  • use bounds and initialization to improve solver behavior;

  • interpret evaluation errors and IPOPT failure modes;

  • compare alternative nonlinear formulations;

  • adapt a reactor-design model from equations to Pyomo code.

import os
import shutil
import subprocess
import sys

import pyomo.environ as pyo
from pyomo.common.errors import ApplicationError

if "google.colab" in sys.modules and shutil.which("ipopt") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "idaes-pse"], check=True)
    subprocess.run(["idaes", "get-extensions"], check=True)
    os.environ["PATH"] = "/root/.idaes/bin:" + os.environ["PATH"]

Example: Rosenbrock function

The Rosenbrock function is defined as below:

f(x,y)=(1−x)2+100(y−x2)2f(x,y) = (1-x)^2 + 100(y-x^2)^2
(1)

Minimize the Rosenbrock function using Pyomo and IPOPT. Initialize xx and yy 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 xx 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 solution for this problem can be found in rosenbrock_init_soln.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.set_value(x_init)
    model.y.set_value(y_init)

    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=5x=5, y=5y=5.

min⁡x,y f(x,y)=(x−1.01)2+y2s.t.  y=x−1.0\begin{aligned} \min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\ \text{s.t.} &~~ y = \sqrt{x - 1.0} \end{aligned}
(2)

  1. Starting with evaluation_error_incomplete.py, 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 locate the domain error. (Hint: error output might be ordered strangely, look up in the console output.) What did you discover? How might you fix this? (A solution can be found in evaluation_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 == 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) )
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'
try:
    solver.solve(model, tee=True)
except ApplicationError as err:
    print(f"Expected solver failure: {err}")
else:
    print( pyo.value(model.x) )
    print( pyo.value(model.y) )
Error evaluating constraint 1: can't evaluate sqrt(-0.752432).

Ipopt 3.14.6: 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.6, 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.

Expected solver failure: Solver (ipopt) did not exit normally

  1. Add bounds x≥1x \geq 1 to fix this problem. Resolve the problem. Comment on the number of iterations and the quality of the solution. (Note: The problem can still occur because x≥1x \geq 1 is enforced within numerical tolerances, so very small violations may remain.) (A solution 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'
try:
    solver.solve(model, tee=True)
except ApplicationError as err:
    print(f"Expected possible solver failure with x >= 1: {err}")
else:
    print( pyo.value(model.x) )
    print( pyo.value(model.y) )
Ipopt 3.14.6: 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.6, 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.

Expected possible solver failure with x >= 1: Solver (ipopt) did not exit normally

  1. Think about other solutions for this problem. (e.g., x≥1.001 x \geq 1.001). (A solution 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.6: 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.6, 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.012

EXIT: Optimal Solution Found.

1.0009999925418367
0.0316226586775465

1.3. Alternative Formulations

Consider the following problem with initial values x=5x=5 and y=5y=5.

min⁡x,y f(x,y)=(x−1.01)2+y2s.t.  x−1y=1\begin{aligned} \min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\ \text{s.t.} &~~ \frac{x-1}{y} = 1 \end{aligned}
(3)

Note that the solution to this problem is x=1.005x=1.005 and y=0.005y=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 formulation and compare the solution quality. What can you learn about problem formulation from these 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):
    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))

Formulation 1. A solution can be found in formulation_1_soln.py.

min⁡x,y f(x,y)=(x−1.01)2+y2s.t.  x−1y=1\begin{aligned} \min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\ \text{s.t.} &~~ \frac{x-1}{y} = 1 \end{aligned}
(4)
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.6: 

******************************************************************************
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.6, 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.7010172e+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.3324854e+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.9193898e+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.1909492e+07 1.78e+00 1.36e+09  -5.7 4.49e+06    -  1.00e+00 1.00e+00f  1
  52r 2.1909492e+07 1.78e+00 9.99e+02   0.2 0.00e+00    -  0.00e+00 4.77e-07R 22
  53r 6.2264652e+07 9.52e-01 6.60e+02   0.2 1.43e+06    -  1.00e+00 2.94e-03f  1
  54  2.9605864e+07 4.56e-02 2.39e+04  -5.7 8.24e+03    -  1.00e+00 5.00e-01f  2
  55  7.5041385e+04 6.23e-01 2.12e+05  -5.7 3.84e+03    -  1.00e+00 1.00e+00f  1
  56  3.8044516e+09 6.20e-01 2.54e+05  -5.7 5.74e+04    -  1.00e+00 1.00e+00h  1
  57  4.8409839e+09 1.71e-01 3.17e+05  -5.7 3.10e+04    -  1.00e+00 1.00e+00h  1
  58  1.4096795e+08 5.45e-01 5.83e+05  -5.7 4.80e+04    -  1.00e+00 1.00e+00f  1
  59  2.0273677e+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.6357816e+10 3.26e-01 1.05e+06  -5.7 5.03e+04    -  1.00e+00 1.00e+00f  1
  61  4.5791473e+09 1.59e-01 2.86e+05  -5.7 5.85e+04    -  1.00e+00 1.00e+00f  1
  62  1.5986896e+09 1.44e-02 6.17e+04  -5.7 4.66e+04    -  1.00e+00 5.00e-01f  2
  63  3.9919790e+08 1.14e-04 2.86e+04  -5.7 2.87e+04    -  1.00e+00 5.00e-01f  2
  64  1.3061186e+00 1.96e+00 6.01e+06  -5.7 1.41e+04    -  1.00e+00 1.00e+00f  1
  65r 1.3061186e+00 1.96e+00 9.99e+02   0.3 0.00e+00    -  0.00e+00 4.77e-07R 22
  66r 2.6019496e+00 9.77e-01 4.87e+02   0.3 1.00e+03    -  1.00e+00 1.94e-03f  1
  67  1.2554949e+00 2.21e-02 4.85e+00  -5.7 1.65e+00    -  1.00e+00 5.00e-01f  2
  68  9.8000074e-04 5.40e-01 7.80e+01  -5.7 7.92e-01    -  1.00e+00 1.00e+00f  1
  69  3.6572624e+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.3470436e+02 1.43e-01 9.63e+01  -5.7 7.69e+00    -  1.00e+00 1.00e+00h  1
  71  1.4059502e+02 2.44e-02 2.93e+01  -5.7 1.44e+01    -  1.00e+00 5.00e-01f  2
  72  3.5786305e+01 5.20e-04 8.68e+00  -5.7 8.51e+00    -  1.00e+00 5.00e-01f  2
  73  7.4391326e-05 3.59e-01 1.54e+02  -5.7 4.23e+00    -  1.00e+00 1.00e+00f  1
  74  5.3687955e+02 3.59e-01 1.78e+02  -5.7 1.95e+01    -  1.00e+00 1.00e+00h  1
  75  5.7064699e+02 7.02e-02 1.09e+02  -5.7 4.95e+00    -  1.00e+00 1.00e+00h  1
  76  1.6408510e+02 6.79e-03 2.64e+01  -5.7 1.67e+01    -  1.00e+00 5.00e-01f  2
  77  4.1148599e+01 3.71e-05 9.08e+00  -5.7 9.10e+00    -  1.00e+00 5.00e-01f  2
  78  5.1696840e-05 3.31e-02 1.12e+01  -5.7 4.54e+00    -  1.00e+00 1.00e+00f  1
  79  1.7828704e-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.4429251e-02 3.91e-03 1.69e+00  -5.7 1.11e-02    -  1.00e+00 1.00e+00h  1
  81  3.7716421e-03 2.69e-04 1.93e-01  -5.7 8.37e-02    -  1.00e+00 5.00e-01h  2
  82  5.0116309e-05 2.32e-03 9.17e-01  -5.7 4.31e-02    -  1.00e+00 1.00e+00h  1
  83  5.0578111e-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.9999999999989716e-05    4.9999999999989716e-05
Dual infeasibility......:   8.3364145447850646e-11    8.3364145447850646e-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.3364145447850646e-11    8.3364145447850646e-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.040

EXIT: Optimal Solution Found.

1.0050000000000423
0.005000000000041267

Formulation 2. A solution can be found in formulation_2_soln.py.

min⁡x,y f(x,y)=(x−1.01)2+y2s.t.  xy+1=1\begin{aligned} \min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\ \text{s.t.} &~~ \frac{x}{y+1} = 1 \end{aligned}
(5)
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.6: 

******************************************************************************
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.6, 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 2.22e-16 4.67e-16  -8.6 4.58e-09    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 9

                                   (scaled)                 (unscaled)
Objective...............:   5.0000000000001622e-05    5.0000000000001622e-05
Dual infeasibility......:   4.6664061503776111e-16    4.6664061503776111e-16
Constraint violation....:   2.2204460492503131e-16    2.2204460492503131e-16
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   4.6664061503776111e-16    4.6664061503776111e-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.008

EXIT: Optimal Solution Found.

1.005
0.005000000000000047

Formulation 3. A solution can be found in formulation_3_soln.py.

min⁡x,y f(x,y)=(x−1.01)2+y2s.t.  y=x−1\begin{aligned} \min_{x,y} &~ f(x,y) = (x - 1.01)^2 + y^2 \\ \text{s.t.} &~~ y = x-1 \\ \end{aligned}
(6)
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.6: 

******************************************************************************
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.6, 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.001

EXIT: Optimal Solution Found.

1.005
0.004999999999999893

Formulation 4. Bounds and initialization can be very helpful when solving nonlinear optimization problems. Starting with formulation_incomplete.py, resolve the original problem with the bound y≥0y \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 in formulation_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.6: 

******************************************************************************
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.6, 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.77e-15 4.02e-13  -8.6 1.58e-06    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 38

                                   (scaled)                 (unscaled)
Objective...............:   5.0000000031647200e-05    5.0000000031647200e-05
Dual infeasibility......:   4.0225288377992996e-13    4.0225288377992996e-13
Constraint violation....:   5.7731597280508140e-15    5.7731597280508140e-15
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5158913589263977e-09    2.5158913589263977e-09
Overall NLP error.......:   2.5158913589263977e-09    2.5158913589263977e-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.043

EXIT: Optimal Solution Found.

1.0050001257911454
0.005000125791145422

1.4. Reactor design problem (Hart et al., 2017; Bequette, 2003)

In this example, we will consider a chemical reactor designed to produce product BB from reactant AA using a reaction scheme known as the Van de Vusse reaction:

A→k1B→k2C2A→k3DA \stackrel{k_1}{\rightarrow} B \stackrel{k_2}{\rightarrow} C \\ 2A \stackrel{k_3}{\rightarrow} D
(7)

Under the appropriate assumptions, FF is the volumetric flowrate through the tank. The concentration of component AA in the feed is cAfc_{Af}, and the concentrations in the reactor are equivalent to the concentrations of each component flowing out of the reactor, given by cAc_A, cBc_B, cCc_C, and cDc_D.

If the reactor is too small, we will not produce sufficient quantity of BB, and if the reactor is too large, much of BB will be further reacted to form the undesired product CC. Therefore, our goal is to solve for the reactor volume that maximizes the outlet concentration for product BB.

The steady-state mole balances for each of the four components are given by,

0=FVcAf−FVcA−k1cA−2k3cA20=−FVcB+k1cA−k2cB0=−FVcC+k2cB0=−FVcD+k3cA2\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}
(8)

The known parameters for the system are,

cAf=10gmolm3k1=56min−1k2=56min−1k3=16000m3mol⋅min.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}}.
(9)

Formulate and solve this optimization problem using Pyomo. Since the volumetric flowrate FF always appears as the numerator over the reactor volume VV, it is common to consider this ratio as a single variable, called the space-velocity SVSV. (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.2588672317142 :  None : False : False : PositiveReals
    cb : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 1072.4372001086317 :  None : False : False : PositiveReals
    cc : Size=1, Index=None
        Key  : Lower : Value            : Upper : Fixed : Stale : Domain
        None :     0 : 1330.09353340888 :  None : False : False : PositiveReals
    cd : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 1861.6051996253866 :  None : False : False : PositiveReals
    sv : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 1.3438117610672788 :  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

References

This notebook adapts material from the public Pyomo Summer Workshop 2018. The Pyomo modeling examples draw on Hart et al., 2017Bynum et al., 2021. The reactor design problem follows Bequette, 2003.

References
  1. Hart, W. E., Laird, C. D., Watson, J.-P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., & Siirola, J. D. (2017). Pyomo–optimization modeling in python (Second, Vol. 67). Springer Science & Business Media.
  2. Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., Watson, J.-P., & Woodruff, D. L. (2021). Pyomo–optimization modeling in python (Third, Vol. 67). Springer Science & Business Media.
  3. Bequette, B. W. (2003). Process Control: Modeling, Design, and Simulation. Prentice Hall.
Pyomo Tutorial
Pyomo Fundamentals
Pyomo Tutorial
Blocks and Model Transformations