Bug when using solve_strongly_connected_components

[andrewlee94]

Summary

I have been attempting to apply the solve_strongly_connected_components method to the IDAES model libraries and have encountered a couple of bugs the underlying Pyomo code. The first issue that comes up appears to be a missing import in pyomo/core/base/external.py, however fixing this unearths a deeper issue.

Steps to reproduce the issue

The best replication I have is to run the IDAES test suite using the following branch: https://github.com/andrewlee94/idaes-pse/tree/bt_solver

One test that reproduces this error is:

pytest idaes/models/properties/modular_properties/examples/tests/test_CO2_bmimPF6_PR.py::TestFlashIntegration::test_initialize

Error Message

The first bug lies in external.py, where it appears that the value method needs to be imported from pyomo.core.expr.numvalue

../../pyomo/pyomo/pyomo/core/base/external.py:214: in evaluate
    args_ = [arg if arg.__class__ in native_types else value(arg)

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
>   args_ = [arg if arg.__class__ in native_types else value(arg) for arg in args]
E   NameError: name 'value' is not defined

../../pyomo/pyomo/pyomo/core/base/external.py:214: NameError

However, fixing this results in a second error from deeper in the code:

idaes/core/base/unit_model.py:527: in initialize
    return solve_strongly_connected_components(
../../pyomo/pyomo/pyomo/contrib/incidence_analysis/util.py:147: in solve_strongly_connected_components
    results = calculate_variable_from_constraint(
../../pyomo/pyomo/pyomo/util/calc_var_value.py:158: in calculate_variable_from_constraint
    expr_deriv = differentiate(expr, wrt=variable,
../../pyomo/pyomo/pyomo/core/expr/calculus/derivatives.py:116: in differentiate
    res = sympy_diff(expr=expr, wrt=wrt, wrt_list=wrt_list)
../../pyomo/pyomo/pyomo/core/expr/calculus/diff_with_sympy.py:47: in differentiate
    objectMap, sympy_expr = sympyify_expression(expr)
../../pyomo/pyomo/pyomo/core/expr/sympy_tools.py:215: in sympyify_expression
    return object_map, visitor.walk_expression(expr)
../../pyomo/pyomo/pyomo/core/expr/visitor.py:256: in walk_expression
    result = self._process_node(root, RECURSION_LIMIT)
../../pyomo/pyomo/pyomo/core/expr/visitor.py:455: in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
../../pyomo/pyomo/pyomo/core/expr/visitor.py:455: in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
../../pyomo/pyomo/pyomo/core/expr/visitor.py:455: in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
../../pyomo/pyomo/pyomo/core/expr/visitor.py:455: in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
../../pyomo/pyomo/pyomo/core/expr/visitor.py:455: in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
../../pyomo/pyomo/pyomo/core/expr/visitor.py:455: in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
../../pyomo/pyomo/pyomo/core/expr/visitor.py:466: in _process_node_bx
    return self.exitNode(node, data)
../../pyomo/pyomo/pyomo/core/expr/sympy_tools.py:152: in exitNode
    return node._apply_operation(values)
../../pyomo/pyomo/pyomo/core/expr/numeric_expr.py:650: in _apply_operation
    return self._fcn.evaluate( result )
../../pyomo/pyomo/pyomo/core/base/external.py:216: in evaluate
    return self._evaluate(args_, None, 0)[0]
../../pyomo/pyomo/pyomo/core/base/external.py:360: in _evaluate
    arglist = _ARGLIST(args, fgh, fixed)
../../pyomo/pyomo/pyomo/core/base/external.py:652: in __init__
    self.ra = (c_double*nr)(*_reals)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
self = 0.120272355044943*x6*(0.6468651916804*x0*x5/x1 + 0.6468651916804*x2*x3/x4)/x7 - 1.0

    def __float__(self):
        # Don't bother testing if it's a number; if it's not this is going
        # to fail, and if it is we still need to check that it evalf'ed to
        # a number.
        result = self.evalf()
        if result.is_Number:
            return float(result)
        if result.is_number and result.as_real_imag()[1]:
            raise TypeError("Cannot convert complex to float")

>       raise TypeError("Cannot convert expression to float")

E       TypeError: Cannot convert expression to float

../../anaconda3/lib/python3.9/site-packages/sympy/core/expr.py:345: TypeError

[andrewlee94]

Looking a bit deeper into the stack trace, it appears to me that the issue is arising when the Newton solver in calculate_variable_from_constraint attempts to get the derivative of a constraint w.r.t. to a variable in some situations.