Hide code cell content
###############################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2023 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved.  Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
###############################################################################

Autothermal Reformer Flowsheet Optimization with PySMO Surrogate Object#

Maintainer: Brandon Paul
Author: Brandon Paul
Updated: 2023-06-01

1. Introduction#

This example demonstrates autothermal reformer optimization leveraging the PySMO Polynomial surrogate trainer. Other than the specific training method syntax, this workflow is identical for PySMO RBF and PySMO Kriging surrogate models. In this notebook, sampled simulation data will be used to train and validate a surrogate model. IDAES surrogate plotting tools will be utilized to visualize the surrogates on training and validation data. Once validated, integration of the surrogate into an IDAES flowsheet will be demonstrated.

2. Problem Statement#

Within the context of a larger NGFC system, the autothermal reformer generates syngas from air, steam and natural gas for use in a solid-oxide fuel cell (SOFC).

2.1. Main Inputs:#

  • Bypass fraction (dimensionless) - split fraction of natural gas to bypass AR unit and feed directly to the power island

  • NG-Steam Ratio (dimensionless) - proportion of natural relative to steam fed into AR unit operation

2.2. Main Outputs:#

  • Steam flowrate (kg/s) - inlet steam fed to AR unit

  • Reformer duty (kW) - required energy input to AR unit

  • Composition (dimensionless) - outlet mole fractions of components (Ar, C2H6, C3H8, C4H10, CH4, CO, CO2, H2, H2O, N2, O2)

from IPython.display import Image
from pathlib import Path


def datafile_path(name):
    return Path("..") / name


Image(datafile_path("AR_PFD.png"))
../../../_images/1a3b4ae911ad91d018a0eda19303684d24f174506309b4831ccbdaefc6dc984b.png

3. Training and Validating Surrogates#

First, let’s import the required Python, Pyomo and IDAES modules:

# Import statements
import os
import numpy as np
import pandas as pd

# Import Pyomo libraries
from pyomo.environ import (
    ConcreteModel,
    SolverFactory,
    value,
    Var,
    Constraint,
    Set,
    Objective,
    maximize,
)
from pyomo.common.timing import TicTocTimer

# Import IDAES libraries
from idaes.core.surrogate.sampling.data_utils import split_training_validation
from idaes.core.surrogate.pysmo_surrogate import PysmoPolyTrainer, PysmoSurrogate
from idaes.core.surrogate.plotting.sm_plotter import (
    surrogate_scatter2D,
    surrogate_parity,
    surrogate_residual,
)
from idaes.core.surrogate.surrogate_block import SurrogateBlock
from idaes.core import FlowsheetBlock

3.1 Importing Training and Validation Datasets#

In this section, we read the dataset from the CSV file located in this directory. 2800 data points were simulated from a rigorous IDAES NGFC flowsheet using a grid sampling method. For simplicity and to reduce training runtime, this example randomly selects 100 data points to use for training/validation. The data is separated using an 80/20 split into training and validation data using the IDAES split_training_validation() method.

# Import Auto-reformer training data
np.set_printoptions(precision=6, suppress=True)

csv_data = pd.read_csv(datafile_path("reformer-data.csv"))  # 2800 data points
data = csv_data.sample(n=100)  # randomly sample points for training/validation
input_data = data.iloc[:, :2]
output_data = data.iloc[:, 2:]

# Define labels, and split training and validation data
# note that PySMO requires that labels are passed as string lists
input_labels = list(input_data.columns)
output_labels = list(output_data.columns)

n_data = data[input_labels[0]].size
data_training, data_validation = split_training_validation(
    data, 0.8, seed=n_data
)  # seed=100

3.2 Training Surrogates with PySMO#

IDAES builds a model class for each type of PySMO surrogate model. In this case, we will call and build the Polynomial Regression class. Regression settings can be directly passed as class arguments, as shown below. In this example, allowed basis terms span a 6th order polynomial as well as a variable product, and data is internally cross-validated using 10 iterations of 80/20 splits to ensure a robust surrogate fit. Note that PySMO uses cross-validation of training data to adjust model coefficients and ensure a more accurate fit, while we separate the validation dataset pre-training in order to visualize the surrogate fits.

Finally, after training the model we save the results and model expressions to a folder which contains a serialized JSON file. Serializing the model in this fashion enables importing a previously trained set of surrogate models into external flowsheets. This feature will be used later.

# capture long output (not required to use surrogate API)
from io import StringIO
import sys

stream = StringIO()
oldstdout = sys.stdout
sys.stdout = stream

# Create PySMO trainer object
trainer = PysmoPolyTrainer(
    input_labels=input_labels,
    output_labels=output_labels,
    training_dataframe=data_training,
)

