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}\\
s.t. \;\;\; y &= \sqrt{x - 1.0}
\end{aligned}\end{split}\]
(a) Below we formulate this Pyomo model and solve using Ipopt. We then get a list of errors from the solver. What did you discover? How might you fix this? (Hint: error output might be ordered strangely, look up in the console output.)
import pyomo.environ as pyo
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))
Error evaluating "var =" definition -1: can't evaluate sqrt(-0.752432).
Ipopt 3.13.2: 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 http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
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[1], line 18
16 solver = pyo.SolverFactory('ipopt')
17 solver.options['halt_on_ampl_error'] = 'yes'
---> 18 solver.solve(model, tee=True)
20 print(pyo.value(model.x))
21 print(pyo.value(model.y))
File ~/CHE498/pyomo_jupyter_book/.venv/lib/python3.10/site-packages/pyomo/opt/base/solvers.py:628, in OptSolver.solve(self, *args, **kwds)
626 elif hasattr(_status, 'log') and _status.log:
627 logger.error("Solver log:\n" + str(_status.log))
--> 628 raise ApplicationError("Solver (%s) did not exit normally" % self.name)
629 solve_completion_time = time.time()
630 if self._report_timing:
ApplicationError: Solver (ipopt) did not exit normally
(b) Add bounds \(x\geq 1\) to fix this problem. Resolve the problem. 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.)
import pyomo.environ as pyo
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.9.1: 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 Common Public License (CPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
NOTE: You are using Ipopt by default with the MUMPS linear solver.
Other linear solvers might be more efficient (see Ipopt documentation).
This is Ipopt version 3.9.1, running with linear solver mumps.
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+001 3.00e+000 8.92e+000 -1.0 0.00e+000 - 0.00e+000 0.00e+000 0
1 1.3964766e+000 9.81e-001 4.57e+000 -1.0 4.20e+000 - 1.00e+000 9.43e-001f 1
2 1.3265288e+000 4.58e-001 9.56e+000 -1.0 3.54e-001 2.0 1.31e-001 1.00e+000f 1
3 3.2528786e-001 1.23e-002 1.25e+000 -1.0 5.70e-001 - 1.00e+000 1.00e+000f 1
4 4.9701812e-003 7.82e-002 1.36e-001 -1.0 3.78e-001 - 1.00e+000 1.00e+000F 1
5 6.0628408e-003 4.53e-002 1.66e+000 -2.5 2.08e-002 - 1.00e+000 1.00e+000h 1
6 6.7063949e-003 4.53e-002 2.24e+000 -2.5 1.05e+000 - 1.53e-002 3.91e-003h 9
7 7.6973110e-003 1.67e-002 5.10e-001 -2.5 6.16e-003 1.5 1.00e+000 1.00e+000h 1
8 3.5855462e-003 1.09e-003 2.51e-001 -2.5 2.81e-002 - 1.00e+000 1.00e+000h 1
9 3.0456986e-003 1.24e-004 1.10e-003 -2.5 4.78e-003 - 1.00e+000 1.00e+000h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 9.2960387e-004 1.58e-002 1.21e+000 -3.8 2.59e-002 - 1.00e+000 1.00e+000h 1
11 2.2722219e-004 1.18e-004 1.42e-001 -3.8 1.75e-002 - 1.00e+000 1.00e+000h 1
12 2.4804811e-004 4.21e-005 4.38e-003 -3.8 8.99e-004 - 1.00e+000 1.00e+000h 1
13 2.5041988e-004 8.02e-007 1.04e-004 -3.8 9.89e-005 - 1.00e+000 1.00e+000h 1
14 1.3928326e-004 4.89e-003 3.54e+000 -5.7 6.12e-003 - 1.00e+000 1.00e+000h 1
15 1.6721345e-004 3.66e-003 4.16e+000 -5.7 1.95e-003 1.0 1.01e-001 1.00e+000h 1
16 1.7756299e-004 1.12e-003 2.28e+000 -5.7 6.51e-004 1.5 1.00e+000 1.00e+000h 1
17 1.5169810e-004 1.68e-005 1.58e-001 -5.7 1.61e-003 1.0 1.00e+000 1.00e+000h 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[1], line 18
16 solver = pyo.SolverFactory('ipopt')
17 solver.options['halt_on_ampl_error'] = 'yes'
---> 18 solver.solve(model, tee=True)
20 print(pyo.value(model.x))
21 print(pyo.value(model.y))
File c:\Users\bmurr\CHE 498\.venv\Lib\site-packages\pyomo\opt\base\solvers.py:628, in OptSolver.solve(self, *args, **kwds)
626 elif hasattr(_status, 'log') and _status.log:
627 logger.error("Solver log:\n" + str(_status.log))
--> 628 raise ApplicationError("Solver (%s) did not exit normally" % self.name)
629 solve_completion_time = time.time()
630 if self._report_timing:
ApplicationError: Solver (ipopt) did not exit normally
(c) Think about other solutions for this problem. (e.g., \(x \geq 1.001\)).
import pyomo.environ as pyo
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.9.1: 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 Common Public License (CPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
NOTE: You are using Ipopt by default with the MUMPS linear solver.
Other linear solvers might be more efficient (see Ipopt documentation).
This is Ipopt version 3.9.1, running with linear solver mumps.
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+001 3.00e+000 8.92e+000 -1.0 0.00e+000 - 0.00e+000 0.00e+000 0
1 1.3985928e+000 9.80e-001 4.51e+000 -1.0 4.20e+000 - 1.00e+000 9.43e-001f 1
2 1.3200328e+000 4.52e-001 1.00e+001 -1.0 3.56e-001 2.0 1.31e-001 1.00e+000f 1
3 3.3527706e-001 1.13e-002 1.27e+000 -1.0 5.60e-001 - 1.00e+000 1.00e+000f 1
4 4.9066403e-003 4.74e-002 2.52e-002 -1.0 3.83e-001 - 1.00e+000 1.00e+000F 1
5 4.4928608e-003 2.36e-002 7.26e-001 -2.5 1.19e-002 - 1.00e+000 1.00e+000h 1
6 5.9528153e-003 7.95e-003 5.53e-001 -2.5 1.04e-002 - 1.00e+000 1.00e+000h 1
7 3.2290130e-003 1.43e-003 4.99e-003 -2.5 2.06e-002 - 1.00e+000 1.00e+000h 1
8 1.5403307e-003 3.85e-003 1.04e-001 -3.8 1.82e-002 - 1.00e+000 1.00e+000h 1
9 1.2315611e-003 2.80e-006 1.34e-003 -3.8 4.28e-003 - 1.00e+000 1.00e+000h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0880355e-003 7.80e-005 2.46e-003 -5.7 2.23e-003 - 1.00e+000 1.00e+000h 1
11 1.0828351e-003 2.74e-010 3.15e-007 -5.7 8.21e-005 - 1.00e+000 1.00e+000h 1
12 1.0809936e-003 1.39e-008 4.39e-007 -8.6 2.96e-005 - 1.00e+000 1.00e+000h 1
13 1.0809927e-003 1.34e-015 6.88e-014 -8.6 1.40e-008 - 1.00e+000 1.00e+000h 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 1.0809926760836025e-003 1.0809926760836025e-003
Dual infeasibility......: 6.8833827526759706e-014 6.8833827526759706e-014
Constraint violation....: 1.3392065234540951e-015 1.3392065234540951e-015
Complementarity.........: 2.5059036424968433e-009 2.5059036424968433e-009
Overall NLP error.......: 2.5059036424968433e-009 2.5059036424968433e-009
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 CPU secs in IPOPT (w/o function evaluations) = 0.005
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
1.001
0.0316226586775465