5.1. Parameter estimation with parmest#
Created by Kanishka Ghosh, Jialu Wang, and Prof. Alex Dowling at the University of Notre Dame.
# This code cell installs packages on Colab
import sys
if "google.colab" in sys.modules:
!wget "https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py"
import helper
helper.install_idaes()
helper.install_ipopt()
helper.download_data(['parmest_20210609_data_exp{:d}.csv'.format(i) for i in range(1,17)])
helper.download_data(['parmest_log_file.csv'])
!pyomo build-extensions
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pyomo.environ import *
from pyomo.dae import *
# Define the directory to save/read the data files
data_dir = https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/'
5.1.1. What is parameter estimation?#
Given a function \(f(x,\theta)\) where \(x\) is the input or array of inputs, \(\theta\) is the vector of unknown model parameters, and \(y\) is the array of observed output, parameter estimation is performed to determine the values of \(\theta\) to minimize the error between \(f(x,\theta)\) and \(y\). Commonly, parameter estimation is set up as a least squares objective problem:
\[\begin{split}
\begin{align}
\begin{split}
\min_{\hat{\theta}} \quad & \sum_{i}^{} (y_i - f(x_i,\hat{\theta}))^2\\
\textrm{s.t.} \quad & \mathrm{bounds \ on} \ \theta\\
& \mathrm{other \ physical \ constraints}\\
\end{split}
\end{align}
\end{split}\]
where \(i\) is used to index the datapoints in a dataset and \(\hat{\theta}\) is the optimal set of parameter values that minimizes the prediction error.
5.1.2. What is parmest?#
parmest is a Python package built on the Pyomo optimization modeling language to support parameter estimation using experimental data along with confidence regions and subsequent creation of scenarios for PySP. parmest supports scenario generation for multiple ‘experiments’ and can be used to characterize estimate uncertainties through, for example, confidence region generations. parmest requires the following positional arguments in order solve the optimization problem:
Function that accepts an ‘experimental’ dataset or a list of ‘experimental’ datasets, each defined as a dictionary, as it’s argument and returns the Pyomo model.
Later in this tutorial, that function is defined above as
create_model()List of datasets where each dataset is a dictionary
Later in this tutorial, the list of datasets is generated using the function
create_data_dict()(defined below) and is stored indata_dict_overallList of parameter names (as they appear in the Pyomo model definition) that are being estimated
later in this tutorial, the list of parameter names to be estimated is defined below by
theta_namesOptional keyword argument to define the verbosity of solver output. Default: False
More information about the parmest package can be found here.
Detailed explanation of the various methods in parmest can be found here.
5.1.3. Example: Reaction Kinetics#
Consider two chemical reactions that convert molecule \(A\) to desired product \(B\) and a less valuable side-product \(C\).
\[A \overset{k_1}{\rightarrow} B \overset{k_2}{\rightarrow} C\]
Our ultimate goals is to design a large-scale continous reactor that maximizes the production of \(B\). This general sequential reactions problem is widely applicable to CO\(_2\) capture and industry more broadly (petrochemicals, pharmasuticals, etc.).
The rate laws for these two chemical reactions are:
\[r_A = -k_1 C_A\]
\[r_B = k_1 C_A - k_2 C_B\]
\[r_C = k_2 C_B\]
Here, \(C_A\), \(C_B\), and \(C_C\) are the concentrations of each species. The rate constants \(k_1\) and \(k_2\) depend on temperature as follows:
\[k_1 = A_1 \exp{\frac{-E_1}{R T}}\]
\[k_2 = A_2 \exp{\frac{-E_2}{R T}}\]
\(A_1, A_2, E_1\), and \(E_2\) are fitted model parameters. \(R\) is the ideal-gas constant and \(T\) is absolute temperature.
5.1.3.1. Batch Reactor#
The concentrations in a batch reactor evolve with time per the following differential equations:
\[ \frac{d C_A}{dt} = r_A = -k_1 C_A \]
\[ \frac{d C_B}{dt} = r_B = k_1 C_A - k_2 C_B \]
\[ \frac{d C_C}{dt} = r_C = k_2 C_B \]
This is a linear system of differential equations. Assuming the feed is only species \(A\), i.e.,
\[C_A(t=0) = C_{A0} \quad C_B(t=0) = 0 \quad C_C(t=0) = 0\]
leads to the following analytic solution:
\[C_A(t) = C_{A,0} \exp(-k_1 t)\]
\[C_B(t) = \frac{k_1}{k_2 - k_1} C_{A,0} \left[\exp(-k_1 t) - \exp(-k_2 t) \right]\]
\[C_C(t) = C_{A,0} - \frac{k_2}{k_2 - k_1} C_{A,0} \exp(-k_1 t) + \frac{k_1}{k_2 - k_1} \exp(-k_2 t) C_{A,0} = C_{A,0} - C_{A}(t) - C_{B}(t)\]
The following Python code simulates and plots this model.
def kinetics(A, E, T):
''' Computes kinetics from Arrhenius equation
Arguments:
A: pre-exponential factor, [1 / hr]
E: activation energy, [kJ / mol]
T: temperature, [K]
Returns:
k: reaction rate coefficient, [1 / hr]
'''
R = 8.31446261815324 # J / K / mole
return A * np.exp(-E*1000/(R*T))
def concentrations(t,k,CA0):
'''
Returns concentrations at time t
Arguments:
t: time, [hr]
k: reaction rate coefficient, [1 / hr]
CA0: initial concentration of A, [mol / L]
Returns:
CA, CB, CC: concentrations of A, B, and C at time t, [mol / L]
'''
CA = CA0 * np.exp(-k[0]*t);
CB = k[0]*CA0/(k[1]-k[0]) * (np.exp(-k[0]*t) - np.exp(-k[1]*t));
CC = CA0 - CA - CB;
return CA, CB, CC
CA0 = 1 # Moles/L
k = [3, 0.7] # 1/hr
t = np.linspace(0,1,51)
CA, CB, CC = concentrations(t,k,CA0)
plt.plot(t, CA, label="$C_{A}$",linestyle="-",color="blue")
plt.plot(t, CB, label="$C_{B}$",linestyle="-.",color="green")
plt.plot(t, CC, label="$C_{C}$",linestyle="--",color="red")
plt.xlabel("Time [hours]")
plt.ylabel("Concentration [mol/L]")
plt.title("Batch Reactor Model")
plt.legend()
plt.show()
plt.close()
5.1.3.2. Experimental Data#
See the notebook Supplementary material: data for parmest tutorial for details on how these experimental data were generated (via simulation).
Experimental data consists of the concentration of species A, B, and C in \(\mathrm{mol/L}\) with respect to time \(t\) in \(\mathrm{hours}\) inside the batch reactor. The experimental data is stored in csv files where the first column records the time \(t\) in the reactor. Next, the temperature \(T\) in \(\mathrm{K}\) at which the reaction was simulated is recorded followed by the initial concentration of species A, \(C_{A,0}\), in \(\mathrm{mol/L}\). Finally, the time-varying species concentrations (\(C_A\)), (\(C_B\)), and (\(C_C\)) are recorded in \(\mathrm{mol/L}\). Following is how the pandas dataframe of a single experiment looks like:
# define function to plot
def plot_exp(k, CA0, data, text):
'''
Plot concentration profiles
Arguments:
k: kinetic parameters
CA0: initial concentration
data: Pandas data frame
text: plot title
'''
# evaluate models
t = np.linspace(0,1,51)
CA, CB, CC = concentrations(t,k,CA0)
# plot model-generated and 'experimental' data
# symbols for 'experimental' data
# solid and dashed lines for model-generated data
plt.plot(t, CA,label="$C_{A}$",linestyle="-",color="blue")
plt.plot(data.time, data.CA, marker='o',linestyle="",color="blue",label=str())
plt.plot(t, CB, label="$C_{B}$",linestyle="-.",color="green")
plt.plot(data.time, data.CB, marker='s',linestyle="",color="green",label=str())
plt.plot(t, CC, label="$C_{C}$",linestyle="--",color="red")
plt.plot(data.time, data.CC, marker='^',linestyle="",color="red",label=str())
plt.xlabel("Time [hours]")
plt.ylabel("Concentration [mol/L]")
plt.title(text)
plt.legend()
plt.show()
plt.close()
5.1.3.3. Pyomo model#
In the following cell, we define a function to define and return the Pyomo model for the kinetic model to be used for parameter estimation.
def create_model(data):
'''
function to create Pyomo model
Argument:
data: a single dictionary of data
Return:
m: Pyomo model
'''
# data
exp_data = data['data']
# This code style matches parmest example found here:
# https://github.com/Pyomo/pyomo/blob/master/pyomo/contrib/parmest/examples/semibatch/semibatch.py
# unpack 'experimental' data into temporary variables
cameastemp = exp_data['CA']
cbmeastemp = exp_data['CB']
ccmeastemp = exp_data['CC']
tmeastemp = exp_data['time']
# create dictionaries for 'experimental' data of CA, CB,and CC indexed by timestep
cameas={}
cbmeas={}
ccmeas={}
for i,j in enumerate(tmeastemp):
cameas[float(j)] = cameastemp[i]
cbmeas[float(j)] = cbmeastemp[i]
ccmeas[float(j)] = ccmeastemp[i]
# define Pyomo model
m = ConcreteModel()
m.T = data['T'] # K
m.CA0 = data['CA0'] # mol/L
# define 'experimental' data timesteps as Pyomo set
m.t = Set(initialize=tmeastemp.tolist())
# define 'experimental' data as Pyomo parameters indexed by timestep set and
# initialized by dictionary of experimental data
m.Ca_meas = Param(m.t, initialize=cameas)
m.Cb_meas = Param(m.t, initialize=cbmeas)
m.Cc_meas = Param(m.t, initialize=ccmeas)
m.R = 8.31446261815324 # J / K / mole
# Kinetic parameters to be fitted defined as Pyomo variables
# Initialized by 'true' values
m.A1 = Var(initialize=200, bounds=(100,300)) # 1/hr
m.A2 = Var(initialize=400, bounds=(300,500)) # 1/hr
m.E1 = Var(initialize=10, bounds=(1,20)) # kJ/mol
m.E2 = Var(initialize=15, bounds=(1,30)) # kJ/mol
# Concentration variables indexed by time
m.CA = Var(m.t, initialize = m.CA0) # mol/L
m.CB = Var(m.t, initialize = 0) # mol/L
m.CC = Var(m.t, initialize = 0) # mol/L
# kinetic rate constants from Arrhenius equation
m.k1 = Expression(rule = m.A1 * exp(-m.E1*1000/(m.R*m.T))) # 1/hr
m.k2 = Expression(rule = m.A2 * exp(-m.E2*1000/(m.R*m.T))) # 1/hr
# Constraints to change concentrations based on kinetics
def conc_A(m,i):
if i == 0:
return Constraint.Skip
else:
return m.CA[i] == m.CA0 * exp(-m.k1*i)
m.CA_rate = Constraint(m.t,rule=conc_A)
def conc_B(m,i):
if i == 0:
return Constraint.Skip
else:
return m.CB[i] == m.k1*m.CA0/(m.k2-m.k1) * (exp(-m.k1*i) - exp(-m.k2*i))
m.CB_rate = Constraint(m.t,rule=conc_B)
def conc_C(m,i):
if i == 0:
return Constraint.Skip
else:
return m.CC[i] == m.CA0 - m.CA[i] - m.CB[i]
m.CC_rate = Constraint(m.t,rule=conc_C)
# Initial Conditions
def _initcon(m):
yield m.CA[m.t.first()] == m.CA0
yield m.CB[m.t.first()] == 0.0
yield m.CC[m.t.first()] == 0.0
m.initcon = ConstraintList(rule=_initcon)
# Objective function
# The objective function for parmest is defined as a 2-stage stochastic optimization objective function
# First stage cost: independent of scenarios ('experiments')
# expression for minimizing fixed realization
# from model. Eg.: reactor temperature, size, etc.
def ComputeFirstStageCost_rule(m):
# In this case, we do not optimize anything besides the kinetic parameters through
# least square fitting realizations at each timestep defined by m.t.
# Hence, the first stage cost is set to 0 here.
return 0
m.FirstStageCost = Expression(rule=ComputeFirstStageCost_rule)
# Second stage cost: Realization at each scenario over which the model is defined
def ComputeSecondStageCost_rule(m):
# In this problem, we want to minimize the sum of squared errors between
# 'experimental' data and the model realization of concentrations of
# A, B, and C over each scenario (here, timesteps defined by m.t)
return sum((m.CA[t] - m.Ca_meas[t]) ** 2 + (m.CB[t] - m.Cb_meas[t]) ** 2
+ (m.CC[t] - m.Cc_meas[t]) ** 2 for t in m.t)
m.SecondStageCost = Expression(rule=ComputeSecondStageCost_rule)
# return the sum of the first-stage and second-stage costs as the objective function
def total_cost_rule(m):
return m.FirstStageCost + m.SecondStageCost
m.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)
return m
5.1.4. Parameter estimation with a single dataset#
Here, we will estimate parameters \(A_1\), \(A_2\), \(E_1\), and \(E_2\) using data generated for a batch ‘experiment’ at 250 C with an inlet concentration of 0.5 mol/L of A.
The parameter estimation problem is solved with the least squares optimization scheme
The parameter estimation problem is solved with the least squares optimization scheme $$
(5.1)#\[\begin{align}
\begin{split}
\min_{\hat{\theta}} \quad & \sum_{i}^{} (y_i - f(x_i,\hat{\theta}))^2\\
\textrm{s.t.} \quad & \mathrm{bounds \ on} \ \theta\\
& \mathrm{other \ physical \ constraints}\\
\end{split}
\end{align}\]
$\( where \)i\( is used to index the datapoints in a dataset and \)\hat{\theta}$ is the optimal set of parameter values that minimizes the prediction error.
# read-in data from csv file
data = pd.read_csv(https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/parmest_20210609_data_exp1.csv',index_col=0)
# create dictionary of data
# format of dictionary needs to be consistent with input provided to create_model() function
data_dict = {}
data_dict['T'] = data['T'].iloc[0]
data_dict['CA0'] = data['CA0'].iloc[0]
data = data.drop(labels=['T','CA0'],axis=1)
data_dict['data'] = data
# create pyomo model instance
model = create_model(data_dict)
# create solver instance for IPOPT
solver = SolverFactory('ipopt')
# solve model
solver.solve(model,tee=True)
print("=== Parameter values ===")
print("A1 = {:0.3f} 1/hr".format(value(model.A1)))
print("A2 = {:0.3f} 1/hr".format(value(model.A2)))
print("E1 = {:0.3f} kJ/mol".format(value(model.E1)))
print("E2 = {:0.3f} kJ/mol".format(value(model.E2)))
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...: 91
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 37
Total number of variables............................: 31
variables with only lower bounds: 0
variables with lower and upper bounds: 4
variables with only upper bounds: 0
Total number of equality constraints.................: 27
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 1.6338341e+00 4.02e-01 6.37e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.9026209e-01 2.34e-04 4.98e-03 -1.0 4.13e-01 - 9.86e-01 1.00e+00f 1
2 1.8657853e-01 5.60e-04 2.39e-04 -1.0 1.99e-01 - 1.00e+00 1.00e+00h 1
3 1.8638673e-01 4.98e-05 1.92e-05 -2.5 7.10e-02 - 1.00e+00 1.00e+00h 1
4 1.8638598e-01 2.10e-07 1.36e-07 -3.8 4.52e-02 - 1.00e+00 1.00e+00h 1
5 1.8638599e-01 2.23e-10 3.93e-11 -5.7 9.66e-03 - 1.00e+00 1.00e+00h 1
6 1.8638599e-01 2.55e-12 2.09e-13 -8.6 1.05e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 6
(scaled) (unscaled)
Objective...............: 1.8638598612196314e-01 1.8638598612196314e-01
Dual infeasibility......: 2.0857837216213979e-13 2.0857837216213979e-13
Constraint violation....: 2.5491275756905907e-12 2.5491275756905907e-12
Complementarity.........: 2.5255218786493285e-09 2.5255218786493285e-09
Overall NLP error.......: 2.5255218786493285e-09 2.5255218786493285e-09
Number of objective function evaluations = 7
Number of objective gradient evaluations = 7
Number of equality constraint evaluations = 7
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 7
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 6
Total CPU secs in IPOPT (w/o function evaluations) = 0.003
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
=== Parameter values ===
A1 = 200.209 1/hr
A2 = 400.037 1/hr
E1 = 9.640 kJ/mol
E2 = 14.939 kJ/mol
5.1.5. Parameter estimation with multiple datasets#
Here, we will estimate parameters \(A_1\), \(A_2\), \(E_1\), and \(E_2\) using data generated for a batch ‘experiment’ at 250 C with an inlet concentration of 0.5 mol/L of A.
The parameter estimation problem is now defined to solve optimization problem where the objective function is the mean of the least square error between observed and calculated data for multiple experiments.
\[\begin{split}
\begin{align}
\begin{split}
\min_{\hat{\theta}} \quad & \frac{1}{N}\sum_{j}^{}\sum_{i}^{} (y_{j,i} - f(x_{j,i},\hat{\theta}))^2\\
\textrm{s.t.} \quad & \mathrm{bounds \ on} \ \theta\\
& \mathrm{other \ physical \ constraints}\\
\end{split}
\end{align}
\end{split}\]
where \(i\) is used to index the datapoints in a dataset, \(j\) is an index on the dataset such that \(j \in [1,N]\) and \(N\) is the number of experiments conducted.
5.1.5.1. Generate list of dataset#
In the following cell, we define a function to generate a list of dictionaries containing the ‘experimental’ data. For this, we read-in the list of file names generated earlier
def create_data_dict(files):
'''
Create a list of dictionaries from multiple datasets
Arguments:
files: pandas dataframe of file names
Return:
data_dict_list: list of dictionaries
'''
data_dict_list = []
for index, file in files.iterrows():
# create a dictionary of 'experimental' data for temperature T and initial concentration of A CA0
data_dict = {}
data = pd.read_csv(str(file.values[0]),index_col=0)
data_dict['T'] = data['T'].iloc[0]
data_dict['CA0'] = data['CA0'].iloc[0]
data = data.drop(labels=['T','CA0'],axis=1)
data_dict['data'] = data
# add dictionary to list to be return
data_dict_list.append(data_dict)
return data_dict_list
# list of temperatures
T_vals = [250,300,350,400] # K
# list of initial concentrations of A
CA0_vals = [0.5,1.0,1.5,2.0] # mol/L
# read-in file names from log file created using gen_data()
file_list = pd.read_csv(data_dir + 'parmest_log_file.csv')
# create a list with all experiments stored separately as dictionaries
data_dict_overall = create_data_dict(file_list)
# printing a dictionary of data for illustrative purposes
print(data_dict_overall[0])
# Plotting a dataset for illustrative purposes
plt.figure()
plt.plot(data_dict_overall[0]['data']['time'], data_dict_overall[0]['data']['CA'], label="$C_{A}$",marker="o",linestyle='',color="blue")
plt.plot(data_dict_overall[0]['data']['time'], data_dict_overall[0]['data']['CB'], label="$C_{B}$",marker="s",linestyle='',color="green")
plt.plot(data_dict_overall[0]['data']['time'], data_dict_overall[0]['data']['CC'], label="$C_{C}$",marker="^",linestyle='',color="red")
plt.xlabel("Time [hours]")
plt.ylabel("Concentration [mol/L]")
plt.title("Batch Reactor Model")
plt.legend()
plt.show()
{'T': 250, 'CA0': 0.5, 'data': time CA CB CC
0 0.000 0.476266 0.000000 0.000000
1 0.125 0.332309 0.133593 0.000000
2 0.250 0.295725 0.362016 0.092655
3 0.375 0.394806 0.131512 0.004404
4 0.500 0.167992 0.271682 0.031808
5 0.625 0.117872 0.313196 0.009108
6 0.750 0.117052 0.582722 0.000000
7 0.875 0.175373 0.224904 0.043844
8 1.000 0.130619 0.360287 0.194820}
5.1.5.2. Parameter estimation with parmest#
In the following cell, we perform parameter estimation using parmest to solve the least squares problem defined in the Pyomo model.
# import parmest
import pyomo.contrib.parmest.parmest as parmest
# defining the names of the parameters in a list
theta_names = ['A1','A2','E1','E2']
# create an object using parmest.Estimator() that stores the Pyomo model realizations for the datasets provided.
# This object which will be used to determined the parameter values that best fit all the datasets
pest = parmest.Estimator(create_model,data_dict_overall,theta_names,tee=True)
# call the method theta_est() for the Estimator() object defined above to solve
# the parameter estimation problem.
# theta_est() returns:
# the overall objective function value
# estimated parameter values (dictionary with keys = parameters names as defined in the Pyomo model)
obj, theta = pest.theta_est()
print('theta:\n',theta)
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...: 1576
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 592
Total number of variables............................: 496
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 492
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 1.7039989e+01 2.00e+00 1.40e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 2.2477764e-01 4.80e-06 1.02e-03 -1.0 2.00e+00 - 9.91e-01 1.00e+00f 1
2 2.2234827e-01 3.02e-04 7.58e-05 -1.7 1.03e-01 - 1.00e+00 1.00e+00h 1
3 2.2228920e-01 7.67e-06 2.49e-06 -3.8 3.78e-01 - 1.00e+00 1.00e+00h 1
4 2.2211618e-01 1.39e-03 1.28e-05 -5.7 1.20e+01 - 9.13e-01 1.00e+00h 1
5 2.2210772e-01 3.75e-05 3.59e-07 -5.7 1.88e+00 - 1.00e+00 1.00e+00h 1
6 2.2210762e-01 7.16e-08 5.59e-10 -5.7 1.77e-01 - 1.00e+00 1.00e+00h 1
7 2.2210762e-01 2.40e-08 2.20e-10 -8.6 4.74e-02 - 1.00e+00 1.00e+00h 1
8 2.2210762e-01 3.00e-12 2.31e-14 -9.0 1.17e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 8
(scaled) (unscaled)
Objective...............: 2.2210762190708977e-01 2.2210762190708977e-01
Dual infeasibility......: 2.3121945602686378e-14 2.3121945602686378e-14
Constraint violation....: 3.0038194154258235e-12 3.0038194154258235e-12
Complementarity.........: 9.0911913344645786e-10 9.0911913344645786e-10
Overall NLP error.......: 9.0911913344645786e-10 9.0911913344645786e-10
Number of objective function evaluations = 9
Number of objective gradient evaluations = 9
Number of equality constraint evaluations = 9
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 9
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 8
Total CPU secs in IPOPT (w/o function evaluations) = 0.004
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
theta:
{'A1': 185.6087678995809, 'A2': 401.1702352092697, 'E1': 9.866878463449424, 'E2': 14.866030991977437}
5.1.5.3. Plotting fitted model simulation with ‘experimental’ data#
Next, we plot the ‘experimental’ data along with the profiles generated using the fitted kinetic model. The symbols represent the ‘experimental’ data and the solid and dashed lines are the profiles generated using the fitted model.
# Parameter values from parameter estimation using parmest
A1 = theta['A1']
E1 = theta['E1']
A2 = theta['A2']
E2 = theta['E2']
A_est1 = [A1, A2]
A_est = np.asarray(A_est1)
E_est1 = [E1, E2]
E_est = np.asarray(E_est1)
ctr = 0
for T in T_vals:
for CA0 in CA0_vals:
# generate concentration profiles using estimated parameter values
k = kinetics(A_est, E_est, T)
# plot model-generated and 'experimental' data
# symbols for 'experimental' data
# solid and dashed lines for model-generated data
plot_exp(k, CA0, data_dict_overall[ctr]['data'], 'Model prediction and experimental value at T = {} K and $C_{}$ = {} mol/L'.format(T,'A0',CA0))
ctr+=1
5.1.6. Using parmest with pyomo.dae#
In contrast to the approach above, we will now try to solve the model without the analytic solution for the concentrations using Pyomo.DAE. To recap, the concentrations in a batch reactor evolve with time per the following differential equations:
\[ \frac{d C_A}{dt} = r_A = -k_1 C_A \]
\[ \frac{d C_B}{dt} = r_B = k_1 C_A - k_2 C_B \]
\[ \frac{d C_C}{dt} = r_C = k_2 C_B \]
This is a linear system of differential equations. Assuming the feed is only species \(A\), i.e.,
\[C_A(t=0) = C_{A0} \quad C_B(t=0) = 0 \quad C_C(t=0) = 0\]
In the following cell, we define a function to define and return the Pyomo DAE model (dynamic mode) for the kinetic model to be used for parameter estimation. In this model, the rate equations are presented in terms of linear differential equations.
def create_model_DAE(data):
'''
function to create Pyomo model
Argument:
data: a single dictionary of data
Return:
m: Pyomo model
'''
# data
exp_data = data['data']
# This code style matches parmest example found here:
# https://github.com/Pyomo/pyomo/blob/master/pyomo/contrib/parmest/examples/semibatch/semibatch.py
# unpack 'experimental' data into temporary variables
cameastemp = exp_data['CA']
cbmeastemp = exp_data['CB']
ccmeastemp = exp_data['CC']
tmeastemp = exp_data['time']
# create dictionaries for 'experimental' data of CA, CB,and CC indexed by timestep
cameas={}
cbmeas={}
ccmeas={}
for i,j in enumerate(tmeastemp):
cameas[float(j)] = cameastemp[i]
cbmeas[float(j)] = cbmeastemp[i]
ccmeas[float(j)] = ccmeastemp[i]
# define Pyomo model
m = ConcreteModel()
m.T = data['T'] # K
m.CA0 = data['CA0'] # mol/L
# define 'experimental' data timesteps as Pyomo ContinuousSet
m.t = ContinuousSet(bounds = (0.0, tmeastemp.iloc[-1]), initialize=tmeastemp.tolist())
# define 'experimental' data as Pyomo parameters indexed by timestep set and
# initialized by dictionary of experimental data
m.Ca_meas = Param(m.t, initialize=cameas)
m.Cb_meas = Param(m.t, initialize=cbmeas)
m.Cc_meas = Param(m.t, initialize=ccmeas)
m.R = 8.31446261815324 # J / K / mole
# Kinetic parameters to be fitted defined as Pyomo variables
# Initialized by 'true' values
m.A1 = Var(initialize=200, bounds=(100,300)) # 1/hr
m.A2 = Var(initialize=400, bounds=(300,500)) # 1/hr
m.E1 = Var(initialize=10, bounds=(1,20)) # kJ/mol
m.E2 = Var(initialize=15, bounds=(1,30)) # kJ/mol
# Concentration variables indexed by time
m.CA = Var(m.t, initialize = m.CA0) # mol/L
m.CB = Var(m.t, initialize = 0) # mol/L
m.CC = Var(m.t, initialize = 0) # mol/L
# Derivatives in the model
#
m.dCA = DerivativeVar(m.CA)
m.dCB = DerivativeVar(m.CB)
m.dCC = DerivativeVar(m.CC)
# kinetic rate constants from Arrhenius equation
m.k1 = Expression(rule = m.A1 * exp(-m.E1*1000/(m.R*m.T))) # 1/hr
m.k2 = Expression(rule = m.A2 * exp(-m.E2*1000/(m.R*m.T))) # 1/hr
# Constraints to change concentrations based on kinetics
def conc_A(m,i):
return m.dCA[i] == - m.k1 * m.CA[i]
m.CA_rate = Constraint(m.t,rule=conc_A)
def conc_B(m,i):
return m.dCB[i] == m.k1 * m.CA[i] - m.k2 * m.CB[i]
m.CB_rate = Constraint(m.t,rule=conc_B)
def conc_C(m,i):
return m.dCC[i] == m.k2 * m.CB[i]
m.CC_rate = Constraint(m.t,rule=conc_C)
# Initial Conditions
def _initcon(m):
yield m.CA[m.t.first()] == m.CA0
yield m.CB[m.t.first()] == 0.0
yield m.CC[m.t.first()] == 0.0
m.initcon = ConstraintList(rule=_initcon)
# Objective function
# The objective function for parmest is defined as a 2-stage stochastic optimization objective function
# First stage cost: independent of scenarios ('experiments')
# expression for minimizing fixed realization
# from model. Eg.: reactor temperature, size, etc.
def ComputeFirstStageCost_rule(m):
# In this case, we do not optimize anything besides the kinetic parameters through
# least square fitting realizations at each timestep defined by m.t.
# Hence, the first stage cost is set to 0 here.
return 0
m.FirstStageCost = Expression(rule=ComputeFirstStageCost_rule)
# Second stage cost: Realization at each scenario over which the model is defined
def ComputeSecondStageCost_rule(m):
# In this problem, we want to minimize the sum of squared errors between
# 'experimental' data and the model realization of concentrations of
# A, B, and C over each scenario (here, timesteps defined by m.t)
return sum((m.CA[t] - m.Ca_meas[t]) ** 2 + (m.CB[t] - m.Cb_meas[t]) ** 2
+ (m.CC[t] - m.Cc_meas[t]) ** 2 for t in m.t)
m.SecondStageCost = Expression(rule=ComputeSecondStageCost_rule)
# return the sum of the first-stage and second-stage costs as the objective function
def total_cost_rule(m):
return m.FirstStageCost + m.SecondStageCost
m.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)
# Discretize model
disc = TransformationFactory('dae.collocation')
# The DAE model is discretized using a collocation scheme with 20 finite elements
# and 4 collocation points per finite element.
# Increasing the number of finite elements will improve the solution quality but,
# it will also increase the solution time. In our example, oweing to the small problem size,
# increasing the number of finite elements will not impact the solution time drastically.
disc.apply_to(m, nfe=20, ncp=4)
return m
5.1.6.1. Parameter estimation with parmest#
In the following cell, we perform parameter estimation using parmest to solve the least squares problem defined in the Pyomo dynamic model.
# import parmest
import pyomo.contrib.parmest.parmest as parmest
# defining the names of the parameters in a list
theta_names = ['A1','A2','E1','E2']
# create an object using parmest.Estimator() that stores the Pyomo model realizations for the datasets provided.
# This object which will be used to determined the parameter values that best fit all the datasets
pest = parmest.Estimator(create_model_DAE,data_dict_overall,theta_names,tee=True)
# call the method theta_est() for the Estimator() object defined above to solve
# the parameter estimation problem.
# theta_est() returns:
# the overall objective function value
# estimated parameter values (dictionary with keys = parameters names as defined in the Pyomo model)
obj, theta = pest.theta_est()
print('theta:\n',theta)
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.7039989e+01 1.98e+01 2.68e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.3426370e-01 2.11e+00 3.16e-01 -1.0 1.74e+01 - 9.58e-01 1.00e+00h 1
2 2.5138069e-01 2.48e-02 3.29e-02 -1.0 8.43e-01 - 1.00e+00 1.00e+00h 1
3 2.2289282e-01 7.55e-02 4.31e-03 -2.5 1.61e+00 - 9.95e-01 1.00e+00h 1
4 2.2225974e-01 6.13e-04 3.19e-06 -2.5 1.79e+00 - 1.00e+00 1.00e+00h 1
5 2.2221660e-01 9.42e-04 5.30e-07 -3.8 2.02e+00 - 1.00e+00 1.00e+00h 1
6 2.2211318e-01 1.91e-02 8.81e-06 -5.7 9.02e+00 - 9.36e-01 1.00e+00h 1
7 2.2210769e-01 7.24e-04 5.48e-07 -5.7 1.66e+00 - 1.00e+00 1.00e+00h 1
8 2.2210762e-01 2.28e-07 4.65e-10 -5.7 1.24e-01 - 1.00e+00 1.00e+00h 1
9 2.2210762e-01 5.67e-07 2.65e-10 -8.6 4.74e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.2210762e-01 1.87e-11 3.84e-14 -8.6 1.13e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 10
(scaled) (unscaled)
Objective...............: 2.2210762193125719e-01 2.2210762193125719e-01
Dual infeasibility......: 3.8434198309153839e-14 3.8434198309153839e-14
Constraint violation....: 1.8748558261449944e-11 1.8748558261449944e-11
Complementarity.........: 2.5059125114858520e-09 2.5059125114858520e-09
Overall NLP error.......: 2.5059125114858520e-09 2.5059125114858520e-09
Number of objective function evaluations = 11
Number of objective gradient evaluations = 11
Number of equality constraint evaluations = 11
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 10
Total CPU secs in IPOPT (w/o function evaluations) = 0.086
Total CPU secs in NLP function evaluations = 0.026
EXIT: Optimal Solution Found.
theta:
{'A1': 185.60880919391812, 'A2': 401.170198669058, 'E1': 9.866878980549787, 'E2': 14.866030768895396}
5.1.6.2. Plotting fitted model simulation with ‘experimental’ data#
Next, we plot the ‘experimental’ data along with the profiles generated using the fitted kinetic model. The symbols represent the ‘experimental’ data and the solid and dashed lines are the profiles generated using the fitted model.
# Parameter values from parameter estimation using parmest
A1 = theta['A1']
E1 = theta['E1']
A2 = theta['A2']
E2 = theta['E2']
A_est1 = [A1, A2]
A_est = np.asarray(A_est1)
E_est1 = [E1, E2]
E_est = np.asarray(E_est1)
ctr = 0
for T in T_vals:
for CA0 in CA0_vals:
# generate concentration profiles using estimated parameter values
k = kinetics(A_est, E_est, T)
# plot model-generated and 'experimental' data
# symbols for 'experimental' data
# solid and dashed lines for model-generated data
plot_exp(k, CA0, data_dict_overall[ctr]['data'], 'Model prediction and experimental value at T = {} K and $C_{}$ = {} mol/L'.format(T,'A0',CA0))
ctr+=1
5.1.7. Local uncertainty analysis#
5.1.7.1. Covariance matrix#
The parameter covariance matrix is calculated using the reduced Hessian approach. Using parmest, the covariance matrix can be calculated by setting optional argument calc_cov to True. More information on this approach can be found here: https://doi.org/10.1002/aic.16242
# defining the names of the parameters in a list
theta_names = ['A1','A2','E1','E2']
# create an object using parmest.Estimator() that stores the Pyomo model realizations for the datasets provided.
# This object which will be used to determined the parameter values that best fit all the datasets
pest = parmest.Estimator(create_model_DAE,data_dict_overall,theta_names,tee=True)
# call the method theta_est() for the Estimator() object defined above to solve
# the parameter estimation problem.
# Argument calc_cov=True is used to return the reduced Hessian matrix for parameter sensitivities
obj, theta, cov = pest.theta_est(calc_cov=True)
print('theta:\n',theta)
Ipopt 3.13.2: bound_relax_factor=0
honor_original_bounds=no
******************************************************************************
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.7039989e+01 1.98e+01 2.68e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.3426370e-01 2.11e+00 3.16e-01 -1.0 1.74e+01 - 9.58e-01 1.00e+00h 1
2 2.5138069e-01 2.48e-02 3.29e-02 -1.0 8.43e-01 - 1.00e+00 1.00e+00h 1
3 2.2289282e-01 7.55e-02 4.31e-03 -2.5 1.61e+00 - 9.95e-01 1.00e+00h 1
4 2.2225974e-01 6.13e-04 3.19e-06 -2.5 1.79e+00 - 1.00e+00 1.00e+00h 1
5 2.2221660e-01 9.42e-04 5.30e-07 -3.8 2.02e+00 - 1.00e+00 1.00e+00h 1
6 2.2211318e-01 1.91e-02 8.81e-06 -5.7 9.02e+00 - 9.36e-01 1.00e+00h 1
7 2.2210769e-01 7.24e-04 5.48e-07 -5.7 1.66e+00 - 1.00e+00 1.00e+00h 1
8 2.2210762e-01 2.28e-07 4.65e-10 -5.7 1.24e-01 - 1.00e+00 1.00e+00h 1
9 2.2210762e-01 5.67e-07 2.65e-10 -8.6 4.74e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.2210762e-01 1.87e-11 3.84e-14 -8.6 1.13e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 10
(scaled) (unscaled)
Objective...............: 2.2210762193125719e-01 2.2210762193125719e-01
Dual infeasibility......: 3.8437985232109095e-14 3.8437985232109095e-14
Constraint violation....: 1.8748558261449944e-11 1.8748558261449944e-11
Complementarity.........: 2.5059125114985724e-09 2.5059125114985724e-09
Overall NLP error.......: 2.5059125114985724e-09 2.5059125114985724e-09
Number of objective function evaluations = 11
Number of objective gradient evaluations = 11
Number of equality constraint evaluations = 11
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 10
Total CPU secs in IPOPT (w/o function evaluations) = 0.087
Total CPU secs in NLP function evaluations = 0.026
EXIT: Optimal Solution Found.
theta:
{'A1': 185.60880919391855, 'A2': 401.1701986690201, 'E1': 9.866878980549789, 'E2': 14.866030768895133}
5.1.7.2. Parameter identifiability#
The Fisher information matrix, \(FIM\), is calculated as the inverse of the parameter covariance matrix: \(FIM = {(Cov.)}^{-1}\)
Eigen decomposition of \(FIM\) following \(FIM = v \lambda {v}^{-1}\) gives eigenvector matrix \(v\) and \(\lambda\) is the diagonal matrix of eigenvalues. Column \(i\) in \(v\) is the eigenvector corresponding to the \(i\)-th eigenvalue. Further, each element of the eigenvector corresponds to the fitted parameters in order in which they appear in in theta_names. The magnitude of each element of an eigenvector denotes the contribution of the parameters on the direciton of the unit vector.
The eigenvector corresponding to the smallest eigenvalue denotes the direction of least variance in the parameter space. The parameter that corresponds to the major contributor in the eigenvector, then, has the lowest impact on model fit quality. Thus, this parameter is considered sloppy and fixing it’s value will not affect overall model behavior.
# Fisher information matrix can be computed using the inverse of the reduced Hessian
fim = np.linalg.inv(cov)
# Eigen decomposition of the Fisher information matrix
eig_values, eig_vectors = np.linalg.eig(fim)
for i,eig in enumerate(eig_values):
print('***************************************************************')
print('\nEigen value: {:0.3e}\n'.format(eig))
print('=== Eigen vector elements with correspondng parameter names ===\n')
print('------------------------------')
print('| Vector element | Parameter |')
print('------------------------------')
for j,theta_name in enumerate(theta_names):
if eig_vectors[i,j] < 0.0:
print('| {:0.3e} | {} |'.format(eig_vectors[i,j],theta_name))
else:
print('| {:0.3e} | {} |'.format(eig_vectors[i,j],theta_name))
print('\n')
***************************************************************
Eigen value: 1.559e+01
=== Eigen vector elements with correspondng parameter names ===
------------------------------
| Vector element | Parameter |
------------------------------
| 1.254e-02 | A1 |
| -1.789e-03 | A2 |
| 9.985e-01 | E1 |
| -5.373e-02 | E2 |
***************************************************************
Eigen value: 8.808e+00
=== Eigen vector elements with correspondng parameter names ===
------------------------------
| Vector element | Parameter |
------------------------------
| 1.519e-03 | A1 |
| 6.805e-03 | A2 |
| 5.373e-02 | E1 |
| 9.985e-01 | E2 |
***************************************************************
Eigen value: 5.745e-05
=== Eigen vector elements with correspondng parameter names ===
------------------------------
| Vector element | Parameter |
------------------------------
| -9.794e-01 | A1 |
| 2.017e-01 | A2 |
| 1.263e-02 | E1 |
| -5.644e-04 | E2 |
***************************************************************
Eigen value: 6.451e-06
=== Eigen vector elements with correspondng parameter names ===
------------------------------
| Vector element | Parameter |
------------------------------
| -2.017e-01 | A1 |
| -9.794e-01 | A2 |
| 1.150e-03 | E1 |
| 6.919e-03 | E2 |
Looking at the eigenvalues and corresponding eigenvectors, we can see that parameters A2 and A1 (in order) have the largest contributions to the direction of the eigenvector corresponding to the smallest eigenvalue. Therefore, as discussed above, we can conclude from this analysis that parameter A2 is the least identifiable parameter in this kinetic model, followed by parameter A1.
5.1.8. Bootstrap resampling#
Bootstrapping is a resampling method by independently sampling with replacement from an existing sample data with same sample size n, and performing inference among these resampled data (link). Bootstrap resampling is often used in parameter estimation problems to determine parameter confidence intervals. More information about bootstrap resampling and confidence interval calculation can be found here.
theta_est_bootstrap() is used to perform resampling with parmest. More information can be found here.
parmest also provides functions to plot bootstrap parameter estimates along with various confidence intervals link.
# create Estimator object
pest = parmest.Estimator(create_model_DAE,data_dict_overall,theta_names,tee=True)
### Parameter estimation with bootstrap resampling
#
bootstrap_theta = pest.theta_est_bootstrap(10)
print(bootstrap_theta.head())
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.5867818e+01 1.98e+01 2.55e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.1317553e-01 2.14e+00 3.19e-01 -1.0 1.74e+01 - 9.57e-01 1.00e+00h 1
2 2.6873829e-01 1.33e-02 3.37e-02 -1.0 8.47e-01 - 1.00e+00 1.00e+00h 1
3 2.4267977e-01 7.48e-02 4.37e-03 -2.5 1.60e+00 - 9.95e-01 1.00e+00h 1
4 2.4236905e-01 7.27e-04 2.31e-06 -2.5 1.82e+00 - 1.00e+00 1.00e+00h 1
5 2.4232473e-01 9.37e-04 4.09e-07 -3.8 2.04e+00 - 1.00e+00 1.00e+00h 1
6 2.4221255e-01 1.64e-02 8.31e-06 -5.7 8.44e+00 - 9.41e-01 1.00e+00h 1
7 2.4219362e-01 1.56e-03 3.30e-06 -5.7 1.02e+01 - 1.00e+00 1.00e+00h 1
8 2.4219368e-01 1.52e-05 3.30e-08 -5.7 9.84e-01 - 1.00e+00 1.00e+00h 1
9 2.4219366e-01 3.52e-06 7.46e-09 -8.6 4.73e-01 - 9.97e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.4219366e-01 2.80e-09 5.91e-12 -8.6 1.33e-02 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 10
(scaled) (unscaled)
Objective...............: 2.4219366015053306e-01 2.4219366015053306e-01
Dual infeasibility......: 5.9064320719730795e-12 5.9064320719730795e-12
Constraint violation....: 2.8011406527639338e-09 2.8011406527639338e-09
Complementarity.........: 2.5083856818837811e-09 2.5083856818837811e-09
Overall NLP error.......: 2.8011406527639338e-09 2.8011406527639338e-09
Number of objective function evaluations = 11
Number of objective gradient evaluations = 11
Number of equality constraint evaluations = 11
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 11
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 10
Total CPU secs in IPOPT (w/o function evaluations) = 0.094
Total CPU secs in NLP function evaluations = 0.027
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.7498853e+01 1.98e+01 2.54e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.5707852e-01 2.24e+00 3.32e-01 -1.0 1.72e+01 - 9.55e-01 1.00e+00h 1
2 2.6441294e-01 1.23e-02 3.69e-02 -1.0 8.86e-01 - 1.00e+00 1.00e+00h 1
3 2.2481185e-01 1.22e-01 3.95e-03 -2.5 2.03e+00 - 9.86e-01 1.00e+00h 1
4 2.2296649e-01 9.62e-04 9.74e-06 -2.5 2.30e+00 - 1.00e+00 1.00e+00h 1
5 2.2285645e-01 2.33e-03 9.46e-07 -3.8 3.14e+00 - 1.00e+00 1.00e+00h 1
6 2.2262528e-01 3.63e-02 1.48e-05 -3.8 1.22e+01 - 1.00e+00 1.00e+00h 1
7 2.2262252e-01 1.41e-04 1.34e-07 -3.8 8.75e-01 - 1.00e+00 1.00e+00h 1
8 2.2260420e-01 2.71e-03 1.17e-06 -5.7 3.19e+00 - 9.85e-01 1.00e+00h 1
9 2.2260208e-01 8.82e-05 1.96e-07 -5.7 2.44e+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 2.2260206e-01 2.22e-08 7.42e-11 -5.7 3.84e-02 - 1.00e+00 1.00e+00h 1
11 2.2260206e-01 7.04e-07 5.64e-10 -8.6 1.31e-01 - 1.00e+00 1.00e+00h 1
12 2.2260206e-01 9.65e-11 2.12e-13 -8.6 2.53e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 12
(scaled) (unscaled)
Objective...............: 2.2260205708195835e-01 2.2260205708195835e-01
Dual infeasibility......: 2.1180235186969325e-13 2.1180235186969325e-13
Constraint violation....: 9.6474828126247303e-11 9.6474828126247303e-11
Complementarity.........: 2.5059753197903154e-09 2.5059753197903154e-09
Overall NLP error.......: 2.5059753197903154e-09 2.5059753197903154e-09
Number of objective function evaluations = 13
Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 12
Total CPU secs in IPOPT (w/o function evaluations) = 0.115
Total CPU secs in NLP function evaluations = 0.035
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.3160828e+01 1.98e+01 2.44e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.8848316e-01 2.15e+00 3.43e-01 -1.0 1.73e+01 - 9.57e-01 1.00e+00h 1
2 2.1823275e-01 2.59e-02 3.26e-02 -1.0 8.32e-01 - 1.00e+00 1.00e+00h 1
3 1.9588731e-01 9.56e-02 3.26e-03 -2.5 1.79e+00 - 9.89e-01 1.00e+00h 1
4 1.9426722e-01 5.67e-03 2.98e-05 -2.5 4.99e+00 - 1.00e+00 1.00e+00h 1
5 1.9381613e-01 9.82e-03 5.75e-06 -3.8 6.41e+00 - 9.93e-01 1.00e+00h 1
6 1.9250198e-01 2.13e-01 1.33e-04 -3.8 2.84e+01 - 1.00e+00 1.00e+00h 1
7 1.9247801e-01 1.81e-03 2.17e-06 -3.8 3.35e+00 - 1.00e+00 1.00e+00h 1
8 1.9248177e-01 1.39e-05 1.14e-08 -3.8 2.16e-01 - 1.00e+00 1.00e+00h 1
9 1.9226978e-01 1.56e-02 1.12e-05 -5.7 1.38e+01 - 9.08e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.9215578e-01 8.01e-03 1.73e-05 -5.7 2.14e+01 - 1.00e+00 1.00e+00h 1
11 1.9214378e-01 6.14e-04 9.61e-07 -5.7 5.66e+00 - 1.00e+00 1.00e+00h 1
12 1.9214299e-01 1.37e-05 2.67e-08 -5.7 8.36e-01 - 1.00e+00 1.00e+00h 1
13 1.9214229e-01 8.20e-05 1.90e-07 -8.6 2.04e+00 - 9.89e-01 1.00e+00h 1
14 1.9214228e-01 4.19e-07 8.58e-10 -8.6 1.45e-01 - 1.00e+00 1.00e+00h 1
15 1.9214228e-01 9.01e-12 1.82e-14 -8.6 6.71e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 15
(scaled) (unscaled)
Objective...............: 1.9214227916791748e-01 1.9214227916791748e-01
Dual infeasibility......: 1.8210932698128006e-14 1.8210932698128006e-14
Constraint violation....: 9.0079055325986701e-12 9.0079055325986701e-12
Complementarity.........: 2.5059130479412850e-09 2.5059130479412850e-09
Overall NLP error.......: 2.5059130479412850e-09 2.5059130479412850e-09
Number of objective function evaluations = 16
Number of objective gradient evaluations = 16
Number of equality constraint evaluations = 16
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 16
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 15
Total CPU secs in IPOPT (w/o function evaluations) = 0.120
Total CPU secs in NLP function evaluations = 0.039
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.6843041e+01 1.48e+01 2.00e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.2887378e-01 1.64e+00 3.10e-01 -1.0 1.30e+01 - 9.57e-01 1.00e+00h 1
2 2.3856072e-01 1.61e-02 3.12e-02 -1.0 6.34e-01 - 1.00e+00 1.00e+00h 1
3 2.0569218e-01 8.58e-02 1.55e-03 -2.5 1.39e+00 - 9.88e-01 1.00e+00h 1
4 2.0397947e-01 2.46e-03 1.97e-05 -2.5 4.27e+00 - 1.00e+00 1.00e+00h 1
5 2.0362946e-01 6.28e-03 3.22e-06 -3.8 5.86e+00 - 9.98e-01 1.00e+00h 1
6 2.0238769e-01 2.52e-01 1.54e-04 -3.8 3.50e+01 - 1.00e+00 1.00e+00h 1
7 2.0239576e-01 3.87e-03 5.28e-06 -3.8 5.15e+00 - 1.00e+00 1.00e+00h 1
8 2.0240503e-01 4.29e-05 3.26e-08 -3.8 4.24e-01 - 1.00e+00 1.00e+00h 1
9 2.0225622e-01 2.69e-02 1.77e-05 -5.7 1.01e+01 - 9.43e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.0223592e-01 3.44e-03 1.90e-06 -5.7 4.49e+00 - 1.00e+00 1.00e+00h 1
11 2.0223489e-01 1.32e-05 1.97e-08 -5.7 8.19e-01 - 1.00e+00 1.00e+00h 1
12 2.0223483e-01 1.82e-05 1.23e-08 -8.6 3.34e-01 - 9.98e-01 1.00e+00h 1
13 2.0223482e-01 1.56e-09 4.26e-12 -8.6 1.22e-02 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 2.0223482293125566e-01 2.0223482293125566e-01
Dual infeasibility......: 4.2614681659941213e-12 4.2614681659941213e-12
Constraint violation....: 1.5566685718226836e-09 1.5566685718226836e-09
Complementarity.........: 2.5071808827341051e-09 2.5071808827341051e-09
Overall NLP error.......: 2.5071808827341051e-09 2.5071808827341051e-09
Number of objective function evaluations = 14
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 14
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 13
Total CPU secs in IPOPT (w/o function evaluations) = 0.105
Total CPU secs in NLP function evaluations = 0.033
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.9317227e+01 1.98e+01 2.60e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 6.1238123e-01 2.24e+00 3.18e-01 -1.0 1.72e+01 - 9.55e-01 1.00e+00h 1
2 2.7818151e-01 2.30e-02 3.74e-02 -1.0 8.98e-01 - 1.00e+00 1.00e+00h 1
3 2.3911692e-01 9.48e-02 5.53e-03 -2.5 1.79e+00 - 9.91e-01 1.00e+00h 1
4 2.3810327e-01 5.27e-04 6.57e-06 -2.5 1.65e+00 - 1.00e+00 1.00e+00h 1
5 2.3805062e-01 8.90e-04 4.20e-07 -3.8 1.97e+00 - 1.00e+00 1.00e+00h 1
6 2.3788233e-01 1.11e-02 5.47e-06 -5.7 1.20e+01 - 9.07e-01 1.00e+00h 1
7 2.3783536e-01 2.86e-03 7.43e-06 -5.7 1.35e+01 - 1.00e+00 1.00e+00h 1
8 2.3783398e-01 3.89e-05 1.03e-07 -5.7 1.53e+00 - 1.00e+00 1.00e+00h 1
9 2.3783397e-01 2.27e-08 5.63e-11 -5.7 3.68e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.3783393e-01 3.86e-06 1.01e-08 -8.6 4.82e-01 - 9.98e-01 1.00e+00h 1
11 2.3783392e-01 1.28e-09 3.27e-12 -8.6 8.73e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 2.3783392446349957e-01 2.3783392446349957e-01
Dual infeasibility......: 3.2731232752983274e-12 3.2731232752983274e-12
Constraint violation....: 1.2825678297190279e-09 1.2825678297190279e-09
Complementarity.........: 2.5079696400494811e-09 2.5079696400494811e-09
Overall NLP error.......: 2.5079696400494811e-09 2.5079696400494811e-09
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.093
Total CPU secs in NLP function evaluations = 0.029
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.2902314e+01 1.98e+01 2.82e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.8525024e-01 1.87e+00 3.20e-01 -1.0 1.77e+01 - 9.63e-01 1.00e+00h 1
2 2.0748122e-01 4.78e-02 2.63e-02 -1.0 7.78e-01 - 1.00e+00 1.00e+00h 1
3 1.8933165e-01 5.96e-02 1.67e-03 -2.5 1.46e+00 - 1.00e+00 1.00e+00h 1
4 1.8872332e-01 3.36e-04 6.26e-06 -2.5 1.34e+00 - 1.00e+00 1.00e+00h 1
5 1.8869519e-01 5.93e-04 2.46e-07 -3.8 1.58e+00 - 1.00e+00 1.00e+00h 1
6 1.8858365e-01 2.95e-02 1.00e-05 -5.7 1.11e+01 - 9.17e-01 1.00e+00h 1
7 1.8856721e-01 4.89e-03 1.49e-06 -5.7 5.78e+00 - 1.00e+00 1.00e+00h 1
8 1.8856667e-01 3.03e-05 3.40e-08 -5.7 1.47e+00 - 1.00e+00 1.00e+00h 1
9 1.8856667e-01 2.58e-09 5.67e-12 -5.7 1.41e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.8856665e-01 4.34e-06 1.81e-09 -8.6 3.65e-01 - 9.98e-01 1.00e+00h 1
11 1.8856665e-01 2.68e-09 2.73e-12 -8.6 1.39e-02 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 1.8856664975547494e-01 1.8856664975547494e-01
Dual infeasibility......: 2.7346159208273144e-12 2.7346159208273144e-12
Constraint violation....: 2.6839197531103309e-09 2.6839197531103309e-09
Complementarity.........: 2.5075979185707942e-09 2.5075979185707942e-09
Overall NLP error.......: 2.6839197531103309e-09 2.6839197531103309e-09
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.093
Total CPU secs in NLP function evaluations = 0.028
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.7354755e+01 1.98e+01 2.98e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.2013292e-01 1.86e+00 2.83e-01 -1.0 1.77e+01 - 9.64e-01 1.00e+00h 1
2 2.7202560e-01 6.97e-02 2.66e-02 -1.0 8.35e-01 - 1.00e+00 1.00e+00h 1
3 2.5654573e-01 2.56e-02 1.19e-03 -2.5 9.75e-01 - 1.00e+00 1.00e+00h 1
4 2.5651613e-01 9.54e-05 1.45e-07 -3.8 7.33e-01 - 1.00e+00 1.00e+00h 1
5 2.5650277e-01 2.08e-03 9.86e-07 -5.7 3.12e+00 - 9.83e-01 1.00e+00h 1
6 2.5650032e-01 3.30e-04 3.07e-07 -5.7 3.02e+00 - 1.00e+00 1.00e+00h 1
7 2.5650025e-01 1.24e-07 1.53e-10 -5.7 6.35e-02 - 1.00e+00 1.00e+00h 1
8 2.5650025e-01 1.95e-07 4.47e-10 -8.6 1.13e-01 - 1.00e+00 1.00e+00h 1
9 2.5650025e-01 9.76e-11 2.27e-13 -8.6 2.55e-03 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 2.5650025180692537e-01 2.5650025180692537e-01
Dual infeasibility......: 2.2664407575036698e-13 2.2664407575036698e-13
Constraint violation....: 9.7613472860302863e-11 9.7613472860302863e-11
Complementarity.........: 2.5059619932443772e-09 2.5059619932443772e-09
Overall NLP error.......: 2.5059619932443772e-09 2.5059619932443772e-09
Number of objective function evaluations = 10
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 10
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 10
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 9
Total CPU secs in IPOPT (w/o function evaluations) = 0.079
Total CPU secs in NLP function evaluations = 0.023
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.5851784e+01 1.98e+01 2.49e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 4.9556238e-01 2.25e+00 3.44e-01 -1.0 1.72e+01 - 9.55e-01 1.00e+00h 1
2 2.4316470e-01 3.57e-02 3.54e-02 -1.0 8.67e-01 - 1.00e+00 1.00e+00h 1
3 2.1006658e-01 1.17e-01 7.41e-03 -2.5 1.97e+00 - 9.85e-01 1.00e+00h 1
4 2.0912807e-01 1.42e-03 3.67e-05 -2.5 1.17e+00 - 1.00e+00 1.00e+00h 1
5 2.0908236e-01 7.60e-04 4.78e-07 -3.8 1.80e+00 - 1.00e+00 1.00e+00h 1
6 2.0873382e-01 1.46e-02 1.46e-05 -5.7 2.17e+01 - 8.37e-01 1.00e+00h 1
7 2.0852924e-01 1.40e-02 2.84e-05 -5.7 2.85e+01 - 1.00e+00 1.00e+00h 1
8 2.0850600e-01 1.42e-03 2.78e-06 -5.7 8.51e+00 - 1.00e+00 1.00e+00h 1
9 2.0850448e-01 6.41e-05 1.28e-07 -5.7 1.78e+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 2.0850444e-01 1.31e-07 2.63e-10 -5.7 8.01e-02 - 1.00e+00 1.00e+00h 1
11 2.0850378e-01 7.84e-05 1.60e-07 -8.6 1.96e+00 - 9.89e-01 1.00e+00h 1
12 2.0850378e-01 3.53e-07 7.10e-10 -8.6 1.31e-01 - 1.00e+00 1.00e+00h 1
13 2.0850378e-01 6.04e-12 1.21e-14 -8.6 5.41e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 2.0850378243005299e-01 2.0850378243005299e-01
Dual infeasibility......: 1.2149736707217614e-14 1.2149736707217614e-14
Constraint violation....: 6.0387250755411515e-12 6.0387250755411515e-12
Complementarity.........: 2.5059099799274794e-09 2.5059099799274794e-09
Overall NLP error.......: 2.5059099799274794e-09 2.5059099799274794e-09
Number of objective function evaluations = 14
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 14
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 13
Total CPU secs in IPOPT (w/o function evaluations) = 0.105
Total CPU secs in NLP function evaluations = 0.034
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 2.3553727e+01 1.98e+01 2.74e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 7.5048522e-01 2.26e+00 2.77e-01 -1.0 1.72e+01 - 9.55e-01 1.00e+00h 1
2 2.5887790e-01 2.47e-02 3.76e-02 -1.0 9.01e-01 - 1.00e+00 1.00e+00h 1
3 2.0714094e-01 1.07e-01 3.43e-03 -2.5 1.90e+00 - 9.89e-01 1.00e+00h 1
4 2.0543545e-01 5.18e-04 6.73e-06 -2.5 1.11e+00 - 1.00e+00 1.00e+00h 1
5 2.0539974e-01 3.70e-04 2.10e-07 -3.8 1.30e+00 - 1.00e+00 1.00e+00h 1
6 2.0519893e-01 8.80e-03 3.97e-06 -5.7 1.54e+01 - 8.82e-01 1.00e+00h 1
7 2.0509916e-01 4.88e-03 5.65e-06 -5.7 1.93e+01 - 1.00e+00 1.00e+00h 1
8 2.0508880e-01 6.00e-04 8.17e-07 -5.7 7.05e+00 - 1.00e+00 1.00e+00h 1
9 2.0508828e-01 9.19e-06 1.25e-08 -5.7 8.82e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.0508805e-01 1.96e-05 2.04e-08 -8.6 1.28e+00 - 9.93e-01 1.00e+00h 1
11 2.0508805e-01 4.71e-08 5.57e-11 -8.6 6.31e-02 - 1.00e+00 1.00e+00h 1
12 2.0508805e-01 1.46e-13 1.31e-16 -8.6 9.03e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 12
(scaled) (unscaled)
Objective...............: 2.0508804642313905e-01 2.0508804642313905e-01
Dual infeasibility......: 1.3064608659185325e-16 1.3064608659185325e-16
Constraint violation....: 9.5923269327613525e-14 1.4566126083082054e-13
Complementarity.........: 2.5059036690133731e-09 2.5059036690133731e-09
Overall NLP error.......: 2.5059036690133731e-09 2.5059036690133731e-09
Number of objective function evaluations = 13
Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 12
Total CPU secs in IPOPT (w/o function evaluations) = 0.100
Total CPU secs in NLP function evaluations = 0.032
EXIT: Optimal Solution Found.
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...: 42648
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5680
Total number of variables............................: 7840
variables with only lower bounds: 0
variables with lower and upper bounds: 64
variables with only upper bounds: 0
Total number of equality constraints.................: 7836
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 1.1977527e+01 1.48e+01 2.02e-02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.8335150e-01 1.43e+00 3.16e-01 -1.0 1.32e+01 - 9.63e-01 1.00e+00h 1
2 2.2909936e-01 2.04e-02 2.50e-02 -1.0 5.67e-01 - 1.00e+00 1.00e+00h 1
3 2.1130723e-01 6.06e-02 2.18e-03 -2.5 1.19e+00 - 9.96e-01 1.00e+00h 1
4 2.1022439e-01 3.08e-03 1.99e-05 -2.5 4.19e+00 - 1.00e+00 1.00e+00h 1
5 2.0999394e-01 4.03e-03 2.25e-06 -3.8 4.70e+00 - 1.00e+00 1.00e+00h 1
6 2.0918544e-01 1.49e-01 8.39e-05 -3.8 2.75e+01 - 1.00e+00 1.00e+00h 1
7 2.0919272e-01 1.78e-03 3.95e-06 -3.8 3.53e+00 - 1.00e+00 1.00e+00h 1
8 2.0919573e-01 7.40e-06 5.64e-09 -3.8 1.85e-01 - 1.00e+00 1.00e+00h 1
9 2.0908787e-01 1.75e-02 1.07e-05 -5.7 8.57e+00 - 9.50e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.0905499e-01 2.36e-03 5.62e-06 -5.7 1.51e+01 - 1.00e+00 1.00e+00h 1
11 2.0905247e-01 1.96e-04 4.92e-07 -5.7 4.55e+00 - 1.00e+00 1.00e+00h 1
12 2.0905240e-01 6.33e-07 1.64e-09 -5.7 2.61e-01 - 1.00e+00 1.00e+00h 1
13 2.0905217e-01 3.26e-05 7.85e-08 -8.6 1.85e+00 - 9.87e-01 1.00e+00h 1
14 2.0905217e-01 2.36e-07 5.69e-10 -8.6 1.59e-01 - 1.00e+00 1.00e+00h 1
15 2.0905217e-01 2.09e-12 5.22e-15 -8.6 4.75e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 15
(scaled) (unscaled)
Objective...............: 2.0905216767485774e-01 2.0905216767485774e-01
Dual infeasibility......: 5.2157307674572543e-15 5.2157307674572543e-15
Constraint violation....: 2.0867751970854442e-12 2.0867751970854442e-12
Complementarity.........: 2.5059046785673633e-09 2.5059046785673633e-09
Overall NLP error.......: 2.5059046785673633e-09 2.5059046785673633e-09
Number of objective function evaluations = 16
Number of objective gradient evaluations = 16
Number of equality constraint evaluations = 16
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 16
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 15
Total CPU secs in IPOPT (w/o function evaluations) = 0.119
Total CPU secs in NLP function evaluations = 0.039
EXIT: Optimal Solution Found.
A1 A2 E1 E2
0 186.769746 382.642388 9.907827 14.726285
1 179.703097 392.070899 9.738017 14.824106
2 156.529846 334.342272 9.464807 14.407121
3 146.617094 406.533938 9.252072 14.878677
4 189.635337 370.602660 9.907778 14.697935
Once the parameter estimates are generated through bootstrap resampling, we can visualize the estimates using the pairwise_plot() function as follows:
# plot parameter estimate
pyomo.contrib.parmest.graphics.pairwise_plot(bootstrap_theta, title='Bootstrap theta estimates')
Confidence regions can be plotted around the bootstrap estimates for various distributions with confidence \(\alpha\).
# plot bootstrap parameter estimates with confidence intervals
pyomo.contrib.parmest.graphics.pairwise_plot(bootstrap_theta, theta, 0.8, ['MVN', 'KDE', 'Rect'],
title='Bootstrap theta with confidence regions')
{'A1': 185.60880919391855, 'A2': 401.1701986690201, 'E1': 9.866878980549789, 'E2': 14.866030768895133}
/Users/kghosh/anaconda3/envs/idaes-pse-dev-pyomo6-1/lib/python3.8/site-packages/pyomo/contrib/parmest/graphics.py:151: UserWarning: No contour levels were found within the data range.
ax.contour(X,Y,Z, levels=[alpha], colors=color)
5.1.9. Nonlinear confidence regions#
The likelihood-ratio test (sometimes called the likelihood-ratio \(\chi^2\) test) is a hypothesis test that helps one choose the “best” model between two modelslink. Basically, the test compares the fit of two models. The null hypothesis is that the first model is the “best” model; It is rejected when the test statistic is large. In other words, if the null hypothesis is rejected, then the second model is a significant improvement over the first model.
In the last part of this notebook, we use the bootstrap parameter estimates to determine the goodness of fit using the likelihood-ratio \(\chi^2\) test. More information can be found here.
from itertools import product
### Likelihood ratio test
# generate arrays of parameter values
A1 = np.arange(180.0, 190.0, 1.0)
A2 = np.arange(395.0, 405.0, 1.0)
E1 = np.arange(5.0, 15.0, 1.0)
E2 = np.arange(10.0, 20.0, 1.0)
# format parameter values into a pandas dataframe to be provided as input to calculate
# corresponding objective function values
# theta_vals = pd.DataFrame(list(product(A1, A2, E1, E2)), columns=theta_names)
theta_vals = bootstrap_theta
obj_at_theta = pest.objective_at_theta(theta_vals)
print(obj_at_theta.head())
A1 A2 E1 E2 obj
0 186.769746 382.642388 9.907827 14.726285 0.222375
1 179.703097 392.070899 9.738017 14.824106 0.222957
2 156.529846 334.342272 9.464807 14.407121 0.224970
3 146.617094 406.533938 9.252072 14.878677 0.225126
4 189.635337 370.602660 9.907778 14.697935 0.222650
LR = pest.likelihood_ratio_test(obj_at_theta, obj, [0.8, 0.85, 0.9, 0.95])
print(LR.head())
A1 A2 E1 E2 obj 0.8 0.85 0.9 \
0 186.769746 382.642388 9.907827 14.726285 0.222375 True True True
1 179.703097 392.070899 9.738017 14.824106 0.222957 True True True
2 156.529846 334.342272 9.464807 14.407121 0.224970 True True True
3 146.617094 406.533938 9.252072 14.878677 0.225126 True True True
4 189.635337 370.602660 9.907778 14.697935 0.222650 True True True
0.95
0 True
1 True
2 True
3 True
4 True
The likelihood ratio test results with confidence \(\alpha\) can be visualized as follows:
pyomo.contrib.parmest.graphics.pairwise_plot(LR, theta, 0.8,
title='LR results within 80% confidence region')
Objective contour plot for A1 A2 slice failed
Objective contour plot for A1 E1 slice failed
Objective contour plot for A2 E1 slice failed
Objective contour plot for A1 E2 slice failed
Objective contour plot for A2 E2 slice failed
Objective contour plot for E1 E2 slice failed
Objective contour plot for A2 A1 slice failed
Objective contour plot for E1 A1 slice failed
Objective contour plot for E2 A1 slice failed
Objective contour plot for E1 A2 slice failed
Objective contour plot for E2 A2 slice failed
Objective contour plot for E2 E1 slice failed
