Handle bounds in `solve_strongly_connected_components`

[Robbybp]

Summary

The algorithm in solve_strongly_connected_components is not robust to the following situation:

1. A subsystem has multiple solutions. We converge to one of them.
2. For this particular local solution (but not some other local solution), a downstream subsystem is infeasible.

We would like the first subsystem to converge to the "right" local solution.

Rationale

Inspired by some problems @dallan-keylogic is having. See #3571.

Description

I propose we do this with FBBT. Here's a sketch of what this looks like in solve_strongly_connected_components:

    scc_input_list = list(
        generate_strongly_connected_components(constraints, variables, igraph=igraph)
    )

    if propagate_bounds:
        orig_bounds = [(var.lb, var.ub) for var in variables]
        for scc, inputs in scc_input_list:
            # Propagate bounds forward. This computes bounds on decision variables
            # based on fixed "degree of freedom" values. This may help us propagate
            # bound backwards in the reverse pass.
            fbbt(scc)
        for scc, inputs in reversed(scc_input_list):
            # Propagate bounds in the "reverse direction". This attempts to strengthen
            # "upstream" bounds based on existing bounds on decision variables.
            fbbt(scc)

    for scc, inputs in scc_input_list:
        # Solve scc...

    if propagate_bounds:
        # Reset bounds to their original values
        for var, (lb, ub) in zip(variables, orig_bounds):
            var.setlb(lb)
            var.setub(ub)

I've implemented this in my scc-fbbt branch. @dallan-keylogic, can you check if this branch helps with your infeasibility issue? Here's an example:

import pyomo.environ as pyo
from pyomo.contrib.incidence_analysis import solve_strongly_connected_components

m = pyo.ConcreteModel()
m.x = pyo.Var([1, 2, 3, 4], initialize=1.0)
m.eq = pyo.Constraint(pyo.PositiveIntegers)
m.eq[1] = m.x[1]**2 == 1.0
m.x[1] = -2.0
m.eq[2] = m.x[2]**2 == 4.0
m.x[2].setlb(1.0)
m.eq[3] = m.x[2] == m.x[3]
m.eq[4] = m.x[1] * m.x[3] * m.x[4] == 1.0
m.x[4].setlb(0.0)
solver = pyo.SolverFactory("ipopt")
solve_strongly_connected_components(m, solver=solver, use_calc_var=False, propagate_bounds=True)
m.x.pprint()

With propagate_bounds=False, this prints:

WARNING: Loading a SolverResults object with a warning status into
model.name="ScalarBlock";
    - termination condition: infeasible
    - message from solver: Ipopt 3.14.11: Converged to a locally infeasible
      point. Problem may be infeasible.
x : Size=4, Index={1, 2, 3, 4}
    Key : Lower : Value                  : Upper : Fixed : Stale : Domain
      1 :  None :     -1.000000000000001 :  None : False :  True :  Reals
      2 :   1.0 :     2.0000000000000013 :  None : False :  True :  Reals
      3 :  None :     2.0000000000000013 :  None : False :  True :  Reals
      4 :   0.0 : -9.958130312559555e-09 :  None : False : False :  Reals

With propagate_bounds=True:

x : Size=4, Index={1, 2, 3, 4}
    Key : Lower : Value : Upper : Fixed : Stale : Domain
      1 :  None :   1.0 :  None : False :  True :  Reals
      2 :   1.0 :   2.0 :  None : False :  True :  Reals
      3 :  None :   2.0 :  None : False :  True :  Reals
      4 :   0.0 :   0.5 :  None : False : False :  Reals

I'm not sure how robust this is in general — maybe fbbt(scc) is not guaranteed to converge? You could probably find a case where this fails to restrict an infeasible solution, but it seems like it should work for simple systems and small SCCs.

Note that this only makes sense if the "subsolver" can solve square problems with bounds, which rules out calculate_variable_from_constraint. Luckily we can turn this off with use_calc_var=False. There's also the question of what to do about inequalities. I think we can just include them when they contain any variable in a particular subsystem. We could also punt on these until they come up.

Another approach would be to use a global solver, enumerate all solutions of each feasibility problem, and build a tree of possible solutions. This is probably out-of-scope.

Additional info

A larger issue here is that solve_strongly_connected_components isn't really a solver, it's just a utility function. It does very little to try to be consistent when different subsolvers are used and doesn't return a results object. If we ever make it a real solver, we should be explicit and consistent about whether bounds and inequalities are considered.

[dallan-keylogic]

@Robbybp As of yet, convergence to a wrong solution is hypothetical. I became aware of the issue only because the cubic complementarity VLE formulation was being solved in an order where the nonnegativity constraints inherent in the formulation were downstream of the main block in which most of the properties were being solved for. If that block had converged to a point for which we could not solve the complementarity constraints, we'd be out of luck.

[Robbybp]

If that block had converged to a point for which we could not solve the complementarity constraints, we'd be out of luck.

I agree. My conjecture is that using the scc-fbbt branch with propagate_bounds=True (and use_calc_var=False) will prevent this potential issue. If you could hack an initial point where this issue becomes a reality, that would be really useful to test this new feature.

[michaelbynum]

The implementation described above makes sense to me and certainly may help in some cases. It is definitely true that fbbt may not converge. In the example below, y should clearly be zero, but fbbt cannot determine that (unless you simplify the expression first). However, there are many examples, where it would work just fine.

In [1]: import pyomo.environ as pe
In [2]: m = pe.ConcreteModel()
In [3]: m.x = pe.Var()
In [4]: m.y = pe.Var()
In [5]: m.c = pe.Constraint(expr=m.x - m.x + m.y == 0)
In [6]: from pyomo.contrib.fbbt.fbbt import fbbt
In [7]: _ = fbbt(m.c)
In [8]: m.x.bounds
Out[8]: (None, None)
In [9]: m.y.bounds
Out[9]: (None, None)