Pyomo contrib and Debugging Tools

Open this notebook in Colab.

This notebook is part of the SECQUOIA Research Group Pyomo Tutorial. It adapts the workshop’s pyomo.contrib orientation and debugging guidance into small exercises that can be run without a commercial solver.

Learning objectives¶

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

locate installed pyomo.contrib packages and recognize when a capability is an extension;
use infeasibility logging to diagnose violated constraints at current variable values;
detect suspicious scaling before handing a model to a nonlinear solver;
check units consistency in a Pyomo model;
add lightweight timing instrumentation while developing model-building code.

import io
import logging
import pkgutil
import time

import pyomo.environ as pyo
import pyomo.contrib
from pyomo.util.check_units import assert_units_consistent
from pyomo.util.infeasible import log_close_to_bounds, log_infeasible_constraints
from pyomo.util.report_scaling import report_scaling

1. Finding extension capabilities¶

The pyomo.contrib namespace contains tools that build on the core modeling objects: alternative solver interfaces, decomposition utilities, diagnostics, GDP-related algorithms, PyNumero interfaces, robust optimization tools, and more. Availability depends on the installed Pyomo version and optional third-party dependencies.

Exercise 1.1. List installed contrib packages¶

Use Python package discovery to see what is installed in the current environment. Pick two names from the output and look them up in the Pyomo documentation before using them in production code.

contrib_names = sorted(name for _, name, _ in pkgutil.iter_modules(pyomo.contrib.__path__))
print(f"installed pyomo.contrib packages: {len(contrib_names)}")
print(", ".join(contrib_names[:18]))

installed pyomo.contrib packages: 34
alternative_solutions, ampl_function_demo, appsi, aslfunctions, benders, community_detection, cp, cspline_external, doe, example, fbbt, fme, gdp_bounds, gdpopt, gjh, iis, incidence_analysis, interior_point

2. Infeasibility checks¶

Solver failures are easier to diagnose when you first evaluate the model at the values currently stored in the variables. The infeasibility utilities report constraints that are violated at those values.

Exercise 2.1. Log infeasible constraints¶

The model below intentionally contains contradictory bounds on x. Capture the Pyomo log output and identify which constraints are violated.

infeasible = pyo.ConcreteModel()
infeasible.x = pyo.Var(initialize=0, bounds=(0, 2))
infeasible.lower_requirement = pyo.Constraint(expr=infeasible.x >= 1)
infeasible.upper_requirement = pyo.Constraint(expr=infeasible.x <= -1)

stream = io.StringIO()
handler = logging.StreamHandler(stream)
logger = logging.getLogger("pyomo.util.infeasible")
old_propagate = logger.propagate
logger.propagate = False
logger.setLevel(logging.INFO)
logger.addHandler(handler)

log_infeasible_constraints(infeasible, log_expression=True, log_variables=True)
log_close_to_bounds(infeasible)

logger.removeHandler(handler)
logger.propagate = old_propagate
print(stream.getvalue())

CONSTR lower_requirement: 1.0 </= 0
  - EXPR: 1.0 </= x
  - VAR x: 0
CONSTR upper_requirement: 0 </= -1.0
  - EXPR: x </= -1.0
  - VAR x: 0
x near LB of 0

3. Scaling checks¶

Poor scaling can make nonlinear solves slow or unreliable. report_scaling uses bounds and expression analysis to flag variables, coefficients, and constraint bodies that may need rescaling.

Exercise 3.1. Report suspicious scaling¶

The following model mixes variables near 1e-7 and 1e6. Use the report to decide which variables or equations should be rescaled.

scaled = pyo.ConcreteModel()
scaled.x = pyo.Var(initialize=1e-7, bounds=(0, 1e-6))
scaled.y = pyo.Var(initialize=1e6, bounds=(0, 1e8))
scaled.balance = pyo.Constraint(expr=1e8 * scaled.x + scaled.y == 1e6)
scaled.obj = pyo.Objective(expr=scaled.y)

stream = io.StringIO()
handler = logging.StreamHandler(stream)
logger = logging.getLogger("pyomo.util.report_scaling")
old_propagate = logger.propagate
logger.propagate = False
logger.setLevel(logging.INFO)
logger.addHandler(handler)

scaling_ok = report_scaling(scaled, too_large=1e4, too_small=1e-5)

logger.removeHandler(handler)
logger.propagate = old_propagate
print(f"scaling check passed: {scaling_ok}")
print(stream.getvalue())

scaling check passed: False


The following variables have large bounds. Please scale them.
          LB          UB    Var
    0.00e+00    1.00e+08    y

The following constraints have potentially large coefficients. Please scale them.
balance
         Coef LB     Coef UB    Var
        1.00e+08    1.00e+08    x

The following constraints have bodies with large bounds. Please scale them.
          LB          UB    Constraint
    0.00e+00    1.00e+08    balance

4. Unit checks¶

Pyomo can attach units to variables, parameters, and expressions. Unit checks catch mistakes before the model is transformed or solved.

Exercise 4.1. Catch inconsistent units¶

The constraint below tries to add a length and a time. Catch the exception and rewrite the expression with consistent units.

units_model = pyo.ConcreteModel()
units_model.length = pyo.Var(initialize=1, units=pyo.units.m)
units_model.time = pyo.Var(initialize=2, units=pyo.units.s)
units_model.bad_balance = pyo.Constraint(
    expr=units_model.length + units_model.time == 3 * pyo.units.m
)

try:
    assert_units_consistent(units_model)
except Exception as err:
    print(type(err).__name__)
    print(str(err).splitlines()[0])

InconsistentUnitsError
Error in units found in expression: length + time: meter not compatible with second.

5. Timing model construction¶

When a model is slow before the solver starts, instrument model construction in small pieces. The exact timing will vary by machine, but relative timing often identifies expensive loops or data transformations.

Exercise 5.1. Time a component-building loop¶

Build a simple indexed constraint and report how many constraints were created per millisecond.

start = time.perf_counter()
timed = pyo.ConcreteModel()
timed.I = pyo.RangeSet(1, 500)
timed.x = pyo.Var(timed.I, bounds=(0, 10))

@timed.Constraint(timed.I)
def upper_envelope(model, i):
    return model.x[i] <= 10 - i / 100

timed.obj = pyo.Objective(expr=sum(timed.x[i] for i in timed.I), sense=pyo.maximize)
elapsed_ms = (time.perf_counter() - start) * 1000
print(f"constraints built: {len(timed.upper_envelope)}")
print(f"elapsed: {elapsed_ms:.2f} ms")

constraints built: 500
elapsed: 1.86 ms

Exercises to continue¶

Run the infeasibility utilities on one of the nonlinear exercise models after setting intentionally poor initial values.
Add units to the reactor-design exercise in Pyomo Nonlinear Problems, then call assert_units_consistent.
Time model construction before and after replacing a Python loop with an indexed Pyomo rule.

References¶

The debugging utilities used here are distributed with Pyomo and complement the modeling guidance in the Pyomo references Bynum et al., 2021.

References¶

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.