Skip to content

Commit

Permalink
Merge pull request #4 from Robbybp/atr_implicit
Browse files Browse the repository at this point in the history 
To debug Implicit ATR Flowsheet
  • Loading branch information
@sbugosen
sbugosen authored Aug 2, 2023
2 parents ca44a3a + 3685cb1 commit d385fad
Show file tree
Hide file tree
Showing 2 changed files with 213 additions and 0 deletions.
172 changes: 172 additions & 0 deletions svi/auto_thermal_reformer/implicit_atr_fsheet.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,172 @@
# ___________________________________________________________________________
#
# Surrogate vs. Implicit: Experiments comparing nonlinear optimization
# formulations
#
# Copyright (c) 2023. Triad National Security, LLC. All rights reserved.
#
# This program was produced under U.S. Government contract 89233218CNA000001
# for Los Alamos National Laboratory (LANL), which is operated by Triad
# National Security, LLC for the U.S. Department of Energy/National Nuclear
# Security Administration. All rights in the program are reserved by Triad
# National Security, LLC, and the U.S. Department of Energy/National Nuclear
# Security Administration. The Government is granted for itself and others
# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license
# in this material to reproduce, prepare derivative works, distribute copies
# to the public, perform publicly and display publicly, and to permit others
# to do so.
#
# This software is distributed under the 3-clause BSD license.
# ___________________________________________________________________________

import pyomo.environ as pyo
from pyomo.environ import (
Constraint,
Var,
ConcreteModel,
Expression,
Objective,
TransformationFactory,
value,
units as pyunits,
)
from pyomo.common.collections import ComponentSet
from pyomo.network import Arc
from idaes.core import FlowsheetBlock
from idaes.models.properties.modular_properties import GenericParameterBlock
from idaes.core.util.model_statistics import degrees_of_freedom

from idaes.models.unit_models import (
Mixer,
Heater,
HeatExchanger,
PressureChanger,
GibbsReactor,
Separator,
Feed,
Product,
)
from idaes.core.util.initialization import propagate_state
from idaes.models_extra.power_generation.properties.natural_gas_PR import get_prop
from idaes.core.solvers import get_solver
from idaes.models.unit_models.pressure_changer import ThermodynamicAssumption
from idaes.models.unit_models.heat_exchanger import delta_temperature_underwood_callback
from pyomo.contrib.incidence_analysis import IncidenceGraphInterface
from pyomo.contrib.pynumero.interfaces.external_pyomo_model import ExternalPyomoModel
from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock
from pyomo.contrib.pynumero.algorithms.solvers.implicit_functions import (
CyIpoptSolverWrapper
)

from svi.external import add_external_function_libraries_to_environment

def make_implicit(m):
########### CREATE EXTERNAL PYOMO MODEL FOR THE AUTOTHERMAL REFORMER ###########
full_igraph = IncidenceGraphInterface(m)
reformer_igraph = IncidenceGraphInterface(m.fs.reformer, include_inequality=False)

# Pressure is not an "external variable" (it is fixed), so the pressure
# linking equation doesn't need to be included as "residual"
to_exclude = ComponentSet([m.fs.REF_OUT_expanded.pressure_equality[0]])
residual_eqns = [
con for con in m.fs.REF_OUT_expanded.component_data_objects(
pyo.Constraint, active=True
)
if con not in to_exclude
]

input_vars = [
m.fs.reformer.inlet.temperature[0],
m.fs.reformer.inlet.mole_frac_comp[0, "H2"],
m.fs.reformer.inlet.mole_frac_comp[0, "CO"],
m.fs.reformer.inlet.mole_frac_comp[0, "H2O"],
m.fs.reformer.inlet.mole_frac_comp[0, "CO2"],
m.fs.reformer.inlet.mole_frac_comp[0, "CH4"],
m.fs.reformer.inlet.mole_frac_comp[0, "C2H6"],
m.fs.reformer.inlet.mole_frac_comp[0, "C3H8"],
m.fs.reformer.inlet.mole_frac_comp[0, "C4H10"],
m.fs.reformer.inlet.mole_frac_comp[0, "N2"],
m.fs.reformer.inlet.mole_frac_comp[0, "O2"],
m.fs.reformer.inlet.mole_frac_comp[0, "Ar"],
m.fs.reformer.inlet.flow_mol[0],
m.fs.reformer.inlet.pressure[0],

m.fs.reformer_recuperator.hot_side_inlet.temperature[0],
m.fs.reformer_recuperator.hot_side_inlet.flow_mol[0],

# NOTE:
# - heat_duty is not necessary as it does not appear elsewhere
# in the model
# - pressure is not an "input" as the pressure linking equation
# is not residual (because the reactor's outlet pressure is fixed)
]
input_vars.extend(m.fs.reformer_recuperator.hot_side_inlet.mole_frac_comp[0, :])