# Set PySMO options
trainer.config.maximum_polynomial_order = 6
trainer.config.multinomials = True
trainer.config.training_split = 0.8
trainer.config.number_of_crossvalidations = 10

# Train surrogate (calls PySMO through IDAES Python wrapper)
poly_train = trainer.train_surrogate()

# create callable surrogate object
xmin, xmax = [0.1, 0.8], [0.8, 1.2]
input_bounds = {input_labels[i]: (xmin[i], xmax[i]) for i in range(len(input_labels))}
poly_surr = PysmoSurrogate(poly_train, input_labels, output_labels, input_bounds)

# save model to JSON
model = poly_surr.save_to_file("pysmo_poly_surrogate.json", overwrite=True)

# revert back to normal output capture
sys.stdout = oldstdout

# display first 50 lines and last 50 lines of output
celloutput = stream.getvalue().split("\n")
for line in celloutput[:50]:
    print(line)
print(".")
print(".")
print(".")
for line in celloutput[-50:]:
    print(line)
---------------------------------------------------------------------------
ApplicationError                          Traceback (most recent call last)
Cell In[5], line 23
     20 trainer.config.number_of_crossvalidations = 10
     22 # Train surrogate (calls PySMO through IDAES Python wrapper)
---> 23 poly_train = trainer.train_surrogate()
     25 # create callable surrogate object
     26 xmin, xmax = [0.1, 0.8], [0.8, 1.2]

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/surrogate/pysmo_surrogate.py:196, in PysmoTrainer.train_surrogate(self)
    194 if hasattr(self, "_input_bounds"):
    195     self._trained.input_bounds = self._input_bounds
--> 196 self._training_main_loop()
    197 return self._trained

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/surrogate/pysmo_surrogate.py:223, in PysmoTrainer._training_main_loop(self)
    221 # Create and train model
    222 model = self._create_model(pysmo_input, output_label)
--> 223 model.training()
    224 # Store results
    225 result = PysmoSurrogateTrainingResult()

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/surrogate/pysmo/polynomial_regression.py:1642, in PolynomialRegression.training(self)
   1638 npe = NumpyEvaluator(cMap)
   1639 additional_data = list(
   1640     npe.walk_expression(term) for term in self.additional_term_expressions
   1641 )
-> 1642 return self.polynomial_regression_fitting(additional_data)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/surrogate/pysmo/polynomial_regression.py:1207, in PolynomialRegression.polynomial_regression_fitting(self, additional_regression_features)
   1205 for poly_order in range(1, self.max_polynomial_order + 1):
   1206     for cv_number in range(1, self.number_of_crossvalidations + 1):
-> 1207         phi, train_error, cv_error = self.polyregression(
   1208             poly_order,
   1209             training_data["training_set_" + str(cv_number)],
   1210             cross_val_data["test_set_" + str(cv_number)],
   1211         )
   1212         if cv_error < best_error:
   1213             best_error = cv_error

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/surrogate/pysmo/polynomial_regression.py:900, in PolynomialRegression.polyregression(self, poly_order, training_data, test_data, additional_x_training_data, additional_x_test_data)
    896     phi_vector = self.bfgs_parameter_optimization(
    897         x_polynomial_data, y_training_data
    898     )
    899 elif self.solution_method == "pyomo":
--> 900     phi_vector = self.pyomo_optimization(x_polynomial_data, y_training_data)
    901 phi_vector = phi_vector.reshape(
    902     phi_vector.shape[0], 1
    903 )  # Pseudo-inverse approach
    905 x_polynomial_data_test = self.polygeneration(
    906     poly_order, self.multinomials, x_test_data, additional_x_test_data
    907 )

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/surrogate/pysmo/polynomial_regression.py:816, in PolynomialRegression.pyomo_optimization(x, y)
    814 opt.options["max_iter"] = 1000
    815 # TODO: Should this be checking the for a feasible solution?
--> 816 opt.solve(instance)
    818 # Convert theta variable into numpy array
    819 phi = np.zeros((len(instance.theta), 1))

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/base/solvers.py:534, in OptSolver.solve(self, *args, **kwds)
    531 def solve(self, *args, **kwds):
    532     """Solve the problem"""
--> 534     self.available(exception_flag=True)
    535     #
    536     # If the inputs are models, then validate that they have been
    537     # constructed! Collect suffix names to try and import from solution.
    538     #
    539     from pyomo.core.base.block import _BlockData

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:139, in SystemCallSolver.available(self, exception_flag)
    137     if exception_flag:
    138         msg = "No executable found for solver '%s'"
