search

Pyomo / pyomo

An object-oriented algebraic modeling language in Python for structured optimization problems.
https://www.pyomo.org
Other
1.98k stars 509 forks source link

FBBT handling of Expr_if expressions #3083

Closed kurbansitterley closed 7 months ago

kurbansitterley commented 8 months ago

Summary

I am a developer working on WaterTAP and a subpackage WaterTAP-REFLO.

We have a model that costs a MgCl2 flow. The constraint used to cost the flow uses Expr_if statements and I am guessing these are causing the DeveloperError to be raised(?). The error is raised during the build.

As of now we are excluding this flow from our costing package.

Steps to reproduce the issue

  1. install watertap-reflo environment
  2. uncomment line 972 in src > watertap_contrib > reflo > costing > units > chemical_softening_zo.py
  3. Run the following
from pyomo.environ import (
    ConcreteModel,
    value,
    assert_optimal_termination,
    units as pyunits,
)

from idaes.core import FlowsheetBlock, UnitModelCostingBlock, MaterialFlowBasis
from idaes.core.util.scaling import calculate_scaling_factors
from idaes.core.solvers import get_solver

from watertap_contrib.reflo.unit_models.zero_order.chemical_softening_zo import (
    ChemicalSofteningZO,
)
from watertap_contrib.reflo.costing import TreatmentCosting
from watertap.property_models.multicomp_aq_sol_prop_pack import (
    MCASParameterBlock,
)

solver = get_solver()

component_list = ["Ca_2+", "Mg_2+", "SiO2", "Alkalinity_2-"]

m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)
m.fs.properties = MCASParameterBlock(
    solute_list=component_list, material_flow_basis=MaterialFlowBasis.mass
)

m.fs.soft = soft = ChemicalSofteningZO(
    property_package=m.fs.properties,
    silica_removal=True,
    softening_procedure_type="excess_lime_soda",
)

ca_in = 1.43 * pyunits.kg / pyunits.m**3  # g/L = kg/m3
mg_in = 0.1814 * pyunits.kg / pyunits.m**3  # g/L = kg/m3
sio2_in = 0.054 * pyunits.kg / pyunits.m**3  # g/L = kg/m3
alk_in = 0.421 * pyunits.kg / pyunits.m**3  # g/L = kg/m3
CO2_in = 0.10844915 * pyunits.kg / pyunits.m**3
q_in = 3785 * pyunits.m**3 / pyunits.day  # m3/d
rho = 1000 * pyunits.kg / pyunits.m**3
prop_in = soft.properties_in[0]

flow_mass_phase_water = pyunits.convert(q_in * rho, to_units=pyunits.kg / pyunits.s)
flow_mass_phase_ca = pyunits.convert(q_in * ca_in, to_units=pyunits.kg / pyunits.s)
flow_mass_phase_mg = pyunits.convert(q_in * mg_in, to_units=pyunits.kg / pyunits.s)
flow_mass_phase_si = pyunits.convert(q_in * sio2_in, to_units=pyunits.kg / pyunits.s)
flow_mass_phase_alk = pyunits.convert(q_in * alk_in, to_units=pyunits.kg / pyunits.s)
prop_in.flow_mass_phase_comp["Liq", "H2O"].fix(flow_mass_phase_water())

prop_in.flow_mass_phase_comp["Liq", "Ca_2+"].fix(value(flow_mass_phase_ca))
prop_in.flow_mass_phase_comp["Liq", "Mg_2+"].fix(value(flow_mass_phase_mg))
prop_in.flow_mass_phase_comp["Liq", "SiO2"].fix(value(flow_mass_phase_si))
prop_in.flow_mass_phase_comp["Liq", "Alkalinity_2-"].fix(value(flow_mass_phase_alk))
prop_in.temperature.fix()
prop_in.pressure.fix()

soft.ca_eff_target.fix(3)
soft.mg_eff_target.fix(0.2)

soft.no_of_mixer.fix(1)
soft.no_of_floc.fix(2)
soft.retention_time_mixer.fix(0.4)
soft.retention_time_floc.fix(25)
soft.retention_time_sed.fix(130)
soft.retention_time_recarb.fix(20)
soft.frac_vol_recovery.fix()
soft.removal_efficiency.fix()
soft.CO2_CaCO3.fix(CO2_in)
soft.vel_gradient_mix.fix(300)
soft.vel_gradient_floc.fix(50)