external_eqns = list(reformer_igraph.constraints)

to_exclude = ComponentSet(input_vars)
# These variables only appear in reformer_igraph due to a bug in how
# coefficients with values of zero were handled. This should be fixed in
# recent Pyomo main, but we leave this explicit filtering in here as we
# often run using the latest release.
to_exclude.add(m.fs.reformer.lagrange_mult[0, "N"])
to_exclude.add(m.fs.reformer.lagrange_mult[0, "Ar"])
external_vars = [var for var in reformer_igraph.variables if var not in to_exclude]

external_var_set = ComponentSet(external_vars)
external_eqn_set = ComponentSet(external_eqns)
residual_eqn_set = ComponentSet(residual_eqns)

# Bounds on variables in the implicit function can lead to
# undefined derivatives. However, removing these bounds may make
# evaluation errors more likely.
for i, var in enumerate(external_vars):
var.setlb(None)
var.setub(None)
var.domain = pyo.Reals

epm = ExternalPyomoModel(
input_vars,
external_vars,
residual_eqns,
external_eqns,
)

########### CONNECT FLOWSHEET TO THE IMPLICIT AUTOTHERMAL REFORMER ###########

m_implicit = ConcreteModel()
m_implicit.egb = ExternalGreyBoxBlock()
m_implicit.egb.set_external_model(epm, inputs=input_vars)

fullspace_cons = [
con
for con in full_igraph.constraints
if con not in residual_eqn_set and con not in external_eqn_set
]

fullspace_vars = [
var for var in full_igraph.variables if var not in external_var_set
]

m_implicit.fullspace_cons = pyo.Reference(fullspace_cons)
m_implicit.fullspace_vars = pyo.Reference(fullspace_vars)
m_implicit.objective = pyo.Reference([m.fs.obj])

return m_implicit


def main():
from svi.auto_thermal_reformer.fullspace_atr_fsheet import make_optimization_model
m = make_optimization_model()
add_external_function_libraries_to_environment(m)
m_implicit = make_implicit(m)
solver = pyo.SolverFactory("cyipopt", options = {"tol": 1e-6, "max_iter": 300})
solver.solve(m_implicit, tee=True)

m.fs.reformer.report()
m.fs.reformer_recuperator.report()
m.fs.product.report()
m.fs.reformer_bypass.split_fraction.display()

if __name__ == "__main__":
main()
41 changes: 41 additions & 0 deletions svi/auto_thermal_reformer/temp_eval_error.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
# These are variables and values that cause a function evaluation error in a
# particular trial run
name_values = [
('fs.reformer.control_volume.properties_out[0.0].temperature', 30.515034422167613),
('fs.reformer.lagrange_mult[0.0,C]', -12363.6731135706),
('fs.reformer.lagrange_mult[0.0,H]', 22345.118341294714),
('fs.reformer.lagrange_mult[0.0,O]', 134176.83881743395),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,CH4]', -5.736412380102407),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,C2H6]', 31.1229565968468),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,C3H8]', -0.3053376828896246),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,C4H10]', 2.011079823130566),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,O2]', -2109.347169893634),
('fs.reformer.control_volume.properties_out[0.0].flow_mol', 3800.0330581516587),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,H2O]', 2.5660306354833556e-16),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,CO]', 0.2677589642847847),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,CO2]', 0.0010799606585477923),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,N2]', 0.26758271932757727),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,Ar]', 0.0034567019034059637),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[H2]', -8.202910038671492e-08),
('fs.reformer.control_volume.properties_out[0.0].flow_mol_phase[Vap]', 3800.0330581516587),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,H2]', -8.202910038671492e-08),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[CO]', 0.2677589642847847),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[H2O]', 2.5660306354833556e-16),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[CO2]', 0.0010799606585477923),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[CH4]', 0.003226314840405902),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[C2H6]', 1.1226724351338853),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[C3H8]', -0.23108282771177824),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[C4H10]', -0.4346941864077285),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[N2]', 0.26758271932757727),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[O2]', -5.761043957622803e-31),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_comp[Ar]', 0.0034567019034059637),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,H2]', -292.7881820210107),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,CO]', 3.1698132265424377),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,H2O]', 1026.7866442938127),
('fs.reformer.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,CO2]', 1658.3151600250535),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,CH4]', 0.003226314840405902),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,C2H6]', 1.1226724351338853),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,C3H8]', -0.23108282771177824),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,C4H10]', -0.4346941864077285),
('fs.reformer.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,O2]', -5.761043957622803e-31),
]

0 comments on commit d385fad

Please sign in to comment.