--> 139         raise ApplicationError(msg % self.name)
    140     return False
    141 return True

ApplicationError: No executable found for solver 'ipopt'

3.3 Visualizing surrogates#

Now that the surrogate models have been trained, the models can be visualized through scatter, parity and residual plots to confirm their validity in the chosen domain. The training data will be visualized first to confirm the surrogates are fit the data, and then the validation data will be visualized to confirm the surrogates accurately predict new output values.

# visualize with IDAES surrogate plotting tools
surrogate_scatter2D(poly_surr, data_training, filename="pysmo_poly_train_scatter2D.pdf")
surrogate_parity(poly_surr, data_training, filename="pysmo_poly_train_parity.pdf")
surrogate_residual(poly_surr, data_training, filename="pysmo_poly_train_residual.pdf")
../../../_images/24fdad58226d2bae1058ce796d8b0ea16072fe3eaff80c6d7be997b8283fcea6.png ../../../_images/5d19029cc5d0f120785b6c757204a9271d8b57f566f481159255d0a37d95b04a.png ../../../_images/c2d25db5f66dd5be2624a8b916848c54f47edc0ad0d7fbb8ab17143d3d73e638.png ../../../_images/8fbf7c98dee1b7c3d823706331c819a801600e6dbe2a92409c72de25e2168067.png ../../../_images/eb2e98fa285c6fbd2074b0bc8343e29e60c4945eabe52dc3cecd5cc769752907.png ../../../_images/8be5a6e1253e4a4f144c2c331e4259f4857bff8ae48cec6e5d77591dfb0c9172.png ../../../_images/5ca527899a934cd2e022d45de1b886e62a545e49d75935bb842b6e5c82dfd490.png ../../../_images/41454b246058d75bc7b3e81ad3b62c91e89b6e9f5d1558971d22a57cf9d7543d.png ../../../_images/a61f8ba1306dfc4507dec73ceb556f52c59da01b3040cc61529c0347d1c75cc8.png ../../../_images/8335de2784051ebe82206d3810826a11723fee8193e48bcdc393fa6d45f5057e.png ../../../_images/4c78fae8dd320e1d66b73a2a9c94b4a190b36fdec285716ec7993aa33b1fdf18.png ../../../_images/dae85dca525ad4ecff34aae4578bc2866d9e6e8a388b03f3dc2baa7344749c2e.png ../../../_images/e7d8d0838860a2483b00a3828ef2f25cbe5a9852a0717d7509455f22dcacb258.png ../../../_images/24a49f7793d4964e0721eb3e5d123c506bdb0258463a1ffc4e4c1d3097a0aa0e.png ../../../_images/402b3d33f254f63e1c75706d1628fabde01d363bf9c5f6f6157878ed11e5a366.png ../../../_images/84a6d4e6420c204619984c4217027300d2b3f25863141eb9a0bffe9c8014dbfa.png ../../../_images/11ddb1799b2b35d007932523ebedd9cc3d3b2a33ec0617d6749f3824630bacde.png ../../../_images/4a899cc23f5bcb6c2a53a9e51559f3ffcd1912d44aff79feb76b836b74d6551c.png ../../../_images/91be969bf16b2f094084ce5014638dec2da6e6c79c193da17b164cf1b8fa9c8c.png ../../../_images/4d6ed4d2adb1bfb472816de65abaf2b8ce9ad221ce2ff7d991eab1102775e086.png ../../../_images/fdce1b581f4661246b416c2d50252aef36719fa4f6d99c702ceb2406f0ed4467.png ../../../_images/9a7f43d0bc70e06c5bb6484343d2fc935fefc2ef884453755b586f9c2de7adcb.png ../../../_images/cc83432f37249d8f5cb3d62cfc9ebdf4dda86f7a8a34accc36ff59050d8bbc40.png ../../../_images/fb5f6ea8ca2560036c9607f3c7e5655f97c4bd803af2b012ec504807c7a22308.png ../../../_images/85e4a39691c3adcc026d8be0e7e94637c8e410cafc0225b161f0f06cef2055c2.png ../../../_images/f8cacaf8d253aeecd842b6440b8b58520094505651d42eda2ed8a404aa6c50d9.png ../../../_images/2285cc41007489bf1ac80cfac19607d3de6cbc8bbfd1b13bca07b8e1e0c152b0.png ../../../_images/0ed329935b6c31e5b6e0d9f3d92157143c4420d5bbeb2afec0144670a8c4ba03.png ../../../_images/1dd4747d1d49e9486f4c391496f12f92f05084aa915d35690d5d93e5158396cb.png ../../../_images/331c29197094e521019d0be7e150f675bff92cec010e71d8c6986e4794c97e42.png ../../../_images/aacdadf79287e5c4435c57b85a077d03b9f4bd99f121b6a5afc2ee161214d307.png ../../../_images/7a071e6f353f29d3a9bac5fd1cc67dfd85157d765921a1cfe4f8455cf4424aff.png ../../../_images/09be0b975d8514c3c9fde3827ba05a0acbeb86479796563c06cdde00e0ac434a.png ../../../_images/44d01efa5dce1c46fd8226055782deaf2a39bcb5007c14ad00ed24e1c9f50c65.png ../../../_images/d0a6d5a63c4874091e360f68d5acee289456ae6d8899e1a03c8dc735ee7881a3.png ../../../_images/3b8d6753a2f27bbdf11da16b0985a3a675dd2c7dee62f723710bfa747ac79b4f.png ../../../_images/76870e17558d36caaf5718c4ab55a01c121cc4892cfdd2c4d7794a9d894af52c.png ../../../_images/fac5d81a797d0573fd93b7f9792bbb16f9b1d0dcd82e5f2813803a6f9fb34863.png ../../../_images/d7a74fe95960c178a276826dfad7af737c086c8d69f755ac29e9c5aed2f7dd85.png ../../../_images/86dcf7887025725010750dcaa298dc86267a7b09ec3b832fbb8ebc6c733aa2b3.png ../../../_images/ff0defb7a4dfda6df3724c04412e3e23d4374dd901770b6886750cb8cd8af50f.png ../../../_images/194e9ece2c761b6906d7d73303f8b75e1122100d45bc4f69d49488a1dceb09d9.png ../../../_images/373c056668c34705f7c0ece3d9a4482e9d7de54877cbf8a72449d9ecb2ac5b1a.png ../../../_images/a969b6818181b39fa6043bebf646e794b0bdd9ba1d3c88755e9186e962adaf7b.png ../../../_images/5c5665cc018fc9b20ad104d12421688bc11427aa0b8264971e7b22c709f32581.png ../../../_images/c7022ef33eab89d507d18616247e0107af525c7152361b04302dd6157ac4ee94.png ../../../_images/bfedb810f1d7cd40682cd6b20ad9e5b8113f98bd86f3091d0335c27a4187fd2e.png ../../../_images/f0c30416b9da47f91cfa0b098729bf5c347492dcac096fce7a9df23a42fcaa97.png ../../../_images/2a0ab975b2c90eeb000f5cca9357e55b0bbed612f02dd00098afff40ebf34011.png ../../../_images/a78413b2a637bdf2d0448de69e5109ad0dc441302a5004ed42dfda21d643ffc8.png ../../../_images/76b3ee537ec749c7de82a54173631040c6342e73858ccd44788a396d5b017f5f.png ../../../_images/7a8f7ac641fca5205827c6f3aa429c8690a237378a44206abdd97cad65f63c40.png ../../../_images/95f657d2168cad06aed62748349bf0cca5b0792fefb323ed5a85df70dd25f61d.png ../../../_images/1bbb60b15240dd3b65c491edfbca57ceb70f3c834a619f7d13605ad817434467.png ../../../_images/163de550596e474892b350746be96387a548cb2c8542511799c1ea54a536d598.png ../../../_images/70f60f0ed4008d2f600d4ff9990a1595f32a7d75320899a5bc9179e20cc72789.png ../../../_images/91a8f186a3b07da05fd010c17b4807a44fcc75f31cfac11cedce57aa43f58eaa.png ../../../_images/90cb17bbdf7f9516d174cc79d48c2f0aefd9ae7d13c8e0dae701c742a5ced464.png ../../../_images/ba29f1f386f2db4f66803eb43003667e4e7323dc8c9a493b0f327a6afadbe2d1.png ../../../_images/c08e65151d27b91f7b712677ec26dded0b823e3b1131ca3840fd3b462c940c8b.png ../../../_images/7b175f04de07dd093a78d491c1b8c1dc84545946198c5a7156d7d0381d6f2813.png ../../../_images/80797f61eba2ead666b4953037320fbd0b1b83197fd4941a4f63e71a056638d4.png ../../../_images/e0854a1fb7a1253a9c41e3981920455964e8c434067eaffff7d30d28932feeae.png ../../../_images/c0628cc3b2d6188f565db24017660747da78fdf907b43b2c66a9d9c9aad02213.png ../../../_images/008a4a7f7aeb4fe23c98c18d9bbe90792b536f813769a30b27cc1f5aaf23c542.png

