###############################################################################
# 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#

Created by Prof. Alex Dowling (adowling@nd.edu) at the University of Notre Dame.

This notebook shows how to use the following Degeneracy Hunter 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 Ipopt performance benefits from removing a single redundant constraint

Setup#

We start by importing Pyomo and Degeneracy Hunter.

import pyomo.environ as pyo

from idaes.core.util.model_diagnostics import DegeneracyHunter

Example 1: Well-Behaved Nonlinear Program#

Consider the following “well-behaved” nonlinear optimization problem.

\[\begin{split}\begin{align*} \min_{\mathbf{x}} \quad & \sum_{i=\{0,...,4\}} x_i^2\\ \mathrm{s.t.} \quad & x_0 + x_1 - x_3 \geq 10 \\ & x_0 \times x_3 + x_1 \geq 0 \\ & x_4 \times x_3 + x_0 \times x_3 + x_4 = 0 \end{align*} \end{split}\]

This problem is feasible, well-initialized, and standard constraint qualifications hold. As expected, we have no trouble solving this problem.

Define the model in Pyomo#

We start by defining the optimization problem in Pyomo.

m = pyo.ConcreteModel()

m.I = pyo.Set(initialize=[i for i in range(5)])

m.x = pyo.Var(m.I, bounds=(-10, 10), initialize=1.0)

m.con1 = pyo.Constraint(expr=m.x[0] + m.x[1] - m.x[3] >= 10)
m.con2 = pyo.Constraint(expr=m.x[0] * m.x[3] + m.x[1] >= 0)
m.con3 = pyo.Constraint(expr=m.x[4] * m.x[3] + m.x[0] * m.x[3] - m.x[4] == 0)

m.obj = pyo.Objective(expr=sum(m.x[i] ** 2 for i in m.I))

m.pprint()

Evaluate the initial point#

Initialization is extremely important for nonlinear optimization problems. By setting the Ipopt option max_iter to zero, we can inspect the initial point.

# Specify Ipopt as the solver
opt = pyo.SolverFactory("ipopt")

# Specifying an iteration limit of 0 allows us to inspect the initial point
opt.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)
opt.solve(m, tee=True)

# Create Degeneracy Hunter object
# The Degeneracy Hunter algorithm needs a MILP solver
# Here we specify CBC, an open source solver
dh = DegeneracyHunter(m, solver=pyo.SolverFactory("cbc"))

We expect the exit status Maximum Number of Iterations Exceeded because we told Ipopt to take zero iterations (only evaluate the initial point).

Identify the constraint residuals larger than 0.1#

When developing nonlinear optimization models, one often wants to know: “what constraints are violated at the initial point (or more generally the point the solver terminated) within a given tolerance?” Degeneracy Hunter makes this very easy by provided a simple interface to several IDAES utility functions.

The following line of code will print out all constraints with residuals larger than 0.1:

dh.check_residuals(tol=0.1)

Important: Ipopt does several preprocessing steps when we executed it with zero iterations. When checking the initial point, it is strongly recommended to call Ipopt with zero iterations first. Otherwise, you will not be analyzing the initial point Ipopt starts with.

Identify all variables within 1 of their bounds#

Another common question when developing optimization models is, “Which variables are within their bounds by a given tolerance?” Below is the syntax:

dh.check_variable_bounds(tol=1.0)

Solve the optimization problem#

Now we can solve the optimization problem. We first set the number of iterations to 50 and then resolve with Ipopt.

opt.options["max_iter"] = 50
opt.solve(m, tee=True)

As expected, Ipopt has no trouble solving this optimization problem.

Check if any constraint residuals are large than 1E-14#

Let’s now inspect the new solution to see which (if any) constraints have residuals larger than 10\(^{-14}\).

dh.check_residuals(tol=1e-14)

As expected, all of the constraints are satisfied, even with this fairly tight tolerance.

Identify all variables within 1E-5 of their bounds#

Finally, let’s check if any of the variables are near their bounds at the new solution.

dh.check_variable_bounds(tol=1e-5)

Great, no variables are near their bounds. If a variable was at its bound, it is important the inspect the model and confirm the bound is physically sensible/what you intended.

Check the rank of the constraint Jacobian at the solution#

The main feature of Degeneracy Hunter is to check if an optimization problem is poorly formulated. Let’s see what happens when we check the rank on a carefully formulated optimization problem:

