Examples
Examples can be found in pyomo/contrib/parmest/examples and include:
Reactor design example [PyomoBookII]
Semibatch example [SemiBatch]
Rooney Biegler example [RooneyBiegler]
Each example includes a Python file that contains the Pyomo model and a Python file to run parameter estimation.
Additional use cases include:
Data reconciliation (reactor design example)
Parameter estimation using data with duplicate sensors and time-series data (reactor design example)
Parameter estimation using mpi4py, the example saves results to a file for later analysis/graphics (semibatch example)
The example below uses the reactor design example. The file reactor_design.py includes a function which returns an populated instance of the Pyomo model. Note that the model is defined to maximize cb and that k1, k2, and k3 are fixed. The _main_ program is included for easy testing of the model declaration.
# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2024
# National Technology and Engineering Solutions of Sandia, LLC
# Under the terms of Contract DE-NA0003525 with National Technology and
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
# rights in this software.
# This software is distributed under the 3-clause BSD License.
# ___________________________________________________________________________
"""
Continuously stirred tank reactor model, based on
pyomo/examples/doc/pyomobook/nonlinear-ch/react_design/ReactorDesign.py
"""
from pyomo.common.dependencies import pandas as pd
import pyomo.environ as pyo
import pyomo.contrib.parmest.parmest as parmest
from pyomo.contrib.parmest.experiment import Experiment
def reactor_design_model():
# Create the concrete model
model = pyo.ConcreteModel()
# Rate constants
model.k1 = pyo.Param(
initialize=5.0 / 6.0, within=pyo.PositiveReals, mutable=True
) # min^-1
model.k2 = pyo.Param(
initialize=5.0 / 3.0, within=pyo.PositiveReals, mutable=True
) # min^-1
model.k3 = pyo.Param(
initialize=1.0 / 6000.0, within=pyo.PositiveReals, mutable=True
) # m^3/(gmol min)
# Inlet concentration of A, gmol/m^3
model.caf = pyo.Param(initialize=10000, within=pyo.PositiveReals, mutable=True)
# Space velocity (flowrate/volume)
model.sv = pyo.Param(initialize=1.0, within=pyo.PositiveReals, mutable=True)
# Outlet concentration of each component
model.ca = pyo.Var(initialize=5000.0, within=pyo.PositiveReals)
model.cb = pyo.Var(initialize=2000.0, within=pyo.PositiveReals)
model.cc = pyo.Var(initialize=2000.0, within=pyo.PositiveReals)
model.cd = pyo.Var(initialize=1000.0, within=pyo.PositiveReals)
# Objective
model.obj = pyo.Objective(expr=model.cb, sense=pyo.maximize)
# Constraints
model.ca_bal = pyo.Constraint(
expr=(
0
== model.sv * model.caf
- model.sv * model.ca
- model.k1 * model.ca
- 2.0 * model.k3 * model.ca**2.0
)
)
model.cb_bal = pyo.Constraint(
expr=(0 == -model.sv * model.cb + model.k1 * model.ca - model.k2 * model.cb)
)
model.cc_bal = pyo.Constraint(
expr=(0 == -model.sv * model.cc + model.k2 * model.cb)
)
model.cd_bal = pyo.Constraint(
expr=(0 == -model.sv * model.cd + model.k3 * model.ca**2.0)
)
return model
class ReactorDesignExperiment(Experiment):
def __init__(self, data, experiment_number):
self.data = data
self.experiment_number = experiment_number
self.data_i = data.loc[experiment_number, :]
self.model = None
def create_model(self):
self.model = m = reactor_design_model()
return m
def finalize_model(self):
m = self.model
# Experiment inputs values
m.sv = self.data_i['sv']
m.caf = self.data_i['caf']
# Experiment output values
m.ca = self.data_i['ca']
m.cb = self.data_i['cb']
m.cc = self.data_i['cc']
m.cd = self.data_i['cd']
return m
def label_model(self):
m = self.model
m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
m.experiment_outputs.update(
[
(m.ca, self.data_i['ca']),
(m.cb, self.data_i['cb']),
(m.cc, self.data_i['cc']),
(m.cd, self.data_i['cd']),
]
)
m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
m.unknown_parameters.update(
(k, pyo.ComponentUID(k)) for k in [m.k1, m.k2, m.k3]
)
return m
def get_labeled_model(self):
m = self.create_model()
m = self.finalize_model()
m = self.label_model()
return m
def main():
# For a range of sv values, return ca, cb, cc, and cd
results = []
sv_values = [1.0 + v * 0.05 for v in range(1, 20)]
caf = 10000
for sv in sv_values:
# make model
model = reactor_design_model()
# add caf, sv
model.caf = caf
model.sv = sv
# solve model
solver = pyo.SolverFactory("ipopt")
solver.solve(model)
# save results
results.append([sv, caf, model.ca(), model.cb(), model.cc(), model.cd()])
results = pd.DataFrame(results, columns=["sv", "caf", "ca", "cb", "cc", "cd"])
print(results)
if __name__ == "__main__":
main()
The file parameter_estimation_example.py uses parmest to estimate values of k1, k2, and k3 by minimizing the sum of squared error between model and observed values of ca, cb, cc, and cd. Additional example files use parmest to run parameter estimation with bootstrap resampling and perform a likelihood ratio test over a range of theta values.
# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2024
# National Technology and Engineering Solutions of Sandia, LLC
# Under the terms of Contract DE-NA0003525 with National Technology and
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
# rights in this software.
# This software is distributed under the 3-clause BSD License.
# ___________________________________________________________________________
from pyomo.common.dependencies import pandas as pd
from os.path import join, abspath, dirname
import pyomo.contrib.parmest.parmest as parmest
from pyomo.contrib.parmest.examples.reactor_design.reactor_design import (
ReactorDesignExperiment,
)
def main():
# Read in data
file_dirname = dirname(abspath(str(__file__)))
file_name = abspath(join(file_dirname, "reactor_data.csv"))
data = pd.read_csv(file_name)
# Create an experiment list
exp_list = []
for i in range(data.shape[0]):
exp_list.append(ReactorDesignExperiment(data, i))
# View one model
# exp0_model = exp_list[0].get_labeled_model()
# exp0_model.pprint()
pest = parmest.Estimator(exp_list, obj_function='SSE')
# Parameter estimation with covariance
obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=17)
print(obj)
print(theta)
The semibatch and Rooney Biegler examples are defined in a similar manner.