incidence_analysis: Blank AssertionError with no error message showing

[ArkoprabhoDasgupta]

Summary

Addition of an error message in traceback for scc_solver

Rationale

Description

`File c:\Users\ad00105\AppData\Local\miniconda3\envs\working_env\Lib\site-packages\pyomo\contrib\incidence_analysis\scc_solver.py:69, in generate_strongly_connected_components(constraints, variables, include_fixed, igraph) 60 if variables is None: 61 variables = list( 62 _generate_variables_in_constraints( 63 constraints, (...) 66 ) 67 ) ---> 69 assert len(variables) == len(constraints) 70 if igraph is None: 71 igraph = IncidenceGraphInterface()

AssertionError: `

I was trying to initialize a solvent extraction unit operation model. During the initialization process, I was passing a model with non-zero degrees of freedom to the scc_solver, and I got a traceback. In the traceback, it is mentioned that there is an AssertionError, but there is no error message present after that, which explains the cause of the error. It would be helpful, if there is an error message shown, after the AssertionError traceback, so that it is easier to troubleshoot the code.

Additional information

[Robbybp]

Can you post some code I can use to reproduce this? If so, I can probably fix it.

[ArkoprabhoDasgupta]

from pyomo.environ import ConcreteModel, SolverFactory, Var, log
import numpy as np
from idaes.core.util.model_statistics import degrees_of_freedom

Elements = ["H", "HSO4", "SO4"]

C_r_ini = np.array([0.1, 0.1, 0])
C_r_ini_e = dict(zip(Elements, C_r_ini))

m = ConcreteModel()

m.C_r = Var(Elements, bounds=(1e-16, 100))
m.alpha = Var(bounds=(0, 100))

@m.Constraint()
def H_balance(m):
    return m.C_r["H"] == C_r_ini_e["H"] + m.alpha 

# @m.Constraint()
# def HSO4_balance(m):
#     return m.C_r["HSO4"] == C_r_ini_e["HSO4"] - m.alpha

# @m.Constraint()
# def SO4_balance(m):
#     return m.C_r["SO4"] == C_r_ini_e["SO4"] + m.alpha

@m.Constraint()
def kinetic_eqm(m):
    return log(m.C_r["H"]) + log(m.C_r["SO4"]) == log(m.C_r["HSO4"]) - 1.99 * log(10)

assert degrees_of_freedom(m) == 0

solver = SolverFactory("ipopt")
solver.solve(m, tee=True)

[ArkoprabhoDasgupta]

This is a small code I developed, if you run this you should see an assertion error with no error message, I hope this will help

[Robbybp]

Running this example, I see an assertion error due to your line: assert degrees_of_freedom(m) == 0. I believe you mentioned above that the assertion error you're seeing is in scc_solver?

[ArkoprabhoDasgupta]

Yes

[Robbybp]

The code you posted doesn't use solve_strongly_connected_components. Do you have an example that does?

[ArkoprabhoDasgupta]

from pyomo.environ import ConcreteModel, SolverFactory, Var, log
import numpy as np
from idaes.core.util.model_statistics import degrees_of_freedom
from pyomo.contrib.incidence_analysis import solve_strongly_connected_components

Elements = ["H", "HSO4", "SO4"]

C_r_ini = np.array([0.1, 0.1, 0])
C_r_ini_e = dict(zip(Elements, C_r_ini))

m = ConcreteModel()

m.C_r = Var(Elements, bounds=(1e-16, 100))
m.alpha = Var(bounds=(0, 100))

@m.Constraint()
def H_balance(m):
    return m.C_r["H"] == C_r_ini_e["H"] + m.alpha 

# @m.Constraint()
# def HSO4_balance(m):
#     return m.C_r["HSO4"] == C_r_ini_e["HSO4"] - m.alpha

# @m.Constraint()
# def SO4_balance(m):
#     return m.C_r["SO4"] == C_r_ini_e["SO4"] + m.alpha

@m.Constraint()
def kinetic_eqm(m):
    return log(m.C_r["H"]) + log(m.C_r["SO4"]) == log(m.C_r["HSO4"]) - 1.99 * log(10)

solver = SolverFactory("ipopt")
solve_strongly_connected_components(m, solver=solver)

[ArkoprabhoDasgupta]

I applied solve_strongly_connected_components to this same model and I saw the same error