How to properly switch the optimizer in IDAES from IPOPT to others (e.g. Couenne)?

[Hanli131]

Thank you in advance！💐

Issue Summary:
In an optimization problem with a clear and well-defined objective, I tried switching the solver from the default IPOPT to Couenne. Although the solver returns a result, the result is clearly incorrect. I would like to ask how to correctly switch the optimizer within the platform.

Detailed Description:
I created a simple example, and the system diagram is shown below.

The optimization objective is to maximize the mixer's outlet temperature. Intuitively, since the inlet temperature is fixed, increasing the flow rate of waterfeed2 should increase the outlet temperature. Therefore, the optimal flow rate should be the upper bound, 150 mol/s. Using IPOPT, I can obtain the correct result as shown below:

However, when I switch to Couenne, the result is as follows, and the flow rate remains at 90 mol/s, which is my initial value:

I am using Windows 10 as the platform, and the specific model code is provided below.
post.py.zip

import pyomo.environ as pyo
from pyomo.environ import Objective, maximize
from pyomo.environ import Constraint, value
from idaes.core import FlowsheetBlock
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.models.properties.modular_properties.base.generic_property import (
    GenericParameterBlock,)
from idaes.models_extra.power_generation.properties.natural_gas_PR import (
    get_prop,
    EosType)
from idaes.models.unit_models import Mixer, MomentumMixingType, Feed
from idaes.models.unit_models.heater import Heater
import idaes.logger as idaeslog

from pyomo.network import Arc, SequentialDecomposition

fuel_comps = ["H2", "H2O"]
oxygen_comps = ["O2", "N2"]

m = pyo.ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False, time_set=[0], time_units=pyo.units.s)   

m.fs.oxygen_params = GenericParameterBlock(
    **get_prop(oxygen_comps, {"Vap"}, eos=EosType.IDEAL), doc="Air-side parameters"
)
m.fs.fuel_params = GenericParameterBlock(
    **get_prop(fuel_comps, {"Vap"}, eos=EosType.IDEAL), doc="Fuel-side parameters"
)

#add units
m.fs.waterfeed = Feed(property_package=m.fs.fuel_params)

m.fs.waterfeed2 = Feed(property_package=m.fs.fuel_params)

m.fs.waterheater = Heater(property_package=m.fs.fuel_params)

m.fs.H2feed = Feed(property_package=m.fs.fuel_params)

m.fs.H2heater = Heater(property_package=m.fs.fuel_params)

m.fs.fuelmixer = Mixer(property_package=m.fs.fuel_params,
                     num_inlets=3,
                     momentum_mixing_type=MomentumMixingType.minimize)


#add arcs
m.fs.wf2wh = Arc(source=m.fs.waterfeed.outlet, destination=m.fs.waterheater.inlet)

m.fs.wh2fm = Arc(source=m.fs.waterheater.outlet, destination=m.fs.fuelmixer.inlet_1)

m.fs.hf2hh = Arc(source=m.fs.H2feed.outlet, destination=m.fs.H2heater.inlet)

m.fs.hh2fm = Arc(source=m.fs.H2heater.outlet, destination=m.fs.fuelmixer.inlet_2)

m.fs.wf22fm = Arc(source=m.fs.waterfeed2.outlet, destination=m.fs.fuelmixer.inlet_3)

pyo.TransformationFactory("network.expand_arcs").apply_to(m)


#set initial inputs
def flow_mol_rule(b):
    return m.fs.H2feed.flow_mol[0] == m.fs.waterfeed.flow_mol[0]/10

m.fs.H2feed.flow_mol_con = Constraint(rule=flow_mol_rule)

m.fs.waterfeed.flow_mol[0].fix(100)
m.fs.waterfeed.pressure[0].fix(101325)
m.fs.waterfeed.temperature[0].fix(400)
m.fs.waterfeed.mole_frac_comp[0, "H2"].fix(1e-6)    
m.fs.waterfeed.mole_frac_comp[0, "H2O"].fix(1-1e-6)

m.fs.waterfeed2.flow_mol[0].fix(90)
m.fs.waterfeed2.pressure[0].fix(101325)
m.fs.waterfeed2.temperature[0].fix(1000)
m.fs.waterfeed2.mole_frac_comp[0, "H2"].fix(1e-6)    
m.fs.waterfeed2.mole_frac_comp[0, "H2O"].fix(1-1e-6)

m.fs.H2feed.pressure.fix(101325)
m.fs.H2feed.temperature[0].fix(400)
m.fs.H2feed.mole_frac_comp[0, "H2"].fix(1-1e-6)    
m.fs.H2feed.mole_frac_comp[0, "H2O"].fix(1e-6)

m.fs.waterheater.outlet.temperature.fix(550)
def H2_heater_rule(b):
    return m.fs.H2heater.outlet.temperature[0] == m.fs.waterheater.outlet.temperature[0] - 100

m.fs.H2heater.H2_heater_con = Constraint(rule=H2_heater_rule)


print(f"initial DOF is {degrees_of_freedom(m)}")
print("before the initialization:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()


#initialization
seq = SequentialDecomposition()
seq.options.select_tear_method = "heuristic"
seq.options.tear_method = "Wegstein"
seq.options.iterLim = 500

G = seq.create_graph(m)
heuristic_tear_set = seq.tear_set_arcs(G, method="heuristic")
order = seq.calculation_order(G)
   
for o in heuristic_tear_set:
    print(o.name) 
   
for o in order:
    print(o[0].name)
    
    
def function(unit):
        unit.initialize(outlvl=idaeslog.INFO)
 
seq.run(m, function)

print("after the initialization:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()


#solve
#solver = pyo.SolverFactory('Couenne')
solver = pyo.SolverFactory('ipopt')
results = solver.solve(
    m, 
    tee=False,
    options={
        "tol": 1e-8,
        "max_iter": 20000,
        "nlp_scaling_method": "user-scaling",
        "constr_viol_tol": 1e-6,
    })  

print("calculation results:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()


#optimization
m.fs.waterfeed2.flow_mol[0].unfix()
m.fs.waterfeed2.flow_mol[0].setlb(50)
m.fs.waterfeed2.flow_mol[0].setub(150)

def mixer_outT_rule(m):
    return m.fs.fuelmixer.outlet.temperature[0]

m.objective = Objective(rule=mixer_outT_rule, sense=maximize)
#solver = pyo.SolverFactory('Couenne')
solver = pyo.SolverFactory('ipopt')
results = solver.solve(
    m, 
    tee=False)

print("optimization results:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()

[dallan-keylogic]

This is more of a Pyomo question than an IDAES question. I suggest you ask in their repository.