calculate_variable_from_constraint is unreliable for exponential functions

[dallan-keylogic]

Summary

Presently, many models in IDAES rely on the block triangularization provided through solve_strongly_connected_components, which relies on calculate_variable_from_constraint as solver for 1x1 blocks. Recently a grad student tried adding a pH property to an existing property package. However, he ran into an OverflowError when the solver (presumably) tried to solve the associated block. Theoretically, Newton's method converges for any convex function, but it's possible to overflow before that happens.

If we initialized a constraint of the form y == 10 ** -x by hand, it would be trivial to set x.value = -log10(y). Would it be possible to have calculate_variable_from_constraint choose from a set of rules for certain simple types of expression instead of immediately falling back on a Newton method? Alternatively, the line search used could cut back the step length if an overflow is detected.

Steps to reproduce the issue

# example.py
import pyomo.environ as pyo
from pyomo.util.calc_var_value import calculate_variable_from_constraint

m = pyo.ConcreteModel()
m.x = pyo.Var(initialize=5)
m.y = pyo.Var(initialize=1)
m.con = pyo.Constraint(expr=m.y == 10 ** (-m.x))

calculate_variable_from_constraint(m.x, m.con)

Error Message

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
Cell In[1], line 9
      6 m.y = pyo.Var(initialize=1)
      7 m.con = pyo.Constraint(expr=m.y == 10 ** (-m.x))
----> 9 calculate_variable_from_constraint(m.x, m.con)

File ~\miniforge3\envs\idaes-fresh\lib\site-packages\pyomo\util\calc_var_value.py:179, in calculate_variable_from_constraint(variable, constraint, eps, iterlim, linesearch, alpha_min, diff_mode)
    177 if slope:
    178     variable.set_value(-intercept / slope, skip_validation=True)
--> 179     body_val = value(body, exception=False)
    180     if body_val.__class__ not in _invalid_types and abs(body_val - upper) < eps:
    181         # Re-set the variable value to trigger any warnings WRT
    182         # the final variable state
    183         variable.set_value(variable.value)

File ~\miniforge3\envs\idaes-fresh\lib\site-packages\pyomo\common\numeric_types.py:403, in value(obj, exception)
    398         raise
    399 else:
    400     #
    401     # Here, we do not try to catch the exception
    402     #
--> 403     return obj(exception=False)

File ~\miniforge3\envs\idaes-fresh\lib\site-packages\pyomo\core\expr\base.py:118, in ExpressionBase.__call__(self, exception)
    103 def __call__(self, exception=True):
    104     """Evaluate the value of the expression tree.
    105
    106     Parameters
   (...)
    116
    117     """
--> 118     return visitor.evaluate_expression(self, exception)

File ~\miniforge3\envs\idaes-fresh\lib\site-packages\pyomo\core\expr\visitor.py:1314, in evaluate_expression(exp, exception, constant)
   1311         clear_active = True
   1313 try:
