Examples

Examples can be found in pyomo/contrib/parmest/examples and include:

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 description 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-2022
#  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
"""
import pandas as pd
from pyomo.environ import (
    ConcreteModel, Param, Var, PositiveReals, Objective, Constraint, maximize,
    SolverFactory
)


def reactor_design_model(data):
    
    # Create the concrete model
    model = ConcreteModel()
    
    # Rate constants
    model.k1 = Param(initialize = 5.0/6.0, within=PositiveReals, mutable=True) # min^-1
    model.k2 = Param(initialize = 5.0/3.0, within=PositiveReals, mutable=True) # min^-1
    model.k3 = Param(initialize = 1.0/6000.0, within=PositiveReals, mutable=True) # m^3/(gmol min)

    # Inlet concentration of A, gmol/m^3
    model.caf = Param(initialize = float(data['caf']), within=PositiveReals)
    
	# Space velocity (flowrate/volume)
    model.sv = Param(initialize = float(data['sv']), within=PositiveReals)
    
    # Outlet concentration of each component
    model.ca = Var(initialize = 5000.0, within=PositiveReals) 
    model.cb = Var(initialize = 2000.0, within=PositiveReals) 
    model.cc = Var(initialize = 2000.0, within=PositiveReals) 
    model.cd = Var(initialize = 1000.0, within=PositiveReals)
    
    # Objective
    model.obj = Objective(expr = model.cb, sense=maximize)
    
    # Constraints
    model.ca_bal = 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 = Constraint(expr=(0 == -model.sv * model.cb \
                     + model.k1 * model.ca - model.k2 * model.cb))
    
    model.cc_bal = Constraint(expr=(0 == -model.sv * model.cc \
                     + model.k2 * model.cb))
    
    model.cd_bal = Constraint(expr=(0 == -model.sv * model.cd \
                     + model.k3 * model.ca ** 2.0))
    
    return model

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:
        model = reactor_design_model({'caf': caf, 'sv': sv})
        solver = SolverFactory('ipopt')
        solver.solve(model)
        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-2022
#  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.
#  ___________________________________________________________________________

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 reactor_design_model 


def main():
    # Vars to estimate
    theta_names = ['k1', 'k2', 'k3']
    
    # Data
    file_dirname = dirname(abspath(str(__file__)))
    file_name = abspath(join(file_dirname, 'reactor_data.csv'))
    data = pd.read_csv(file_name) 
    
    # Sum of squared error function
    def SSE(model, data): 
        expr = (float(data['ca']) - model.ca)**2 + \
               (float(data['cb']) - model.cb)**2 + \
               (float(data['cc']) - model.cc)**2 + \
               (float(data['cd']) - model.cd)**2
        return expr
    
    # Create an instance of the parmest estimator
    pest = parmest.Estimator(reactor_design_model, data, theta_names, SSE)
    
    # Parameter estimation
    obj, theta = pest.theta_est()
    
    # Assert statements compare parameter estimation (theta) to an expected value 
    k1_expected = 5.0/6.0
    k2_expected = 5.0/3.0
    k3_expected =  1.0/6000.0
    relative_error = abs(theta['k1'] - k1_expected)/k1_expected
    assert relative_error < 0.05
    relative_error = abs(theta['k2'] - k2_expected)/k2_expected
    assert relative_error < 0.05
    relative_error = abs(theta['k3'] - k3_expected)/k3_expected
    assert relative_error < 0.05
    
if __name__ == "__main__":
    main()
    

The semibatch and Rooney Biegler examples are defined in a similar manner.