Trough surrogate model infeasible solves with WaterTAP solver

[kurbansitterley]

Description

The trough surrogate model does not consistently solve with solver = get_solver() but consistently solves with solver = SolverFactory("ipopt"). Sometimes the infeasibilities are with only the trough on the flowsheet and sometimes when the trough costing package is added.

Steps to Reproduce

The most straightforward way to reproduce is to replace the solver in the trough test file with the WaterTAP solver. The tests with and without costing will fail. The current trough test file does not have scaling, but if the following scaling is added to the test frame, the costing test will pass but the validation test will still fail:

set_scaling_factor(m.fs.trough.electricity, 1e-4)
set_scaling_factor(m.fs.trough.heat, 1e-5)
set_scaling_factor(m.fs.trough.heat_annual, 1e-9)
set_scaling_factor(m.fs.trough.electricity_annual, 1e-8)

Here is a specific build with LT-MED included on the flowsheet that will results in an infeasible solve:

from pyomo.environ import (
    ConcreteModel,
    assert_optimal_termination,
    units as pyunits,
    SolverFactory,
)

from idaes.core.solvers import get_solver
from idaes.core.util.scaling import *
from idaes.core import FlowsheetBlock, UnitModelCostingBlock
from idaes.core.util.model_statistics import *

from watertap.property_models.seawater_prop_pack import SeawaterParameterBlock
from watertap.property_models.water_prop_pack import WaterParameterBlock
from watertap_contrib.seto.costing import (
    EnergyCosting,
    SETOSystemCosting,
    TreatmentCosting,
)

from watertap_contrib.seto.solar_models.surrogate.trough import TroughSurrogate
from watertap_contrib.seto.unit_models.surrogate import LTMEDSurrogate
from watertap.core.util.model_diagnostics.infeasible import *

solver_wt = get_solver()
solver_ipopt = SolverFactory("ipopt")

trough_scale_cost_vars = ['fs.energy.trough.costing.capital_cost',
 'fs.energy.trough.costing.variable_operating_cost',
 'fs.energy.trough.costing.fixed_operating_cost',
 'fs.energy.trough.costing.direct_cost',
 'fs.energy.costing.aggregate_capital_cost',
 'fs.energy.costing.aggregate_fixed_operating_cost',
 'fs.energy.costing.aggregate_variable_operating_cost',
 'fs.energy.costing.aggregate_flow_electricity',
 'fs.energy.costing.aggregate_flow_heat',
 'fs.energy.costing.total_capital_cost',
 'fs.energy.costing.total_operating_cost',
 ]

def build_trough(m=None, add_trough_costing=False):
    if m is None:
        m = ConcreteModel()
        m.fs = FlowsheetBlock(dynamic=False)
    m.fs.energy = Block()
    m.fs.energy.trough = trough = TroughSurrogate()
    m.fs.energy.trough.heat_load.fix(500)
    m.fs.energy.trough.hours_storage.fix(12)
    if add_trough_costing:
        m.fs.energy.costing = EnergyCosting()
        m.fs.energy.trough.costing = UnitModelCostingBlock(
            flowsheet_costing_block=m.fs.energy.costing
        )
        m.fs.energy.costing.cost_process()
        for v in m.fs.energy.component_objects(Var):
            if v.name in trough_scale_cost_vars:
                set_scaling_factor(v, 1e-5)

    set_scaling_factor(trough.electricity, 1e-4)
    set_scaling_factor(trough.heat, 1e-5)
    set_scaling_factor(trough.heat_annual, 1e-9)
    set_scaling_factor(trough.electricity_annual, 1e-8)

    return m