3.4 Model Validation#

# visualize with IDAES surrogate plotting tools
surrogate_scatter2D(poly_surr, data_validation, filename="pysmo_poly_val_scatter2D.pdf")
surrogate_parity(poly_surr, data_validation, filename="pysmo_poly_val_parity.pdf")
surrogate_residual(poly_surr, data_validation, filename="pysmo_poly_val_residual.pdf")
../../../_images/2db40c0740d4a553508e64dd663eb7c4be9d6de72bfe1c6bd58d4bea40af91e9.png ../../../_images/13fa956edc6519671249cd736f0c2927952e974815cee132f12df27cca9c050c.png ../../../_images/c6504f7c19e7502b4bdf80608feb5f9e5494b82169ff4b2754331d569019cea2.png ../../../_images/bbad5d2c7a057b12737bf1213c5eb4004ba4ce80c9d070310e19d905eab97a37.png ../../../_images/f1cd66e8055ecfd6448449a244297a08ad6ecdd103f56531f3fcba05518f455d.png ../../../_images/dba402eebb11cacc1ec74ddfe6ae1c51cdd1a99e1b6bd510cbc77ad4e3497ad3.png ../../../_images/7863b59cab4d9d170aa65333c7eee4c704e9d06ffd0bd1081119f2db3733122f.png ../../../_images/946389babc29cb345d77c28ebeee4d5b21aa4b149f80bf7442e7c62b578ee43d.png ../../../_images/f917bd022bd636275a5c822f38dac91c33e41109795abebc11729ce715e04a1c.png ../../../_images/d7d2acf396060df683f1a8f33b0393bd03da845dc1be1165d324c9fb45209c4b.png ../../../_images/7aab11909dfd1bb8b90057c2ecfd5972ae22974517df77a19dcff1e69a3eeaa5.png ../../../_images/d218575bdc0eb05409f19d15ed864d1415661472f6c79e6f3dd924361c4468f4.png ../../../_images/413e660a7eb2c102578c04524ce0168e71482a267703148f3aeef61ae682442b.png ../../../_images/b9d5f938ecef5e2bd2ea624bfde74d063eb214680a11fdf1eb988a589f0ed899.png ../../../_images/4866bbc080d0dbaa81fb1fd701ab4b91ee4950d8f3718974ad179679873baa24.png ../../../_images/a70a8e41555d37e4207f7997b26e443365aeeab19708a3f7ca87088fa8961a42.png ../../../_images/d314db1ec77f580a07c29f547d453c50566c81507b3b98f68740865a1178e146.png ../../../_images/6bf26d54aa52c9d86530da6317f1774e2ff21485365d7f42bb60214bbbf938dc.png ../../../_images/a6526d6a9abda7d01382f7f25a6b12336c4445c0a7de0906d4d0a58cd96de23c.png ../../../_images/ba17ce886d8b990e9f5a8f0442c184e4b60b826dd80ea0da2bc3d4fdabe1e738.png ../../../_images/44821c0debe8b8455050b85a315daea2fbfb6cef70985bb05b6e79ca5b936ec5.png ../../../_images/0cf2d6213a31f952f9d3e3dad40e816522d23327e2690cf869edda598404f527.png ../../../_images/144231c7a5d6d493f731c088acc1e18d6770b33a43f19bd232acbe4cc3924b47.png ../../../_images/45b60640352d358adb43ce35d2bf760e9b35677ea981cfb486adc715ab811194.png ../../../_images/3013cfb4ffae1753f46f88a0cddea329fbaea82305b82d20ac82287e040309f7.png ../../../_images/feeb79921283fa0be834f041cdd24ea8a01f3db69260a96b6846bb13cc448143.png ../../../_images/ae3b9b024c138adb8372fdd7f86d4b801183c618b13f64154e81df728513540e.png ../../../_images/fbcca0d238441d9984b0300202d329f04d395fd125e5be448e18a4fc615e19be.png ../../../_images/4b9f9dbffd10c11b0f92000ba32c67573adba2f4c7784ee39c30cafc076e9612.png ../../../_images/7c18d9388f464b1b5b9d038431916087307c13a93cfd0b70c7b3ed3cf877349d.png ../../../_images/f512d0ca370a962a628287b707fd203bf3295e9710f6d0bb2aa52ffdf760eae2.png ../../../_images/04f79effd05c29f7f71ff1e646c6c21927ebabb62de7460f1794fb33bae3df32.png ../../../_images/384959b22ccbb6cd138b1f53f7ef3c146616e00a6d7c7967e225839306034e88.png ../../../_images/f994bd7a3b4e37e978cb5e4fcfac76fc2ac804f55a5c5b12004eaab20edba626.png ../../../_images/3adda051b35f6193472804b1c6754d82bf996c8f6793d4c377807b327e4262b9.png ../../../_images/3b56197cbd795d8c8a4a9f1973d17464bd4c15ac362f1484fcd32ecb3f0f13dc.png ../../../_images/385665a668b4dfa9ffc6aa25bb6f96a482c702a685df711fbadf88c0b4a09235.png ../../../_images/5fa637d704ebde910f8696756b13e1cc9e7a47dcaec86f0a939ba82f72b20d8f.png ../../../_images/63d6c9163eafc005cec28ac41601c3bdb41064bc0b4674d432c2e5e1566bd226.png ../../../_images/f5515e1b1aaa9408610326cf55b91d76e77e76fbc3487fd2a81da13e28829801.png ../../../_images/1ee207bec84a21c0cad82298b40a7c6b911b5bf2ec847b5eed6075fc48703712.png ../../../_images/9920d1fc6187094fc7164acf868a111e2790d46e36224e588ae2d9e6e6bb7df2.png ../../../_images/e930fe69bf3d7c08fb20a77b2d5d96c3604d058f5a96a4f29b90c020ea463350.png ../../../_images/588d8b01eab5d79c9b281f2ea15522d05d53293577e8d6a692b905aadcf8c01c.png ../../../_images/daa09467e0a3a8161aa9bbe3184e40a1227e7dc9a0b984d47cda6da807a19485.png ../../../_images/609d1bae57c6743a3809aadee5af4a1e74e4a26c9546e1d49bde39e6df610812.png ../../../_images/3cff81a48a33bd451eba382e1d78d118ad49406e660ffc055c88995e8ae925cb.png ../../../_images/cb64d9d6b5f4c479b30dac26e8d0c62c58175da578e03262cc62e70cffbe398a.png ../../../_images/1229b633c183669faae839c0f8b752f32cfcde3b0b1cc6dbf03d83356fad7a40.png ../../../_images/2bf148284cef6f11fdd9a8334535291c1592c9c44d8344e986ef001cd0f51ceb.png ../../../_images/6bc53f9fbb5ab1ee3c14f385d4c10cb5b2cbac708e9add68508c8d5c7fbb42b4.png ../../../_images/16cdc7828613d45a99ef4fab2944d4eabefefe5cd76b2a88028baa45e7435036.png ../../../_images/5a907d9d0d51810c35eee3ce4e58e6ce6e69cc15bb95b0010a20b13687d23c80.png ../../../_images/54ebc93a021bbfda156de3322f3de9ff3e81d1725a585afe005ee06f79cc80b8.png ../../../_images/7ddb3244d036bda34857af98810522550ffb9d150910ea37478bb8c5be606db3.png ../../../_images/215df276d11e21896e6217b69153416cd6bb4ae75d7d7b499e8e3a0d0c34f9a0.png ../../../_images/8b7b4c4c1e17a18825b29dbe94ef73315ca4a9506eee81d5734ab05b742d507e.png ../../../_images/31698956865547a07e477f517741b5f0a400e9dc4fa5a9e9bdbe10af93251cf2.png ../../../_images/3a3dafe78c12229ecf2552a42e31c5795c01236ef26c25f67e9f6a5dc3744fde.png ../../../_images/8e955a38a89c3c2067cbb44a5f578e971eff0dad6255d0d35d7fe75121e4b015.png ../../../_images/166cac8e3ef514c9e91c4bce2d603fc1cc10633024b7c030ec19a342a5fca48a.png ../../../_images/f80eee74ad2827b1d8e26f86fc9a137d1d7ce029545ad64ba2a5e93ebad43e73.png ../../../_images/4fb9bc2c4e71d4c99dbe4ed17510af8e7b9fa018f93ec17b3a67d8523b887a48.png ../../../_images/5c40c186942643d08791197ebc6e8b71010988e70a7023e68ac4ae26010f3e1e.png ../../../_images/02957c75db60e767be9647bdb3e1422e4851ef35412e687ac02ba5e129ecb9b8.png

