Constraints on chained expressions cause RecursionError in generate_standard_repn

[mfripp]

Summary

I have a model that defines the value of a Pyomo Expression for each time step as the value of the same expression in the previous time step plus the value of a decision variable. This is a common structure for models of energy systems with storage, e.g., the state of charge of a battery or the temperature of a water heater.

The Expression seems to be defined OK, but if there are many time steps and I define a Constraint that uses this expression for a late time step, it causes a RecursionError during presolve. This happens when constraint_generator calls generate_standard_repn(constraint_data.body). That sets off a recursive sequence that goes through the same four functions for each time step back to the start. If there are more than 250 time steps, this results in a 1000-layer call stack, which creates a RecursionError. (You need 750 steps to reproduce this in Jupyter, which has a 3000 call recursion limit.)

Steps to reproduce the issue

Running the code below will generate the error RecursionError: maximum recursion depth exceeded.

import pyomo.environ as pyo
from pyomo.repn import generate_standard_repn

m = pyo.ConcreteModel()
m.STEPS = pyo.Set(initialize=range(1000), ordered=True)

m.Rate = pyo.Var(within=pyo.NonNegativeReals)

m.Level = pyo.Expression(
    m.STEPS,
    rule=lambda m, t: (0.5 if t == 0 else m.Level[m.STEPS.prev(t)]) + m.Rate
)

m.Min_Level = pyo.Constraint(rule=lambda m: m.Level[m.STEPS.last()] >= 0.5)

m.Cost = pyo.Objective(rule=lambda m: m.Rate)

opt = pyo.SolverFactory("glpk")
opt.solve(m)

# direct path to the error (called during opt.solve(m)):
# repn = generate_standard_repn(m.Min_Level.body)

Error Message