-> 1314     ans = visitor.dfs_postorder_stack(exp)
   1315 except (
   1316     TemplateExpressionError,
   1317     ValueError,
   (...)
   1330     #   TypeError: This can be raised in Python3 when evaluating a
   1331     #      operation returns a complex number (e.g., sqrt(-1))
   1332     if exception:

File ~\miniforge3\envs\idaes-fresh\lib\site-packages\pyomo\core\expr\visitor.py:949, in ExpressionValueVisitor.dfs_postorder_stack(self, node)
    945         _result = []
    946 #
    947 # Process the current node
    948 #
--> 949 ans = self.visit(_obj, _result)
    950 if _stack:
    951     #
    952     # "return" the recursion by putting the return value on
    953     # the end of the results stack
    954     #
    955     _stack[-1][-1].append(ans)

File ~\miniforge3\envs\idaes-fresh\lib\site-packages\pyomo\core\expr\visitor.py:1200, in _EvaluationVisitor.visit(self, node, values)
   1198 def visit(self, node, values):
   1199     """Visit nodes that have been expanded"""
-> 1200     return node._apply_operation(values)

File pyomo\\core\\expr\\numeric_expr.pyx:952, in pyomo.core.expr.numeric_expr.PowExpression._apply_operation()

OverflowError: (34, 'Result too large')

Information on your system

Pyomo version: 6.9.0
Python version: 3.10.15
Operating system: Windows 10
How Pyomo was installed (PyPI, conda, source): Probably pip?
Solver (if applicable): n/a

[jsiirola]

calculate_variable_from_constraint is meant to be a fast, lightweight solver. As such it does no interrogation of the underlying expression. My feeling is that adding that kind of processing would be both expensive and fragile (recognizing a particular expression type is effectively as complex as expression simplification - so not cheap - plus I could see the list of "special" structures growing to be quite large).

As an alternative, #3541 just traps the OverflowError and treats it as a general evaluation error. As these were happening in an optional heuristic and in the line search, this allows the core Newton's Method to progress. This is not be a perfect fix (you can always provide an initialization that is so far off you will almost certainly get yourself into an unsolvable state), but it is likely to catch and work around the bulk of these kinds of issues for "well posed" models.

[michaelbynum]

For things like this, we can also start calculate_variable_from_constraint with a call to FBBT. It won't always work, but it should any time the variable being solved for only appears once in the constraint. The nice thing about this approach is that it does not depend on initialization.

import pyomo.environ as pyo
from pyomo.util.calc_var_value import calculate_variable_from_constraint
from pyomo.contrib.fbbt.fbbt import fbbt

m = pyo.ConcreteModel()
m.x = pyo.Var(initialize=5)
m.y = pyo.Var(initialize=1)
m.con = pyo.Constraint(expr=m.y == 10 ** (-m.x))

m.y.fix()
fbbt(m.con)
m.x.pprint()

Output:

x : Size=1, Index=None
    Key  : Lower : Value : Upper : Fixed : Stale : Domain
    None :  -0.0 :     5 :  -0.0 : False : False :  Reals

[michaelbynum]

I should elaborate slightly. Here is what I am suggesting (assuming we are solving for x):

1. Fix all variables except x
2. Remove the bounds on x
3. Run FBBT
4. If x is initialized outside the bounds found from FBBT, move the initial guess to inside the bounds
5. run calculate_variable_from_constraint
6. restore the bounds on x
7. unfix any variables we fixed in step 1

[jsiirola]

@michaelbynum: That seems like a reasonable enhancement, but for performance reasons shouldn't we only do that if an initial call to calculate_variable_from_constraint fails?

[michaelbynum]

Possibly. I'm not sure. FBBT would require walking the expression twice (although that walker is currently slow). If it gives the answer right away, that might not be bad. How many times do we walk the expression for each step of Newton's method?

[jsiirola]

compute_variable_from_constraint has an extra trick before the Newton solve: since most of the time is is actually called on an expression where x appears linearly, it does a spot-check assuming that is the case (and therefore can avoid the entire Newton solve). When this works, the expression tree is only evaluated 2 or 3 times in total (2 if the coefficient on x is 1, and 3 if x appears linearly). According to @andrewlee94, that led to this method being significantly faster than using other tools (I think he was comparing it with something from scipy?).

FBBT would need to walk the tree 3 times, right (the first one would be to identify / fix the other variables)? This still seems like a good idea - especially when either there wasn't an initial guess for x or the guess was "very bad". We could also have a specialized version of FBBT that treats all variables (except for x) as fixed - thereby avoiding the walk implied by step 1.

[michaelbynum]

Oh, nice! That makes perfect sense.

Yes, you are right - it would be three times. Sounds like some careful integration (certainly after the linear trick) could be useful.

[michaelbynum]

For the example in this problem, FBBT + identify_variables is more than 100 times faster than newton's method.

[jsiirola]

Oooo... that's really good. ...And i suppose not surprising - I am pretty sure the cost in the Newton solve is getting the symbolic derivative. Did you (or can you easily) try all three diff_modes?

[michaelbynum]

These are from single runs, so take them with a grain of salt:

• calculate_variable_from_constraint with sympy: 0.3047 seconds
• calculate_variable_from_constraint with reverse_symbolic: 0.0004 seconds
• calculate_variable_from_constraint with reverse_numeric: 0.0005 seconds
• FBBT + identify_variables: 0.0003 seconds

[michaelbynum]

So FBBT is not really faster if you use reverse_symbolic or reverse_numeric, but it could still help with bad initialization.

[Robbybp]

@michaelbynum: That seems like a reasonable enhancement, but for performance reasons shouldn't we only do that if an initial call to calculate_variable_from_constraint fails?

Or have a propagate_bounds=False option that we can use in cases where the extra processing is warranted.