1.2 Evaluation errors:

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