4. IDAES Flowsheet Integration#

4.1 Build and Run IDAES Flowsheet#

Next, we will build an IDAES flowsheet and import the surrogate model object. Each output variable has a unique PySMO model expression, and the surrogate expressions may be added to the model via an indexed Constraint() component.

# create the IDAES model and flowsheet
m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)

# create flowsheet input variables
m.fs.bypass_frac = Var(
    initialize=0.80, bounds=[0.1, 0.8], doc="natural gas bypass fraction"
)
m.fs.ng_steam_ratio = Var(
    initialize=0.80, bounds=[0.8, 1.2], doc="natural gas to steam ratio"
)

# create flowsheet output variables
m.fs.steam_flowrate = Var(initialize=0.2, doc="steam flowrate")
m.fs.reformer_duty = Var(initialize=10000, doc="reformer heat duty")
m.fs.AR = Var(initialize=0, doc="AR fraction")
m.fs.C2H6 = Var(initialize=0, doc="C2H6 fraction")
m.fs.C3H8 = Var(initialize=0, doc="C3H8 fraction")
m.fs.C4H10 = Var(initialize=0, doc="C4H10 fraction")
m.fs.CH4 = Var(initialize=0, doc="CH4 fraction")
m.fs.CO = Var(initialize=0, doc="CO fraction")
m.fs.CO2 = Var(initialize=0, doc="CO2 fraction")
m.fs.H2 = Var(initialize=0, doc="H2 fraction")
m.fs.H2O = Var(initialize=0, doc="H2O fraction")
m.fs.N2 = Var(initialize=0, doc="N2 fraction")
m.fs.O2 = Var(initialize=0, doc="O2 fraction")

