Show code cell content
###############################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2023 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
###############################################################################
Degeneracy Hunter Examples#
Author: Prof. Alex Dowling (adowling@nd.edu), University of Notre Dame
Maintainer: Andrew Lee
Updated: 2023-06-01
Alexander W. Dowling and Lorenz T. Biegler (2015). Degeneracy Hunter: An Algorithm for Determining Irreducible Sets of Degenerate Constraints in Mathematical Programs. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering. Ed. by Krist V. Gernaey, Jakob K. Huusom, and Raqul Gani. Computer-Aided Chemical Engineering, Vol. 37, p. 809 – 814. link to PDF
This notebook shows how to use the Degeneracy Hunter tool which is part of the IDAES Diagnostics Toolbox features using two motivating examples:
Inspect constraint violations and bounds of a Pyomo model
Compute the Irreducible Degenerate Set (IDS) for a Pyomo model
Demonstrates the performance benefits in Ipopt from removing a single redundant constraint
What does Degeneracy Mean?#
In simple terms, a degenerate model contains one or more constraints that do not add any additional information that is not already included in other constraints in the model. The simplest example of this is the case where a constraint in the model is duplicated. For example, consider the model below:
Notice the two equality constraints are identical and thus redundant. This means the constraint qualifications that solvers rely on (e.g., Linear Independence Constraint Qualification or LICQ) do not hold which has three important implications:
The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)
The calculations performed by the optimization solver may become numerically poorly scaled
Theoretical convergence properties of optimization algorithms may not hold
Put another way, for a deterministic (square) problem we required that the number of equality constraints be equal to the number of free variables. In this case, there appear to be two equality constraints but in reality there is effectively only one as the two constraints are the same. Thus, the problem is not well defined and there exist a family of solutions that can satisfy the constraints given; any combination of values for x1 and x2 that sum to one will satisfy the constraints as written.
In practice, degeneracies are rarely as straight-forward as a duplicated constraint. A more common occurrence is a set of linearly dependent constraints; that is a system of constraints where one constraint in the set can be formed by a linear combination of other constraints in the set. For example:
Here, the second constraint can be formed by multiplying the first constraint by two. Another common example of linearly dependent constraints is the combination of a set of material balances for all components in a system and an overall material balance. In this case, the overall material balances can be formed by adding all the individual component balances together. Depending on the number of components being modeled, this system of constraints can be quite large, but the end result is the same; the model effectively has one less constraint than it would appear to have.
The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the purpose of Degeneracy Hunter. Let’s see it in action.
Setup#
We start by importing Pyomo and the IDAES Diagnostics Toolbox.
import pyomo.environ as pyo
from idaes.core.util.model_diagnostics import DiagnosticsToolbox
Example: Linear Program with Redundant Equality Constraints#
Let’s apply these tools to a poorly formulated optimization problem:
As you can see, this is similar to our example above with duplicated equality constraints.
The cell below shows how to construct this model in Pyomo.
def build_example(include_redundant_constraint=True):
m = pyo.ConcreteModel()
m.I = pyo.Set(initialize=[i for i in range(1, 4)])
m.x = pyo.Var(m.I, bounds=(0, 5), initialize=1.0)
m.con1 = pyo.Constraint(expr=m.x[1] + m.x[2] >= 1)
m.con2 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)
m.con3 = pyo.Constraint(expr=m.x[2] - 2 * m.x[3] <= 1)
m.con4 = pyo.Constraint(expr=m.x[1] + m.x[3] >= 1)
if include_redundant_constraint:
m.con5 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)
m.obj = pyo.Objective(expr=sum(m.x[i] for i in m.I))
return m
m = build_example()
Check for Structural Issues#
Before we try to solve the model, the first thing we should do is use the Diagnostics Toolbox to check for any structural issues in our model. The cell below creates an instance of the Diagnostics Toolbox and calls the report_structural_issues method.
dt = DiagnosticsToolbox(m)
dt.report_structural_issues()
====================================================================================
Model Statistics
Activated Blocks: 1 (Deactivated: 0)
Free Variables in Activated Constraints: 3 (External: 0)
Free Variables with only lower bounds: 0
Free Variables with only upper bounds: 0
Free Variables with upper and lower bounds: 3
Fixed Variables in Activated Constraints: 0 (External: 0)
Activated Equality Constraints: 2 (Deactivated: 0)
Activated Inequality Constraints: 3 (Deactivated: 0)
Activated Objectives: 1 (Deactivated: 0)
------------------------------------------------------------------------------------
2 WARNINGS
WARNING: 1 Degree of Freedom
WARNING: Structural singularity found
Under-Constrained Set: 3 variables, 2 constraints
Over-Constrained Set: 0 variables, 0 constraints
------------------------------------------------------------------------------------
0 Cautions
No cautions found!
------------------------------------------------------------------------------------
Suggested next steps:
display_underconstrained_set()
====================================================================================
As you can see, the toolbox is reporting one degree of freedom and a structural singularity. In this case, these are both expected as we are trying to solve an optimization problem with one degree of freedom so we can move on to the next step.
Evaluate the initial point#
Before we try to solve our degenerate model, it can often be useful to check the state of the model at the initial step seen by the solver. The cell below creates an instance of Ipopt and sets the iteration limit to 0 so that it terminates at the initial condition.
solver = pyo.SolverFactory("ipopt")
# Specifying an iteration limit of 0 allows us to inspect the initial point
solver.options["max_iter"] = 0
# "Solving" the model with an iteration limit of 0 load the initial point and applies
# any preprocessors (e.g., enforces bounds)
solver.solve(m, tee=True)
WARNING: Could not locate the 'ipopt' executable, which is required for solver
ipopt
---------------------------------------------------------------------------
ApplicationError Traceback (most recent call last)
Cell In[5], line 8
4 solver.options["max_iter"] = 0
6 # "Solving" the model with an iteration limit of 0 load the initial point and applies
7 # any preprocessors (e.g., enforces bounds)
----> 8 solver.solve(m, tee=True)
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/base/solvers.py:534, in OptSolver.solve(self, *args, **kwds)
531 def solve(self, *args, **kwds):
532 """Solve the problem"""
--> 534 self.available(exception_flag=True)
535 #
536 # If the inputs are models, then validate that they have been
537 # constructed! Collect suffix names to try and import from solution.
538 #
539 from pyomo.core.base.block import _BlockData
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:139, in SystemCallSolver.available(self, exception_flag)
137 if exception_flag:
138 msg = "No executable found for solver '%s'"
--> 139 raise ApplicationError(msg % self.name)
140 return False
141 return True
ApplicationError: No executable found for solver 'ipopt'
Identify constraints with large residuals at the initial point#
With the solution at the initial point, we can then use the Diagnostics Toolbox to check the residuals of all constraints at the initial point.
# Check for large residuals
dt.display_constraints_with_large_residuals()
====================================================================================
The following constraint(s) have large residuals (>1.0E-05):
con2: 2.00000E+00
con5: 2.00000E+00
====================================================================================
We can also check for variables which are close to (or outside of) their bounds.
dt.display_variables_near_bounds()
====================================================================================
The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04):
====================================================================================
Solve the optimization problem and extract the solution#
There appear to be no significant issues at the initial point, so let us move on to solving the optimization problem.
solver.options["max_iter"] = 50
solver.solve(m, tee=True)
Ipopt 3.13.2: max_iter=50
******************************************************************************
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...: 6
Number of nonzeros in inequality constraint Jacobian.: 6
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 3
variables with only lower bounds: 0
variables with lower and upper bounds: 3
variables with only upper bounds: 0
Total number of equality constraints.................: 2
Total number of inequality constraints...............: 3
inequality constraints with only lower bounds: 2
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 3.0000000e+00 2.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.3934810e+00 3.93e-01 2.17e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01h 1
2 1.0184419e+00 1.84e-02 1.26e-02 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1
3 1.0011246e+00 1.12e-03 3.82e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1
4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1
5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1
6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 4.99e-04 - 1.00e+00 5.82e-01h 1
7 1.0000788e+00 7.88e-05 4.33e+01 -2.5 1.88e-04 - 1.00e+00 2.94e-01f 2
8 1.0000153e+00 1.53e-05 2.19e+01 -2.5 6.85e-05 - 1.00e+00 8.08e-01h 1
9 1.0000118e+00 1.18e-05 2.78e+02 -2.5 3.47e-05 - 1.00e+00 2.45e-01f 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0000014e+00 1.44e-06 2.83e-08 -2.5 7.79e-06 - 1.00e+00 1.00e+00h 1
11 1.0000000e+00 1.39e-09 1.84e-11 -5.7 2.56e-06 - 1.00e+00 1.00e+00h 1
12 1.0000000e+00 1.39e-09 1.24e+02 -8.6 1.16e-08 - 1.00e+00 9.54e-07h 21
13 9.9999999e-01 9.82e-09 1.14e-13 -8.6 1.33e-08 - 1.00e+00 1.00e+00s 22
14 9.9999999e-01 9.96e-09 6.46e-14 -8.6 2.29e-10 - 1.00e+00 1.00e+00s 22
15 9.9999999e-01 9.96e-09 1.82e+02 -9.0 4.00e-11 - 1.00e+00 9.54e-07h 21
16 9.9999999e-01 1.00e-08 1.52e-13 -9.0 5.89e-11 - 1.00e+00 1.00e+00s 22
Number of Iterations....: 16
(scaled) (unscaled)
Objective...............: 9.9999999000238915e-01 9.9999999000238915e-01
Dual infeasibility......: 1.5188693260016487e-13 1.5188693260016487e-13
Constraint violation....: 9.9976108502985994e-09 9.9976108502985994e-09
Complementarity.........: 9.2217032157601989e-10 9.2217032157601989e-10
Overall NLP error.......: 9.9976108502985994e-09 9.9976108502985994e-09
Number of objective function evaluations = 111
Number of objective gradient evaluations = 17
Number of equality constraint evaluations = 111
Number of inequality constraint evaluations = 111
Number of equality constraint Jacobian evaluations = 17
Number of inequality constraint Jacobian evaluations = 17
Number of Lagrangian Hessian evaluations = 16
Total CPU secs in IPOPT (w/o function evaluations) = 0.000
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.04887509346008301}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled ls. This is the number of line search evaluations. For iterations 0 to 11, ls is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, ls is greater than 20. This means Ipopt is struggling (a little) to converge to the solution.
The cell below shows the values of the variables at the solution.
m.x.display()
x : Size=3, Index=I
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 0 : 1.0000000099975996 : 5 : False : False : Reals
2 : 0 : 0.0 : 5 : False : False : Reals
3 : 0 : 0.0 : 5 : False : False : Reals
Check the rank of the Jacobian of the equality constraints#
One way to check for the presence of degenerate equations is to calculate the rank of the Jacobian matrix of the equality constraints. A singular value near 0 indicates the Jacobian is rank deficient, and for each near-zero singular value there is likely one degenerate constraint.
The Diagnostics Toolbox contains a Singular Value Decomposition (SVD) toolbox which can be used to determine the number of small singular values in the Jacobian for a model as shown below.
svd = dt.prepare_svd_toolbox()
svd.display_rank_of_equality_constraints()
====================================================================================
Number of Singular Values less than 1.0E-06 is 1
====================================================================================
As can be seen above, there is one singular value less than 1e-6, suggesting we have one degenerate constraint.
Finding Irreducible Degenerate Sets (IDS)#
Next, we can use Degeneracy Hunter to identify irreducible degenerate sets; that is, the smallest set of constraints that contain a redundant constraint.
Degeneracy Hunter first identifies candidate degenerate equations by solving a mixed integer linear program (MILP). It then iterates through the candidate equations and solves a second MILP for each to compute the irreducible degenerate set that must contain the candidate equation.
The cell below shows how to create an instance of Degeneracy Hunter and generate a report of the irreducible degenerate sets.
dh = dt.prepare_degeneracy_hunter()
dh.report_irreducible_degenerate_sets()
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Candidate Equations
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model.
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving Candidates MILP model.
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Irreducible Degenerate Sets
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model to compute irreducible degenerate set.
Solving MILP 1 of 2.
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con2.
Solving MILP 2 of 2.
2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con5.
====================================================================================
Irreducible Degenerate Sets
Irreducible Degenerate Set 0
nu Constraint Name
1.0 con2
-1.0 con5
Irreducible Degenerate Set 1
nu Constraint Name
-1.0 con2
1.0 con5
====================================================================================
As can be seen above, Degeneracy Hunter identified two constraints as potentially redundant (the duplicate constraints con2 and con5). Degeneracy Hunter then reports an Irreducible Degenerate Set for each of these constraints, which in the case as identical. More complex models may have partially overlapping IDS’s.
These results show us that con2 and con5 are mutually redundant; i.e., we can eliminate one of these with no effect on the model solution.
Reformulate Model to Remove Redundant Constraint#
Now let’s reformulate the model by removing the redundant equality constraint con5 and solve the reformulated model.
m2 = build_example(include_redundant_constraint=False)
solver.solve(m2, tee=True)
Ipopt 3.13.2: max_iter=50
******************************************************************************
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...: 3
Number of nonzeros in inequality constraint Jacobian.: 6
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 3
variables with only lower bounds: 0
variables with lower and upper bounds: 3
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 3
inequality constraints with only lower bounds: 2
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 3.0000000e+00 2.00e+00 6.30e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.3934810e+00 3.93e-01 1.34e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01f 1
2 1.0184419e+00 1.84e-02 9.15e-03 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1
3 1.0011246e+00 1.12e-03 3.88e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1
4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1
5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1
6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 5.00e-04 - 1.00e+00 5.81e-01h 1
7 1.0000788e+00 7.88e-05 4.34e+01 -2.5 1.88e-04 - 1.00e+00 2.93e-01f 2
8 1.0000154e+00 1.54e-05 2.26e+01 -2.5 6.89e-05 - 1.00e+00 8.04e-01h 1
9 1.0000118e+00 1.18e-05 2.98e+02 -2.5 3.64e-05 - 1.00e+00 2.33e-01f 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0000020e+00 2.04e-06 1.33e+02 -2.5 9.72e-06 - 1.00e+00 8.28e-01h 1
11 1.0000016e+00 1.61e-06 2.26e+03 -2.5 4.79e-06 - 1.00e+00 2.12e-01f 2
12 1.0000002e+00 2.46e-07 8.72e+02 -2.5 1.20e-06 - 1.00e+00 8.47e-01h 1
13 1.0000002e+00 1.95e-07 1.71e+04 -2.5 7.02e-07 - 1.00e+00 2.09e-01f 2
14 1.0000000e+00 1.54e-08 3.23e+03 -2.5 1.50e-07 - 1.00e+00 9.21e-01h 1
15 1.0000000e+00 1.15e-08 9.99e+04 -2.5 6.89e-08 - 1.00e+00 2.54e-01f 2
16 1.0000000e+00 2.22e-16 2.83e-08 -2.5 8.21e-09 - 1.00e+00 1.00e+00h 1
17 1.0000000e+00 2.22e-16 4.14e-11 -8.6 8.25e-16 - 1.00e+00 1.00e+00 0
Number of Iterations....: 17
(scaled) (unscaled)
Objective...............: 9.9999999999999978e-01 9.9999999999999978e-01
Dual infeasibility......: 4.1425156707686206e-11 4.1425156707686206e-11
Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16
Complementarity.........: 2.6658012082875325e-09 2.6658012082875325e-09
Overall NLP error.......: 2.6658012082875325e-09 2.6658012082875325e-09
Number of objective function evaluations = 23
Number of objective gradient evaluations = 18
Number of equality constraint evaluations = 23
Number of inequality constraint evaluations = 23
Number of equality constraint Jacobian evaluations = 18
Number of inequality constraint Jacobian evaluations = 18
Number of Lagrangian Hessian evaluations = 17
Total CPU secs in IPOPT (w/o function evaluations) = 0.002
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 4, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.006506204605102539}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
First, let us check to see if the optimal solution has changed.
m2.x.display()
x : Size=3, Index=I
Key : Lower : Value : Upper : Fixed : Stale : Domain
1 : 0 : 1.0000000005924994 : 5 : False : False : Reals
2 : 0 : 0.0 : 5 : False : False : Reals
3 : 0 : 4.55844318120227e-11 : 5 : False : False : Reals
We get the same answer as before, but careful inspection of the Ipopt output reveals a subtle improvement. Notice ls is only 1 or 2 for all of the iterations, in contrast to more than 20 for the original model. This means Ipopt is taking (nearly) full steps for all iterations.
Let’s also compare the number of function evaluations.
Original model (using Ipopt 3.13.2 with ma27):
Number of objective function evaluations = 111
Number of objective gradient evaluations = 17
Number of equality constraint evaluations = 111
Number of inequality constraint evaluations = 111
Number of equality constraint Jacobian evaluations = 17
Number of inequality constraint Jacobian evaluations = 17
Number of Lagrangian Hessian evaluations = 16
Reformulated model (using Ipopt 3.13.2 with ma27):
Number of objective function evaluations = 23
Number of objective gradient evaluations = 18
Number of equality constraint evaluations = 23
Number of inequality constraint evaluations = 23
Number of equality constraint Jacobian evaluations = 18
Number of inequality constraint Jacobian evaluations = 18
Number of Lagrangian Hessian evaluations = 17
Removing a single redundant constraint reduced the number of objective and constraint evaluations by a factor of 5!
Often degenerate equations have a much worse impact on large-scale problems; for example, degenerate equations can cause Ipopt to require many more iterations or terminate at an infeasible point.