def build_trough_ltmed(
    add_ltmed_costing=False,
    add_trough=False,
    add_trough_costing=False,
):
    m = ConcreteModel()
    m.fs = FlowsheetBlock(dynamic=False)
    m.fs.water_prop = SeawaterParameterBlock()
    m.fs.steam_prop = WaterParameterBlock()
    m.fs.treatment = Block()
    m.fs.treatment.lt_med = lt_med = LTMEDSurrogate(
        property_package_liquid=m.fs.water_prop,
        property_package_vapor=m.fs.steam_prop,
    )

    dist = lt_med.distillate_props[0]
    steam = lt_med.steam_props[0]

    # LT-MED system specification
    lt_med.recovery_vol_phase[0, "Liq"].fix(0.5)
    feed_flow = pyunits.convert(
        50 * pyunits.Mgallons / pyunits.day, to_units=pyunits.m**3 / pyunits.s
    )
    lt_med.feed_props.calculate_state(
        var_args={
            ("flow_vol_phase", "Liq"): feed_flow,
            ("conc_mass_phase_comp", ("Liq", "TDS")): 60,
            ("temperature", None): 25 + 273.15,
            ("pressure", None): 101325,
        },
        hold_state=True,
    )

    steam.temperature.fix(80 + 273.15)

    #  Add LT-MED costing
    if add_ltmed_costing:
        m.fs.treatment.costing = TreatmentCosting()
        m.fs.treatment.lt_med.costing = UnitModelCostingBlock(
            flowsheet_costing_block=m.fs.treatment.costing
        )
        m.fs.treatment.costing.cost_process()
        m.fs.treatment.costing.add_LCOW(dist.flow_vol_phase["Liq"])

    # Build energy block
    if add_trough:
        m = build_trough(
            m=m, add_trough_costing=add_trough_costing
        )
        m.fs.energy.trough.heat_load.unfix()

        if add_trough_costing:
            # Add system costing
            m.fs.sys_costing = SETOSystemCosting()
            m.fs.sys_costing.add_LCOW(dist.flow_vol_phase["Liq"])

        # Linking constraint
        m.fs.energy.lt_med_heat_demand_constr = Constraint(
            expr=m.fs.energy.trough.heat_load
            >= pyunits.convert(lt_med.thermal_power_requirement, to_units=pyunits.MW)
        )

    m.fs.water_prop.set_default_scaling(
        "flow_mass_phase_comp", 1e-3, index=("Liq", "H2O")
    )
    m.fs.water_prop.set_default_scaling(
        "flow_mass_phase_comp", 1e-2, index=("Liq", "TDS")
    )
    m.fs.steam_prop.set_default_scaling(
        "flow_mass_phase_comp", 1e-3, index=("Liq", "H2O")
    )
    m.fs.steam_prop.set_default_scaling(
        "flow_mass_phase_comp", 1e-3, index=("Vap", "H2O")
    )

    calculate_scaling_factors(m)
    lt_med.initialize()

    return m

'''
Only Trough
Optimal with WaterTAP solver
'''
print("Trough no costing")
m = build_trough()
results = solver_wt.solve(m)
assert_optimal_termination(results)

'''
Fails with WaterTAP solver, 
unless use calculate_variable_from_constraint prior to solve
'''
print("Trough with costing")
m = build_trough(add_trough_costing=True)
# calculate_variable_from_constraint(
#     m.fs.energy.trough.heat_annual,
#     m.fs.energy.trough.surrogate_blk.pysmo_constraint["heat_annual"],
# )
# m.fs.energy.trough.heat_annual.fix()
# m.fs.energy.trough.heat_load.unfix()
# calculate_variable_from_constraint(
#     m.fs.energy.trough.electricity_annual,
#     m.fs.energy.trough.surrogate_blk.pysmo_constraint["electricity_annual"],
# )
# results = solver_wt.solve(m)
results = solver_ipopt.solve(m)
assert_optimal_termination(results)

'''
Flowsheet with only LT-MED
Optimal with WaterTAP solver
'''
print("LT-MED no costing")
m = build_trough_ltmed()
results = solver_wt.solve(m)
assert_optimal_termination(results)

'''
Flowsheet with only LT-MED + costing
Optimal with WaterTAP solver
'''
print("LT-MED with costing")
m = build_trough_ltmed(add_ltmed_costing=True)
results = solver_wt.solve(m)
assert_optimal_termination(results)

'''
Flowsheet with LT-MED + Trough
No trough costing
Optimal with WaterTAP solver
'''
print("LT-MED + Trough, no Trough costing")
m = build_trough_ltmed(add_ltmed_costing=True, add_trough=True)
results = solver_wt.solve(m)
assert_optimal_termination(results)

'''
Flowsheet with LT-MED + Trough
Add trough costing
Optimal with WaterTAP solver
'''
print("LT-MED + Trough, with Trough costing")
m = build_trough_ltmed(
    add_ltmed_costing=True,
    add_trough=True,
    add_trough_costing=True,
)
results = solver_wt.solve(m)
assert_optimal_termination(results)

In this demo I have also scaled the costing variables for the trough. If they aren't scaled, the last flowsheet (LT-MED + Trough with costing) will be infeasible.

For the second flowsheet, using the commented code with calculate_variable_from_constraint() will lead to an optimal solve.

[kurbansitterley]

@bknueven interested in any thoughts you have on this issue.

[bknueven]

So, using the WaterTAP solver applies variable scaling by default, whereas get_solver("ipopt") will not apply variable scaling.

There are also a few other differences I do not want to enumerate here.

@kurbansitterley if you can share the ipopt logs with me either here or elsewhere that would be helpful.

[kurbansitterley]

Here is the ipopt log for the failure with the WaterTAP solver:

ipopt-watertap: Ipopt with user variable scaling and IDAES jacobian constraint scaling
Ipopt 3.13.2: halt_on_ampl_error=yes
tol=1e-08
constr_viol_tol=1e-08
nlp_scaling_method=user-scaling
bound_relax_factor=0.0
alpha_for_y=bound-mult

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk/.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:       36
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:       18
                     variables with only lower bounds:        3
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       18
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 2.61e+09 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 4.77e-07 1.00e+00  -1.0 3.07e+07    -  9.90e-07 1.00e+00h  1
   2  0.0000000e+00 9.28e-08 1.00e-02  -1.0 7.06e-09    -  9.90e-01 1.00e+00   0
   3  0.0000000e+00 5.96e-08 1.01e-04  -1.0 1.74e-09    -  9.90e-01 1.00e+00   0
   4  0.0000000e+00 9.28e-08 2.00e-06  -1.0 1.74e-09    -  9.90e-01 1.00e+00   0
   5  0.0000000e+00 5.96e-08 1.01e-06  -1.0 1.74e-09    -  9.93e-01 1.00e+00T  0
   6  0.0000000e+00 9.28e-08 2.00e-07  -1.7 1.74e-09    -  1.00e+00 1.00e+00T  0
   7  0.0000000e+00 5.96e-08 1.50e-09  -3.8 1.74e-09    -  1.00e+00 1.00e+00T  0
   8  0.0000000e+00 9.28e-08 1.84e-11  -5.7 1.74e-09    -  1.00e+00 1.00e+00T  0
   9  0.0000000e+00 5.96e-08 2.51e-14  -8.6 1.74e-09    -  1.00e+00 1.00e+00T  0
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  0.0000000e+00 9.28e-08 9.09e-15  -9.0 1.74e-09    -  1.00e+00 1.00e+00T  0

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   9.0909090909253664e-15    9.0909090909090953e-15
Constraint violation....:   1.2417634328206383e-10    9.2786489744867137e-08
Complementarity.........:   9.0909090909090900e-10    9.0909090909090900e-10
Overall NLP error.......:   1.2417634328206383e-10    9.2786489744867137e-08

Number of objective function evaluations             = 11
Number of objective gradient evaluations             = 11
Number of equality constraint evaluations            = 11
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 11
Total CPU secs in IPOPT (w[/o](https://file+.vscode-resource.vscode-cdn.net/o) function evaluations)   =      0.001
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Search Direction is becoming Too Small.

However, @MuktaHardikar is able to run her treatment + energy_surrogate flowsheets with the WaterTAP solver, so I suspect it is something I did wrong with my build... However, we are still not able to run the tests of the surrogates with the WaterTAP solver

[bknueven]

Thanks @kurbansitterley; this log is helpful.

I think the issue here is that your problem seems to have large orders of magnitude (based on the variable scaling factors), whereas the watertap-ipopt solver is tuned by default to handle models with small magnitudes.

You could use the constraint_scaling_transform to scale the constraints that don't need to be enforced to 1e-08. But I'm going to revisit some of the assumptions ipopt-watertap makes based on this example.

[adam-a-a]

I would also add that you should probably be initializing your costing block if you aren't already. I did a quick ctrl+f on your pasted code, @kurbansitterley , so may have missed that if you do include initialization for costing as well.

[adam-a-a]

@kurbansitterley and team - I copied your script and was able to get the solve to work to some extent by commenting out the following lines:

        # for v in m.fs.energy.component_objects(Var):
        #     if v.name in trough_scale_cost_vars:
        #         set_scaling_factor(v, 1e-5)

    # set_scaling_factor(trough.electricity, 1e-4)
    # set_scaling_factor(trough.heat, 1e-5)
    # set_scaling_factor(trough.heat_annual, 1e-9)
    # set_scaling_factor(trough.electricity_annual, 1e-8)

Now, the first WT solve that failed will pass, and I make it up to the follow section before the WT solve fails again. ''' Flowsheet with LT-MED + Trough No trough costing Optimal with WaterTAP solver ''' Will continue probing ...

[kurbansitterley]

@adam-a-a your trial here is sort of illustrative of my point: ideally, for a given surrogate under given conditions, you shouldn't have to selectively scale some variables and not scale others based on how it is used in a flowsheet (right?). Like I don't scale certain variables for RO model depending on if I am using it with or without costing.

Looking at this now in light of our conversation today, I think an initialization routine that makes use of calculate_variable_from_constraint should help here. In the issue, I even mention that it will solve with get_solver() if you do that:

Fails with WaterTAP solver, 
unless use calculate_variable_from_constraint prior to solve
'''
print("Trough with costing")
m = build_trough(add_trough_costing=True)
# calculate_variable_from_constraint(
#     m.fs.energy.trough.heat_annual,
#     m.fs.energy.trough.surrogate_blk.pysmo_constraint["heat_annual"],
# )
# m.fs.energy.trough.heat_annual.fix()
# m.fs.energy.trough.heat_load.unfix()
# calculate_variable_from_constraint(
#     m.fs.energy.trough.electricity_annual,
#     m.fs.energy.trough.surrogate_blk.pysmo_constraint["electricity_annual"],
# )
# results = solver_wt.solve(m)

So that is something that I will try. But I appreciate any additional effort you put into this.

[adam-a-a]

Definitely agree with your points- I probed a bit more and the trouble seems to start once including trough and ltmed, without costing enabled on either model. I'm interested and want to help, so I'll check a couple more things out and report back with any findings (although you'll probably figure out a solution to the issue quickly).

[adam-a-a]

62 should also be relevant to this issue. We should confirm whether this issue is resolved by #62 before closing. From a very quick skim through the aforementioned PR, I only noticed usage of "ipopt" rather than "ipopt-watertap".

[adam-a-a]

This was resolved by PR #62.