# create input and output variable object lists for flowsheet
inputs = [m.fs.bypass_frac, m.fs.ng_steam_ratio]
outputs = [
    m.fs.steam_flowrate,
    m.fs.reformer_duty,
    m.fs.AR,
    m.fs.C2H6,
    m.fs.C4H10,
    m.fs.C3H8,
    m.fs.CH4,
    m.fs.CO,
    m.fs.CO2,
    m.fs.H2,
    m.fs.H2O,
    m.fs.N2,
    m.fs.O2,
]

# create the Pyomo/IDAES block that corresponds to the surrogate
# PySMO

# capture long output (not required to use surrogate API)
stream = StringIO()
oldstdout = sys.stdout
sys.stdout = stream

surrogate = PysmoSurrogate.load_from_file("pysmo_poly_surrogate.json")
m.fs.surrogate = SurrogateBlock(concrete=True)
m.fs.surrogate.build_model(surrogate, input_vars=inputs, output_vars=outputs)

# revert back to normal output capture - don't need to print PySMO load output
sys.stdout = oldstdout

# fix input values and solve flowsheet
m.fs.bypass_frac.fix(0.5)
m.fs.ng_steam_ratio.fix(1)

solver = SolverFactory("ipopt")
results = solver.solve(m, tee=True)
2023-11-02 10:23:37 [INFO] idaes.core.surrogate.pysmo_surrogate: Decode surrogate. type=poly
Ipopt 3.13.2: 