m.fs.properties.set_default_scaling(
    "flow_mass_phase_comp", 1 / flow_mass_phase_water(), index=("Liq", "H2O")
)
m.fs.properties.set_default_scaling(
    "flow_mass_phase_comp", 1 / flow_mass_phase_ca(), index=("Liq", "Ca_2+")
)
m.fs.properties.set_default_scaling(
    "flow_mass_phase_comp", 1 / flow_mass_phase_mg(), index=("Liq", "Mg_2+")
)
m.fs.properties.set_default_scaling(
    "flow_mass_phase_comp", 1 / flow_mass_phase_si(), index=("Liq", "SiO2")
)
m.fs.properties.set_default_scaling(
    "flow_mass_phase_comp",
    1 / flow_mass_phase_alk(),
    index=("Liq", "Alkalinity_2-"),
)

calculate_scaling_factors(m)
m.fs.costing = TreatmentCosting()
m.fs.soft.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing)
m.fs.costing.cost_process()
m.fs.costing.add_LCOW(prop_in.flow_vol)

results = solver.solve(m)
assert_optimal_termination(results)

Error Message

Traceback (most recent call last):
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/idaes/core/base/process_block.py", line 41, in _rule_default
    b.build()
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/idaes/core/base/costing_base.py", line 605, in build
    method(self, **self.config.costing_method_arguments)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/watertap/costing/util.py", line 86, in add_costing_parameter_block
    return func(blk, *args, **kwargs)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/watertap/costing/util.py", line 86, in add_costing_parameter_block
    return func(blk, *args, **kwargs)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/watertap/costing/util.py", line 86, in add_costing_parameter_block
    return func(blk, *args, **kwargs)
  [Previous line repeated 2 more times]
  File "/Users/ksitterl/Documents/Python/watertap-reflo/watertap-reflo/src/watertap_contrib/reflo/costing/units/chemical_softening_zo.py", line 972, in cost_chemical_softening
    blk.costing_package.cost_flow(blk.mgcl2_dosing, "mgcl2")
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/idaes/core/base/costing_base.py", line 283, in cost_flow
    ebounds = compute_bounds_on_expr(flow_expr)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/contrib/fbbt/fbbt.py", line 1506, in compute_bounds_on_expr
    ).walk_expression(expr)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/core/expr/visitor.py", line 268, in walk_expression
    result = self._process_node(root, RECURSION_LIMIT)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/core/expr/visitor.py", line 470, in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/core/expr/visitor.py", line 470, in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/core/expr/visitor.py", line 470, in _process_node_bx
    data.append(self._process_node(child, recursion_limit))
  [Previous line repeated 2 more times]
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/core/expr/visitor.py", line 481, in _process_node_bx
    return self.exitNode(node, data)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/contrib/fbbt/expression_bounds_walker.py", line 283, in exitNode
    return _operator_dispatcher[node.__class__](self, node, *data)
  File "/Users/ksitterl/opt/anaconda3/envs/watertap-reflo/lib/python3.10/site-packages/pyomo/repn/util.py", line 419, in register_dispatcher
    raise DeveloperError(
pyomo.common.errors.DeveloperError: Internal Pyomo implementation error:
        "Unexpected expression node type 'InequalityExpression' found while
        walking expression tree."
    Please report this to the Pyomo Developers.
ERROR: Constructing component 'fs.soft.costing' from data=None failed:
DeveloperError: Internal Pyomo implementation error:
            "Unexpected expression node type 'InequalityExpression' found
            while walking expression tree."
        Please report this to the Pyomo Developers.

Information on your system

Pyomo version: 6.7.0 Python version: 3.10.13 Operating system: MacOS Ventura 13.6.3 (22G436) How Pyomo was installed (PyPI, conda, source): conda Solver (if applicable): watertap-ipopt

Additional information

blnicho commented 8 months ago

Here is a MWE that reproduces the error:

import pyomo.environ as pe
from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr

m = pe.ConcreteModel()

m.x = pe.Var()
m.y = pe.Var()
e = pe.Expr_if(m.x <= 0, m.y + m.x == 0, m.y - m.x == 0)

ebounds = compute_bounds_on_expr(e)

Seems like the issue is that FBBT is not failing gracefully when it encounters unknown expression types like Expr_if. @michaelbynum this seems counter to your comment in #846.

michaelbynum commented 8 months ago

Thanks. I'll take a look.

michaelbynum commented 8 months ago

My comment in #846 was correct except that compute_bounds_on_expr now uses a different walker than fbbt. I just opened a PR to fix the new walker for compute_bounds_on_expr.