Parameter Estimation Nrtl Using State Block Solution¶
Parameter Estimation Using the NRTL State Block¶
In this module, we use Pyomo's parmest
tool in conjunction with IDAES models for parameter estimation. We demonstrate these tools by estimating the parameters associated with the NRTL property model for a benzene-toluene mixture. The NRTL model has 2 sets of parameters: the non-randomness parameter (alpha_ij
) and the binary interaction parameter (tau_ij
), where i
and j
are the pure component species. In this example, we only estimate the binary interaction parameter (tau_ij
) for a given dataset. When estimating parameters associated with the property package, IDAES provides the flexibility of doing the parameter estimation by just using the state block or by using a unit model with a specified property package. This module will demonstrate parameter estimation by using only the state block.
We will complete the following tasks:
- Set up a method to return an initialized model
- Set up the parameter estimation problem using
parmest
- Analyze the results
- Demonstrate advanced features using
parmest
Key links to documentation:¶
# Todo: import ConcreteModel from pyomo.environ
from pyomo.environ import ConcreteModel, value
# Todo: import FlowsheetBlock from idaes.core
from idaes.core import FlowsheetBlock
In the next cell, we import the parameter block used in this module and the idaes logger.
from idaes.models.properties.activity_coeff_models.\
BTX_activity_coeff_VLE import BTXParameterBlock
import idaes.logger as idaeslog
In the next cell, we import parmest
from Pyomo and the pandas
package. We need pandas
as parmest
uses pandas.dataframe
for handling the input data and the results.
import pyomo.contrib.parmest.parmest as parmest
import pandas as pd
Setting up an initialized model¶
We need to provide a method that returns an initialized model to the parmest
tool in Pyomo.
def NRTL_model(data):
#Todo: Create a ConcreteModel object
m = ConcreteModel()
#Todo: Create FlowsheetBlock object
m.fs = FlowsheetBlock(dynamic=False)
#Todo: Create a properties parameter object with the following options:
# "valid_phase": ('Liq', 'Vap')
# "activity_coeff_model": 'NRTL'
m.fs.properties = BTXParameterBlock(valid_phase=('Liq', 'Vap'),
activity_coeff_model='NRTL')
m.fs.state_block = m.fs.properties.build_state_block(
defined_state=True)
# Fix the state variables on the state block
# hint: state variables exist on the state block i.e. on m.fs.state_block
m.fs.state_block.flow_mol.fix(1)
m.fs.state_block.temperature.fix(368)
m.fs.state_block.pressure.fix(101325)
m.fs.state_block.mole_frac_comp["benzene"].fix(0.5)
m.fs.state_block.mole_frac_comp["toluene"].fix(0.5)
# Fix NRTL specific parameters.
# non-randomness parameter - alpha_ij (set at 0.3, 0 if i=j)
m.fs.properties.\
alpha["benzene", "benzene"].fix(0)
m.fs.properties.\
alpha["benzene", "toluene"].fix(0.3)
m.fs.properties.\
alpha["toluene", "toluene"].fix(0)
m.fs.properties.\
alpha["toluene", "benzene"].fix(0.3)
# binary interaction parameter - tau_ij (0 if i=j, else to be estimated later but fixing to initialize)
m.fs.properties.\
tau["benzene", "benzene"].fix(0)
m.fs.properties.\
tau["benzene", "toluene"].fix(-0.9)
m.fs.properties.\
tau["toluene", "toluene"].fix(0)
m.fs.properties.\
tau["toluene", "benzene"].fix(1.4)
# Initialize the flash unit
m.fs.state_block.initialize(outlvl=idaeslog.INFO_LOW)
# Fix at actual temperature
m.fs.state_block.temperature.fix(float(data["temperature"]))
# Set bounds on variables to be estimated
m.fs.properties.\
tau["benzene", "toluene"].setlb(-5)
m.fs.properties.\
tau["benzene", "toluene"].setub(5)
m.fs.properties.\
tau["toluene", "benzene"].setlb(-5)
m.fs.properties.\
tau["toluene", "benzene"].setub(5)
# Return initialized flash model
return m
Parameter estimation using parmest¶
In addition to providing a method to return an initialized model, the parmest
tool needs the following:
- List of variable names to be estimated
- Dataset
- Expression to compute the sum of squared errors
In this example, we only estimate the binary interaction parameter (tau_ij
). Given that this variable is usually indexed as tau_ij = Var(component_list, component_list)
, there are 2*2=4 degrees of freedom. However, when i=j, the binary interaction parameter is 0. Therefore, in this problem, we estimate the binary interaction parameter for the following variables only:
- fs.properties.tau['benzene', 'toluene']
- fs.properties.tau['toluene', 'benzene']
# Todo: Create a list of vars to estimate
variable_name = ["fs.properties.tau['benzene', 'toluene']",
"fs.properties.tau['toluene', 'benzene']"]
Pyomo's parmest
tool supports the following data formats:
- pandas dataframe
- list of dictionaries
- list of json file names.
Please see the documentation for more details.
For this example, we load data from the csv file BT_NRTL_dataset.csv
. The dataset consists of fifty data points which provide the mole fraction of benzene in the vapor and liquid phase as a function of temperature.
# Load data from csv
data = pd.read_csv('BT_NRTL_dataset.csv')
# Display the dataset
display(data)
temperature | liq_benzene | vap_benzene | |
---|---|---|---|
0 | 365.500000 | 0.480953 | 0.692110 |
1 | 365.617647 | 0.462444 | 0.667699 |
2 | 365.735294 | 0.477984 | 0.692441 |
3 | 365.852941 | 0.440547 | 0.640336 |
4 | 365.970588 | 0.427421 | 0.623328 |
5 | 366.088235 | 0.442725 | 0.647796 |
6 | 366.205882 | 0.434374 | 0.637691 |
7 | 366.323529 | 0.444642 | 0.654933 |
8 | 366.441176 | 0.427132 | 0.631229 |
9 | 366.558824 | 0.446301 | 0.661743 |
10 | 366.676471 | 0.438004 | 0.651591 |
11 | 366.794118 | 0.425320 | 0.634814 |
12 | 366.911765 | 0.439435 | 0.658047 |
13 | 367.029412 | 0.435655 | 0.654539 |
14 | 367.147059 | 0.401350 | 0.604987 |
15 | 367.264706 | 0.397862 | 0.601703 |
16 | 367.382353 | 0.415821 | 0.630930 |
17 | 367.500000 | 0.420667 | 0.640380 |
18 | 367.617647 | 0.391683 | 0.598214 |
19 | 367.735294 | 0.404903 | 0.620432 |
20 | 367.852941 | 0.409563 | 0.629626 |
21 | 367.970588 | 0.389488 | 0.600722 |
22 | 368.000000 | 0.396789 | 0.612483 |
23 | 368.088235 | 0.398162 | 0.616106 |
24 | 368.205882 | 0.362340 | 0.562505 |
25 | 368.323529 | 0.386958 | 0.602680 |
26 | 368.441176 | 0.363643 | 0.568210 |
27 | 368.558824 | 0.368118 | 0.577072 |
28 | 368.676471 | 0.384098 | 0.604078 |
29 | 368.794118 | 0.353605 | 0.557925 |
30 | 368.911765 | 0.346474 | 0.548445 |
31 | 369.029412 | 0.350741 | 0.556996 |
32 | 369.147059 | 0.362347 | 0.577286 |
33 | 369.264706 | 0.362578 | 0.579519 |
34 | 369.382353 | 0.340765 | 0.546411 |
35 | 369.500000 | 0.337462 | 0.542857 |
36 | 369.617647 | 0.355729 | 0.574083 |
37 | 369.735294 | 0.348679 | 0.564513 |
38 | 369.852941 | 0.338187 | 0.549284 |
39 | 369.970588 | 0.324360 | 0.528514 |
40 | 370.088235 | 0.310753 | 0.507964 |
41 | 370.205882 | 0.311037 | 0.510055 |
42 | 370.323529 | 0.311263 | 0.512055 |
43 | 370.441176 | 0.308081 | 0.508437 |
44 | 370.558824 | 0.308224 | 0.510293 |
45 | 370.676471 | 0.318148 | 0.528399 |
46 | 370.794118 | 0.308334 | 0.513728 |
47 | 370.911765 | 0.317937 | 0.531410 |
48 | 371.029412 | 0.289149 | 0.484824 |
49 | 371.147059 | 0.298637 | 0.502318 |
We need to provide a method to return an expression to compute the sum of squared errors that will be used as the objective in solving the parameter estimation problem. For this problem, the error will be computed for the mole fraction of benzene in the vapor and liquid phase between the model prediction and data.
# Create method to return an expression that computes the sum of squared error
def SSE(m, data):
# Todo: Add expression for computing the sum of squared errors in mole fraction of benzene in the liquid
# and vapor phase. For example, the squared error for the vapor phase is:
# (float(data["vap_benzene"]) - m.fs.state_block.mole_frac_phase_comp["Vap", "benzene"])**2
expr = ((float(data["vap_benzene"]) -
m.fs.state_block.mole_frac_phase_comp["Vap", "benzene"])**2 +
(float(data["liq_benzene"]) -
m.fs.state_block.mole_frac_phase_comp["Liq", "benzene"])**2)
return expr*1E4
We are now ready to set up the parameter estimation problem. We will create a parameter estimation object called pest
. As shown below, we pass the method that returns an initialized model, data, variable_name, and the SSE expression to the Estimator method. tee=True
will print the solver output after solving the parameter estimation problem.
import logging; idaeslog.getIdaesLogger("core.property_meta").setLevel(logging.ERROR)
# Initialize a parameter estimation object
pest = parmest.Estimator(NRTL_model, data, variable_name, SSE, tee=True)
# Run parameter estimation using all data
obj_value, parameters = pest.theta_est()
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...: 3746 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 2200 Total number of variables............................: 1100 variables with only lower bounds: 0 variables with lower and upper bounds: 300 variables with only upper bounds: 0 Total number of equality constraints.................: 1098 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 6.0671019e+01 3.15e+00 4.84e+01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 5.2961050e+00 1.76e+03 5.05e+00 -1.0 1.37e+04 - 9.74e-01 1.00e+00f 1 2 5.2586169e+00 4.01e+02 1.09e+00 -1.0 5.15e+02 - 1.00e+00 1.00e+00h 1 3 5.1450958e+00 7.04e+01 2.27e-01 -1.0 4.11e+01 - 1.00e+00 1.00e+00h 1 4 5.0748980e+00 1.25e+02 2.08e-01 -1.7 5.74e+02 - 1.00e+00 1.00e+00h 1 5 5.0775194e+00 7.87e+00 1.92e-01 -1.7 8.44e+01 - 1.00e+00 1.00e+00h 1 6 5.0726692e+00 1.37e+01 1.90e-01 -2.5 1.38e+02 - 1.00e+00 1.00e+00h 1 7 5.0750377e+00 2.85e+00 2.60e-02 -2.5 6.99e+01 - 1.00e+00 1.00e+00h 1 8 5.0749670e+00 7.36e-02 2.81e-03 -3.8 9.72e+00 - 1.00e+00 1.00e+00h 1 9 5.0749687e+00 4.51e-04 4.80e-06 -3.8 1.01e+00 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 5.0749686e+00 2.91e-04 1.36e-06 -5.7 5.81e-01 - 1.00e+00 1.00e+00h 1 11 5.0749686e+00 4.78e-08 2.18e-10 -8.6 7.65e-03 - 1.00e+00 1.00e+00h 1 Number of Iterations....: 11 (scaled) (unscaled) Objective...............: 5.0749685783046434e+00 5.0749685783046434e+00 Dual infeasibility......: 2.1827409324437497e-10 2.1827409324437497e-10 Constraint violation....: 1.6625508263665860e-10 4.7832145355641842e-08 Complementarity.........: 2.5076274461651402e-09 2.5076274461651402e-09 Overall NLP error.......: 2.5076274461651402e-09 4.7832145355641842e-08 Number of objective function evaluations = 12 Number of objective gradient evaluations = 12 Number of equality constraint evaluations = 12 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 12 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 11 Total CPU secs in IPOPT (w/o function evaluations) = 0.029 Total CPU secs in NLP function evaluations = 0.009 EXIT: Optimal Solution Found.
You will notice that the resulting parameter estimation problem will have 1102 variables and 1100 constraints. Let us display the results by running the next cell.
print("The SSE at the optimal solution is %0.6f" % (obj_value*1e-4))
print()
print("The values for the parameters are as follows:")
for k,v in parameters.items():
print(k, "=", v)
The SSE at the optimal solution is 0.000507 The values for the parameters are as follows: fs.properties.tau[benzene,toluene] = -0.8987624036283798 fs.properties.tau[toluene,benzene] = 1.4104861099366137
Using the data that was provided, we have estimated the binary interaction parameters in the NRTL model for a benzene-toluene mixture. Although the dataset that was provided was temperature dependent, in this example we have estimated a single value that fits best for all temperatures.
Advanced options for parmest: bootstrapping¶
Pyomo's parmest
tool allows for bootstrapping where the parameter estimation is repeated over n
samples with resampling from the original data set. Parameter estimation with bootstrap resampling can be used to identify confidence regions around each parameter estimate. This analysis can be slow given the increased number of model instances that need to be solved. Please refer to https://pyomo.readthedocs.io/en/stable/contributed_packages/parmest/driver.html for more details.
For the example above, the bootstrapping can be run by uncommenting the code in the following cell:
# Run parameter estimation using bootstrap resample of the data (10 samples),
# plot results along with confidence regions
# Uncomment the following code:
# bootstrap_theta = pest.theta_est_bootstrap(4)
# display(bootstrap_theta)