******************************************************************************
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...:       13
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:       13
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:       13
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 1.11e+04 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 0.00e+00 0.00e+00  -1.0 1.11e+04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.


Let’s print some model results:

print("Steam flowrate = ", value(m.fs.steam_flowrate))
print("Reformer duty = ", value(m.fs.reformer_duty))
print("Mole Fraction Ar = ", value(m.fs.AR))
print("Mole Fraction C2H6 = ", value(m.fs.C2H6))
print("Mole Fraction C3H8 = ", value(m.fs.C3H8))
print("Mole Fraction C4H10 = ", value(m.fs.C4H10))
print("Mole Fraction CH4 = ", value(m.fs.CH4))
print("Mole Fraction CO = ", value(m.fs.CO))
print("Mole Fraction CO2 = ", value(m.fs.CO2))
print("Mole Fraction H2 = ", value(m.fs.H2))
print("Mole Fraction H2O = ", value(m.fs.H2O))
print("Mole Fraction N2 = ", value(m.fs.N2))
print("Mole Fraction O2 = ", value(m.fs.O2))
Steam flowrate =  0.6059308497978191
Reformer duty =  21066.948701158035
Mole Fraction Ar =  0.0036793744229938544
Mole Fraction C2H6 =  0.004187279676589231
Mole Fraction C3H8 =  0.000523414598169937
Mole Fraction C4H10 =  0.0009159732540583096
Mole Fraction CH4 =  0.12786005023329045
Mole Fraction CO =  0.09697157382062967
Mole Fraction CO2 =  0.046010703278916036
Mole Fraction H2 =  0.2938730753304199
Mole Fraction H2O =  0.11952683194799225
Mole Fraction N2 =  0.30644275497350865
Mole Fraction O2 =  -5.551115123125783e-17

4.2 Optimizing the Autothermal Reformer#

Extending this example, we will unfix the input variables and optimize hydrogen production. We will restrict nitrogen below 34 mol% of the product stream and leave all other variables unfixed.

Above, variable values are called in reference to actual objects names; however, as shown below this may be done much more compactly by calling the list objects we created earlier.

