Open mfripp opened 3 years ago
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.
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.
@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.
@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.]
@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.
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 aConstraint
that uses this expression for a late time step, it causes aRecursionError
during presolve. This happens whenconstraint_generator
callsgenerate_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 aRecursionError
. (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
.Error Message
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: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.