$ python chain_test.py
Traceback (most recent call last):
  File "chain_test.py", line 21, in <module>
    opt.solve(m)
  File ".../pyomo/opt/base/solvers.py", line 575, in solve
    self._presolve(*args, **kwds)
  File ".../pyomo/opt/solver/shellcmd.py", line 198, in _presolve
    OptSolver._presolve(self, *args, **kwds)
  File ".../pyomo/opt/base/solvers.py", line 672, in _presolve
    self._convert_problem(args,
  File ".../pyomo/opt/base/solvers.py", line 742, in _convert_problem
    return convert_problem(args,
  File ".../pyomo/opt/base/convert.py", line 105, in convert_problem
    problem_files, symbol_map = converter.apply(*tmp, **tmpkw)
  File ".../pyomo/solvers/plugins/converter/model.py", line 84, in apply
    instance.write(
  File ".../pyomo/core/base/block.py", line 1811, in write
    (filename, smap) = problem_writer(self,
  File ".../pyomo/repn/plugins/cpxlp.py", line 161, in __call__
    symbol_map = self._print_model_LP(
  File ".../pyomo/repn/plugins/cpxlp.py", line 610, in _print_model_LP
    for constraint_data, repn in yield_all_constraints():
  File ".../pyomo/repn/plugins/cpxlp.py", line 593, in constraint_generator
    repn = generate_standard_repn(constraint_data.body)
  File ".../pyomo/repn/standard_repn.py", line 372, in generate_standard_repn
    return _generate_standard_repn(expr,
  File ".../pyomo/repn/standard_repn.py", line 983, in _generate_standard_repn
    ans = _collect_standard_repn(expr, 1, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 950, in _collect_standard_repn
    return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 862, in _collect_identity
    return _collect_standard_repn(exp.expr, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 950, in _collect_standard_repn
    return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 477, in _collect_sum
    res_ = _collect_standard_repn(e_, multiplier, idMap,
  File ".../pyomo/repn/standard_repn.py", line 950, in _collect_standard_repn
    return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 862, in _collect_identity
    return _collect_standard_repn(exp.expr, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 950, in _collect_standard_repn
    return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 477, in _collect_sum
    res_ = _collect_standard_repn(e_, multiplier, idMap,

... many repetitions of the previous four calls ...

  File ".../pyomo/repn/standard_repn.py", line 950, in _collect_standard_repn
    return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 477, in _collect_sum
    res_ = _collect_standard_repn(e_, multiplier, idMap,
  File ".../pyomo/repn/standard_repn.py", line 950, in _collect_standard_repn
    return fn(exp, multiplier, idMap, compute_values, verbose, quadratic)
  File ".../pyomo/repn/standard_repn.py", line 855, in _collect_identity
    if exp._args_[0].__class__ in native_numeric_types:
RecursionError: maximum recursion depth exceeded

Information on your system

Pyomo version: 5.7.3
Python version: CPython 3.8.0
Operating system: macOS 10.14.6, Darwin 18.7.0
How Pyomo was installed (PyPI, conda, source): conda or source
Solver (if applicable): glpk

Additional information

As noted in a comment at the end of the code, you can also reproduce the error by calling generate_standard_repn directly on the body of the constraint.

It is possible to workaround this by calculating the Expression directly as a sum of the initial level plus all adjustments since then (without reference to prior levels), like this:

m.Level = pyo.Expression(
    m.STEPS,
    rule=lambda m, t: 0.5 + sum(m.Rate for t_ in m.STEPS if t_ < t)
)

But this is less clear and more error-prone than chaining levels from one step to the next. I had to introduce a dummy time-step variable and make the (possibly incorrect) assumption that members of the STEPS set are sorted.

[mfripp]

After working on this a little more, I've realized that using an Expression for this type of component is a bad idea if there are a lot of time steps, because it ends up creating very long expressions for the later time steps. Model construction time ends up increasing quadratically with the number of time steps, so it gets too slow pretty soon.

Instead, it's better to define m.Level as a Var and then use a Constraint to make sure it has the right value at each time step. This is faster to construct and doesn't lead to the RecursionError, because Pyomo doesn't need to recurse back through all the time steps when it is constructing the model. Here's what that code looks like:

m.Level = pyo.Var(m.STEPS, within=pyo.NonNegativeReals)

m.Set_Level = pyo.Constraint(
    m.STEPS,
    rule=lambda m, t: 
        m.Level[t] == (0.5 if t == m.STEPS.first() else m.Level[m.STEPS.prev(t)]) + m.Rate
)

So for my models, fixing the RecursionError isn't urgent. It may not be important for other people either, since any model with 250+ chained steps in an Expression is going to be pretty big and slow to construct, i.e., the Expression object will have elements with 1, 2, 3, ..., 250+ terms in the sum. That size of object will get too slow to handle pretty soon, even if it is constructed without recursion.

[jsiirola]

Thanks for the detailed analysis! The root problem is that generate_standard_repn() is one of the last remaining recursive expression routines in Pyomo (the rest have been converted to leverage a non-recursive expression walker). For "normal" expressions, long summations are automatically converted to a single "n-ary" SumExpression node. However, by using Expression, you are forcing the expression system to keep each of the sub-trees, resulting in a completely unbalanced expression graph (of alternating unary _ExpressionData nodes and binary SumExpression nodes). The deep tree is also a lot less efficient to traverse, resulting in significantly longer run times.

We have a couple ideas / draft implementations for a new generate_standard_repn() implementation (e.g., see jsiirola/pyomo:standard_repn_rework) that will resolve the recursion error. However, deep unbalanced trees will probably always be measurably slower to work with than shallow / balanced trees.

That said, I am somewhat surprised that the model construction time scales quadratically - that is probably a bug.

EDIT: the quadratic runtime makes perfect sense for your modified model, as at each time step, t, the Expression has to construct a summation that is O(t) terms long.

An alternative approach is to declare a variable indexed by the number of steps; e.g.:

m.Level = Var(m.STEPS)

@m.Constraint(m.STEPS)
def Set_Level(m, t):
    if t == m.STEPS.first():
        return m.Level[t] == 0.5
    else:
        return m.Level[t] == m.Level[m.STEPS.prev(t)] + m.Rate

Assuming that your model is a (MI)LP, you can rely on the solver's presolve step to sort things out and perform variable aggregation as appropriate.

[mfripp]

@jsiirola Thanks for the update on this. I was going to mention that the quadratic issue was in my model itself, not in the Expression code per se, but you already got there in the edit. I realized after I posted the issue that I was trying to use Expression to construct a type of constraint that I shouldn't, where the limit on Level at each step is defined as a sum over all previous steps. That framing has a quadratic dependency on the number of steps.

I went back and looked at our other code (the StateOfCharge calculation in the battery storage component of Switch) and realized that we had already adopted the approach you recommended. That doesn't have the quadratic scaling in Pyomo, and the solver does a fine job of managing those step-by-step linkages without needing to sum from start to end.

By the way, I like the decorator notation for defining Constraint's and Expression's. I'm thinking of switching over to that, but would need to be systematic for a big codebase. Do you find the decorators are easier to explain to new users, or about the same as the old way? When teaching Pyomo, I've settled on just defining a rule function every time, rather than using lambdas, because that seems to be easier for new users to understand. But I think explaining decorators might be tough.

[jsiirola]

@mfripp Most of our tutorial / training material all still uses the explicit rule notation. It is verbose, but has the benefit of being the "least magical." Lambdas and decorators are just a little too magical for new users -- and really don't help them to understand what is going on. The big advantages of the decorator notation is that the declaration is more "linear" (you don't have to declare the rule, then use it to declare a constraint), and more importantly the "same" name doesn't appear in three places (the rule name, the component name, and passing the rule into the component constructor) - which always led to copy-and-paste errors when I would take a constraint, copy it, and modify it -- but forget to update one of the three names. [The syntax was actually developed on a "dare" during the Pyomo workshop at the 2017 INFORMS Computing Society meeting in Austin, where a user was complaining about exactly that problem.]

[mfripp]

@jsiirola Sounds like we're grappling with the same issues (no surprise). I switched for a while to just calling the rule rule, to avoid the repetition. But then I have to explain to new users that it's OK to define a function called rule, pass it to a constructor, then define a different function called rule and use that for something else. And it can make error messages hard to interpret -- I have sometimes seen errors in the rule functions reported using only the name of the function, not the model component or line number.

In Switch, we only use .csv files for data input, so I'm thinking of switching over to using a concrete model with a helper function that can load data directly into parameters and sets as they are defined (possibly subclassing Param and Set to make this an argument to their __init__'s). I think it might be easier for people to understand and debug the code if it runs in one pass, instead of running all the callbacks during create_instance(). This could use either the decorator or explicit rule for the constraints, but they would be run on the spot instead of deferred till later. That would be easier to explain and debug.