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 2017 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, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory

def reactor_design_model(data):
    
    # Create the concrete model
    model = ConcreteModel()
    
    # Rate constants
    model.k1 = Var(initialize = 5.0/6.0, within=PositiveReals) # min^-1
    model.k2 = Var(initialize = 5.0/3.0, within=PositiveReals) # min^-1
    model.k3 = Var(initialize = 1.0/6000.0, within=PositiveReals) # m^3/(gmol min)
    model.k1.fixed = True
    model.k2.fixed = True
    model.k3.fixed = True
    
    # Inlet concentration of A, gmol/m^3
    model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
    model.caf.fixed = True
    
	# Space velocity (flowrate/volume)
    model.sv = Var(initialize = float(data['sv']), within=PositiveReals)
    model.sv.fixed = True
    
    # 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

if __name__ == "__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)
    

The file parmest_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. The file also uses 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 2017 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 numpy as np
import pandas as pd
from itertools import product
import pyomo.contrib.parmest.parmest as parmest
from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model 

### Parameter estimation

# Vars to estimate
theta_names = ['k1', 'k2', 'k3']

# Data
data = pd.read_excel('reactor_data.xlsx') 

# 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

pest = parmest.Estimator(reactor_design_model, data, theta_names, SSE)
obj, theta = pest.theta_est()
print(obj)
print(theta)

### Parameter estimation with bootstrap resampling

bootstrap_theta = pest.theta_est_bootstrap(50)
print(bootstrap_theta.head())

parmest.graphics.pairwise_plot(bootstrap_theta, title='Bootstrap theta estimates')
parmest.graphics.pairwise_plot(bootstrap_theta, theta, 0.8, ['MVN', 'KDE', 'Rect'], 
                      title='Bootstrap theta with confidence regions')

### Likelihood ratio test

k1 = np.arange(0.78, 0.92, 0.02)
k2 = np.arange(1.48, 1.79, 0.05)
k3 = np.arange(0.000155, 0.000185, 0.000005)
theta_vals = pd.DataFrame(list(product(k1, k2, k3)), columns=theta_names)

obj_at_theta = pest.objective_at_theta(theta_vals)
print(obj_at_theta.head())

LR = pest.likelihood_ratio_test(obj_at_theta, obj, [0.8, 0.85, 0.9, 0.95])
print(LR.head())

parmest.graphics.pairwise_plot(LR, theta, 0.8, 
                      title='LR results within 80% confidence region')

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