dh.check_rank_equality_constraints()

Example 2: Linear Program with Redundant Equality Constraints#

Now let’s apply Degeneracy Hunter to a poorly formulated optimization problem:

\[\begin{split}\begin{align*} \min_{\mathbf{x}} \quad & \sum_{i=\{1,...,3\}} x_i \\ \mathrm{s.t.}~~& x_1 + x_2 \geq 1 \\ & x_1 + x_2 + x_3 = 1 \\ & x_2 - 2 x_3 \leq 1 \\ & x_1 + x_3 \geq 1 \\ & x_1 + x_2 + x_3 = 1 \\ \end{align*} \end{split}\]

Notice the two equality constraints are redundant. This means the constraint qualifications (e.g., LICQ) do not hold which has three important implications:

  1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)

  2. The calculations performed by the optimization solver may become numerically poorly scaled

  3. Theoretical convergence properties of optimization algorithms may not hold

The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the primary purpose of Degeneracy Hunter. Let’s see it in action.

Define the model in Pyomo#

def example2(with_degenerate_constraint=True):
    """Create the Pyomo model for Example 2

    Arguments:
        with_degenerate_constraint: Boolean, if True, include the redundant linear constraint

    Returns:
        m2: Pyomo model
    """

    m2 = pyo.ConcreteModel()

    m2.I = pyo.Set(initialize=[i for i in range(1, 4)])

    m2.x = pyo.Var(m2.I, bounds=(0, 5), initialize=1.0)

    m2.con1 = pyo.Constraint(expr=m2.x[1] + m2.x[2] >= 1)
    m2.con2 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)
    m2.con3 = pyo.Constraint(expr=m2.x[2] - 2 * m2.x[3] <= 1)
    m2.con4 = pyo.Constraint(expr=m2.x[1] + m2.x[3] >= 1)

    if with_degenerate_constraint:
        m2.con5 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)

    m2.obj = pyo.Objective(expr=sum(m2.x[i] for i in m2.I))

    m2.pprint()

    return m2


# Create the Pyomo model for Example 2 including the redundant constraint
m2 = example2()

Evaluate the initial point#

# Specifying an iteration limit of 0 allows us to inspect the initial point
opt.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)
opt.solve(m2, tee=True)

# Create Degeneracy Hunter object
# The Degeneracy Hunter algorithm needs a MILP solver
# Here we specify CBC, an open source solver
dh2 = DegeneracyHunter(m2, solver=pyo.SolverFactory("cbc"))

Identify constraints with residuals greater than 0.1 at the initial point#

dh2.check_residuals(tol=0.1)

Solve the optimization problem and extract the solution#

Now let’s solve the optimization problem.

opt.options["max_iter"] = 50
opt.solve(m2, tee=True)

for i in m2.I:
    print("x[", i, "]=", m2.x[i]())

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 linesearch 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.

Check the rank of the Jacobian of the equality constraints#

n_deficient = dh2.check_rank_equality_constraints()

A singular value near 0 indicates the Jacobian of the equality constraints is rank deficient. For each near-zero singular value, there is likely one degenerate constraint.

Identify candidate degenerate constraints#

Degeneracy Hunter first identifies candidate degenerate equations.

ds2 = dh2.find_candidate_equations(verbose=True, tee=True)

Find irreducible degenerate sets (IDS)#

Next, Degeneracy Hunter enumerates through the candidate equations. For each candidate equation, Degenerate Hunter solves a MILP to compute the irreducible degenerate set that must contain the candidate equation.

ids = dh2.find_irreducible_degenerate_sets(verbose=True)

Reformulate Example 2#

Now let’s reformulate the model by skipping/removing the redundant equality constraint:

\[\begin{split}\begin{align*} \min_{\mathbf{x}} \quad & \sum_{i=\{1,...,3\}} x_i \\ \mathrm{s.t.}~~& x_1 + x_2 \geq 1 \\ & x_1 + x_2 + x_3 = 1 \\ & x_2 - 2 x_3 \leq 1 \\ & x_1 + x_3 \geq 1 \end{align*} \end{split}\]
m2b = example2(with_degenerate_constraint=False)

Solve the reformulated model#

opt.options["max_iter"] = 50
opt.solve(m2b, tee=True)

for i in m2b.I:
    print("x[", i, "]=", m.x[i]())

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.