Skip to content

Commit

Permalink
Refine some scaling in demo flowsheet
Browse files Browse the repository at this point in the history
  • Loading branch information
@MarcusHolly
MarcusHolly committed Sep 8, 2023
1 parent 5abe575 commit afde76a
Show file tree
Hide file tree
Showing 2 changed files with 87 additions and 43 deletions.
121 changes: 79 additions & 42 deletions idaes/models/flowsheets/demo_flowsheet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
Constructs a basic flowsheet for a benzene-toluene system with a mixer, heater
and flash unit.
"""

import pyomo.environ as pyo
from pyomo.environ import ConcreteModel, TransformationFactory
from pyomo.network import Arc

Expand All @@ -31,6 +31,10 @@
from idaes.models.unit_models import Mixer, Heater, Flash

import idaes.logger as idaeslog
from idaes.core.util.model_diagnostics import DiagnosticsToolbox
from watertap.core.util.model_diagnostics.infeasible import *
from idaes.core.util.model_statistics import degrees_of_freedom, total_constraints_set, large_residuals_set, variables_near_bounds_generator



def build_flowsheet():
Expand All @@ -57,15 +61,21 @@ def build_flowsheet():

def set_scaling(m):
"""Set scaling for demo flowsheet"""
m.fs.BT_props.set_default_scaling(
"flow_mol_phase_comp", 1e2, index=("Liq", "benzene")
)
m.fs.BT_props.set_default_scaling(
"flow_mol_phase_comp", 1e2, index=("Liq", "toluene")
)

iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_comp["benzene"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_comp["toluene"], 1e4)
# m.fs.BT_props.set_default_scaling(
# "flow_mol_phase_comp", 1e2, index=("Liq", "benzene")
# )
# m.fs.BT_props.set_default_scaling(
# "flow_mol_phase_comp", 1e2, index=("Liq", "toluene")
# )

# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_comp["benzene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_comp["toluene"], 1e5)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].flow_mol_phase["Liq"], 1e6)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].flow_mol_phase["Vap"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Liq", "toluene"], 1e5)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].mole_frac_phase_comp["Vap", "toluene"], 1e6)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].pressure_sat_comp["benzene"], 1e-5)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].pressure_sat_comp["toluene"], 1e-4)
iscale.set_scaling_factor(
m.fs.M01.inlet_1_state[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4
)
Expand All @@ -78,21 +88,29 @@ def set_scaling(m):
iscale.set_scaling_factor(
m.fs.M01.inlet_1_state[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4
)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_total, 1)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_sum_mol_frac, 1)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_comp["benzene"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_comp["toluene"], 1)
iscale.set_scaling_factor(
m.fs.M01.inlet_1_state[0].eq_phase_equilibrium["benzene"], 1
)
iscale.set_scaling_factor(
m.fs.M01.inlet_1_state[0].eq_phase_equilibrium["toluene"], 1
)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_P_vap["benzene"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_P_vap["toluene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_total, 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_sum_mol_frac, 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_comp["benzene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_comp["toluene"], 1)
# iscale.set_scaling_factor(
# m.fs.M01.inlet_1_state[0].eq_phase_equilibrium["benzene"], 1
# )
# iscale.set_scaling_factor(
# m.fs.M01.inlet_1_state[0].eq_phase_equilibrium["toluene"], 1
# )
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_P_vap["benzene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_1_state[0].eq_P_vap["toluene"], 1)
#
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_comp["benzene"], 1e5)
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_comp["toluene"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].flow_mol_phase["Vap"], 1e1)

iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Liq", "benzene"], 1e5)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_phase_comp["Vap", "benzene"], 1e5)

iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].pressure_sat_comp["benzene"], 1e-6)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].pressure_sat_comp["toluene"], 1e-6)

iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_comp["benzene"], 1e4)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].mole_frac_comp["toluene"], 1)
iscale.set_scaling_factor(
m.fs.M01.inlet_2_state[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4
)
Expand All @@ -105,18 +123,18 @@ def set_scaling(m):
iscale.set_scaling_factor(
m.fs.M01.inlet_2_state[0].enth_mol_phase_comp["Vap", "toluene"], 1e-4
)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_total, 1)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_sum_mol_frac, 1)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_comp["benzene"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_comp["toluene"], 1)
iscale.set_scaling_factor(
m.fs.M01.inlet_2_state[0].eq_phase_equilibrium["benzene"], 1
)
iscale.set_scaling_factor(
m.fs.M01.inlet_2_state[0].eq_phase_equilibrium["toluene"], 1
)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_P_vap["benzene"], 1)
iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_P_vap["toluene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_total, 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_sum_mol_frac, 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_comp["benzene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_comp["toluene"], 1)
# iscale.set_scaling_factor(
# m.fs.M01.inlet_2_state[0].eq_phase_equilibrium["benzene"], 1
# )
# iscale.set_scaling_factor(
# m.fs.M01.inlet_2_state[0].eq_phase_equilibrium["toluene"], 1
# )
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_P_vap["benzene"], 1)
# iscale.set_scaling_factor(m.fs.M01.inlet_2_state[0].eq_P_vap["toluene"], 1)

iscale.set_scaling_factor(
m.fs.M01.mixed_state[0].enth_mol_phase_comp["Liq", "benzene"], 1e-4
Expand Down Expand Up @@ -412,18 +430,18 @@ def set_scaling(m):

iscale.calculate_scaling_factors(m)


def set_dof(m):
"""Set degrees of freedom for demo flowsheet"""
eps = 1e-5
m.fs.M01.inlet_1.flow_mol.fix(1.0)
m.fs.M01.inlet_1.mole_frac_comp[:, "benzene"].fix(1)
m.fs.M01.inlet_1.mole_frac_comp[:, "toluene"].fix(1e-5)
m.fs.M01.inlet_1.mole_frac_comp[:, "benzene"].fix(1 - eps)
m.fs.M01.inlet_1.mole_frac_comp[:, "toluene"].fix(eps)
m.fs.M01.inlet_1.pressure.fix(101325)
m.fs.M01.inlet_1.temperature.fix(370)

m.fs.M01.inlet_2.flow_mol.fix(1.0)
m.fs.M01.inlet_2.mole_frac_comp[:, "benzene"].fix(1e-5)
m.fs.M01.inlet_2.mole_frac_comp[:, "toluene"].fix(1)
m.fs.M01.inlet_2.mole_frac_comp[:, "benzene"].fix(eps)
m.fs.M01.inlet_2.mole_frac_comp[:, "toluene"].fix(1 - eps)
m.fs.M01.inlet_2.pressure.fix(1.3e5)
m.fs.M01.inlet_2.temperature.fix(380)

Expand All @@ -442,12 +460,31 @@ def initialize_flowsheet(m):
propagate_state(m.fs.s02)
m.fs.F03.initialize(outlvl=idaeslog.WARNING)


def solve_flowsheet(m, stee=False):
"""Solve demo flowsheet"""
solver = get_solver()
solver.solve(m, tee=stee)

dt = DiagnosticsToolbox(model=m)
dt.report_numerical_issues()
dt.display_constraints_with_large_residuals()
dt.display_variables_at_or_outside_bounds()

badly_scaled_var_list = iscale.badly_scaled_var_generator(m, large=1e2, small=1e-2)
print("---------------- badly_scaled_var_list ----------------")
for x in badly_scaled_var_list:
print(f"{x[0].name}\t{x[0].value}\tsf: {iscale.get_scaling_factor(x[0])}")
print("---------------- variables_near_bounds_list ----------------")
variables_near_bounds_list = variables_near_bounds_generator(m)
for x in variables_near_bounds_list:
print(f"{x.name}\t{x.value}")
print("---------------- total_constraints_set_list ----------------")
total_constraints_set_list = total_constraints_set(m)
for x in total_constraints_set_list:
residual = abs(pyo.value(x.body) - pyo.value(x.lb))
if residual > 1e-8:
print(f"{x}\t{residual}")


def display_results(m):
"""Display results for flowsheet"""
Expand Down
9 changes: 8 additions & 1 deletion idaes/models/flowsheets/tests/test_demo_flowsheet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
from idaes.models.flowsheets.demo_flowsheet import (
build_flowsheet,
set_dof,
set_scaling,
initialize_flowsheet,
solve_flowsheet,
)
Expand All @@ -34,6 +35,7 @@
)
from idaes.models.unit_models import Mixer, Heater, Flash
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.util.model_diagnostics import DiagnosticsToolbox


@pytest.fixture(scope="module")
Expand All @@ -59,6 +61,11 @@ def test_build_flowsheet(model):
assert degrees_of_freedom(model) == 13


@pytest.mark.unit
def test_set_scaling(model):
set_scaling(model)


@pytest.mark.unit
def test_set_dof(model):
set_dof(model)
Expand Down Expand Up @@ -95,7 +102,7 @@ def test_solve_flowsheet(model):
0.5, 1e-4
)
assert model.fs.M01.outlet.pressure[0].value == pytest.approx(101325, 1e-4)
assert model.fs.M01.outlet.temperature[0].value == pytest.approx(368.2, 1e-4)
assert model.fs.M01.outlet.temperature[0].value == pytest.approx(368.12, 1e-4)

assert model.fs.H02.outlet.flow_mol[0].value == pytest.approx(2.0, 1e-4)
assert model.fs.H02.outlet.mole_frac_comp[0, "benzene"].value == pytest.approx(
Expand Down

0 comments on commit afde76a

Please sign in to comment.