# unfix input values and add the objective/constraint to the model
m.fs.bypass_frac.unfix()
m.fs.ng_steam_ratio.unfix()
m.fs.obj = Objective(expr=m.fs.H2, sense=maximize)
m.fs.con = Constraint(expr=m.fs.N2 <= 0.34)

# solve the model
tmr = TicTocTimer()
status = solver.solve(m, tee=True)
solve_time = tmr.toc("solve")

# print and check results
assert abs(value(m.fs.H2) - 0.33) <= 0.01
assert value(m.fs.N2 <= 0.4 + 1e-8)
print("Model status: ", status)
print("Solve time: ", solve_time)
for var in inputs:
    print(var.name, ": ", value(var))
for var in outputs:
    print(var.name, ": ", value(var))
Ipopt 3.13.2: 

******************************************************************************
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...:       39
Number of nonzeros in inequality constraint Jacobian.:        1
Number of nonzeros in Lagrangian Hessian.............:        3

Total number of variables............................:       15
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        2
                     variables with only upper bounds:        0
Total number of equality constraints.................:       13
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -2.9387308e-01 2.91e-10 2.31e-02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -2.9586866e-01 1.29e+00 1.91e-03  -1.7 4.88e+02    -  1.00e+00 1.00e+00f  1
   2 -3.1973235e-01 1.17e+02 9.90e-03  -2.5 6.00e+03    -  8.80e-01 1.00e+00h  1
   3 -3.2537868e-01 1.23e+02 3.60e-03  -2.5 4.80e+03    -  1.00e+00 1.00e+00h  1
   4 -3.2610698e-01 2.48e+01 3.26e-04  -2.5 1.90e+03    -  1.00e+00 1.00e+00h  1
   5 -3.2615002e-01 1.98e-02 2.27e-06  -2.5 9.78e+01    -  1.00e+00 1.00e+00h  1
   6 -3.3107576e-01 8.06e+01 5.01e-03  -3.8 3.89e+03    -  9.20e-01 1.00e+00h  1
   7 -3.3136311e-01 1.95e+00 8.59e-04  -3.8 6.41e+02    -  1.00e+00 1.00e+00h  1
   8 -3.3128690e-01 2.96e-01 1.24e-05  -3.8 2.18e+02    -  1.00e+00 1.00e+00h  1
   9 -3.3157177e-01 2.17e-01 1.01e-04  -5.7 1.63e+02    -  9.95e-01 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -3.3157586e-01 4.17e-04 1.45e-08  -5.7 9.24e+00    -  1.00e+00 1.00e+00h  1
  11 -3.3157954e-01 4.46e-05 7.07e-09  -8.6 1.84e+00    -  1.00e+00 1.00e+00h  1
  12 -3.3157954e-01 1.34e-09 4.26e-14  -8.6 1.29e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:  -3.3157953843145921e-01   -3.3157953843145921e-01
Dual infeasibility......:   4.2581604889876294e-14    4.2581604889876294e-14
Constraint violation....:   3.5549293562270731e-12    1.3387762010097504e-09
Complementarity.........:   2.5059038878794521e-09    2.5059038878794521e-09
Overall NLP error.......:   2.5059038878794521e-09    2.5059038878794521e-09


Number of objective function evaluations             = 13
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 13
Number of inequality constraint evaluations          = 13
Number of equality constraint Jacobian evaluations   = 13
Number of inequality constraint Jacobian evaluations = 13
Number of Lagrangian Hessian evaluations             = 12
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.001

EXIT: Optimal Solution Found.

[+   0.07] solve
Model status:  
Problem: 
- Lower bound: -inf
  Upper bound: inf
  Number of objectives: 1
  Number of constraints: 14
  Number of variables: 15
  Sense: unknown
Solver: 
- Status: ok
  Message: Ipopt 3.13.2\x3a Optimal Solution Found
  Termination condition: optimal
  Id: 0
  Error rc: 0
  Time: 0.058537960052490234
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

Solve time:  0.074246599979233
fs.bypass_frac :  0.10000006743045317
fs.ng_steam_ratio :  1.1111587695778178
fs.steam_flowrate :  1.211913588633006
fs.reformer_duty :  38820.801141863
fs.AR :  0.004107083638191985
fs.C2H6 :  0.0005392232921631147
fs.C4H10 :  0.00011795111658708467
fs.C3H8 :  6.741873287192471e-05
fs.CH4 :  0.016806518961534973
fs.CO :  0.10494647112247085
fs.CO2 :  0.05346576317928386
fs.H2 :  0.3315795384314592
fs.H2O :  0.14839740391436385
fs.N2 :  0.3400000043352795
fs.O2 :  8.2852873934491285e-16