Uncertainty Quantification - Pyomo.DoE Workshop and Tutorials

Uncertainty quantification is required to ensure that the estimates of the parameters are close to their true values. This parameter accuracy is measured by computing the covariance matrix. The leading diagonal of this matrix contains the variance of the estimated parameters. The uncertainty in the parameters is taken as the standard deviation of the variance from the covariance matrix.

Assuming Gaussian (correlated or independent and identically distributed) measurement errors, the covariance matrix can be computed using three methods in ParmEst.

Reduced Hessian Method
Finite Difference Method
Automatic Differentiation Method

These methods are described in more details in the sections below.

Computing the Covariance Matrix with ParmEst¶

This section describes the steps in computing the covariance matrix in ParmEst.

Import the Necessary Packages¶

import sys
import pandas as pd
import numpy as np

# If running on Google Colab, install Pyomo and Ipopt via IDAES
on_colab = "google.colab" in sys.modules
if on_colab:
    !wget "https://raw.githubusercontent.com/dowlinglab/pyomo-doe/main/notebooks/tclab_pyomo.py"

# import TCLab model, simulation, and data analysis functions
from tclab_pyomo import (
    TC_Lab_data,
    TC_Lab_experiment,
    extract_results,
    extract_plot_results,
    extract_trace_covariance,
)

Load the Experimental Data¶

# load experimental data
if on_colab:
    file = "https://raw.githubusercontent.com/dowlinglab/pyomo-doe/main/data/tclab_sine_test_5min_period.csv"
else:
    file = "../data/tclab_sine_test_5min_period.csv"
df = pd.read_csv(file)

# store in custom data class
tc_data = TC_Lab_data(
    name="Sine Wave Test for Heater 1",
    time=df["Time"].values,
    T1=df["T1"].values,
    u1=df["Q1"].values,
    P1=200,
    TS1_data=None,
    T2=df["T2"].values,
    u2=df["Q2"].values,
    P2=200,
    TS2_data=None,
    Tamb=df["T1"].values[0],
)

# set default number of states in the TCLab model
number_tclab_states = 2

Get the Parameter Estimates¶

# import parmest
import pyomo.contrib.parmest.parmest as parmest

# define an Experiment object within parmest
TC_Lab_sine_exp = TC_Lab_experiment(data=tc_data, number_of_states=number_tclab_states)

# define an Estimator object using the weighted sum of squared errors objective function
pest = parmest.Estimator([TC_Lab_sine_exp], obj_function="SSE_weighted", tee=True)

# estimate the paramters
obj, theta = 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...:    15309
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     7203

Total number of variables............................:     3610
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     1804
                     variables with only upper bounds:        0
Total number of equality constraints.................:     3606
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  3.7689956e+05 3.12e-01 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  6.9950237e+04 1.34e-01 6.09e+03  -1.0 1.00e+01    -  4.72e-01 5.76e-01f  1
   2  4.4812180e+02 2.63e-03 1.60e+03  -1.0 4.33e+00    -  9.31e-01 1.00e+00f  1
   3  4.3029017e+02 1.98e-05 8.77e+01  -1.0 5.40e-02    -  1.00e+00 1.00e+00f  1
   4  4.3019471e+02 6.02e-03 1.14e+01  -1.0 6.03e-01    -  1.00e+00 1.00e+00h  1
   5  4.3019369e+02 5.29e-03 9.95e-01  -1.0 4.64e-01    -  1.00e+00 1.00e+00h  1
   6  4.3019125e+02 1.21e-04 2.63e-02  -1.7 4.73e-02    -  1.00e+00 1.00e+00h  1
   7  4.3019194e+02 7.05e-08 2.41e-04  -3.8 1.18e-03    -  1.00e+00 1.00e+00h  1
   8  4.3019194e+02 7.87e-11 7.30e-06  -5.7 3.38e-03    -  1.00e+00 1.00e+00H  1
   9  4.3019194e+02 4.01e-11 6.21e-08  -8.6 3.67e-04    -  1.00e+00 1.00e+00H  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  4.3019194e+02 7.12e-11 6.25e-08  -8.6 2.69e-01    -  1.00e+00 1.22e-04h 14
  11  4.3019194e+02 2.09e-03 3.34e-02  -8.6 2.69e-01    -  1.00e+00 1.00e+00h  1
  12  4.3019194e+02 2.09e-03 3.34e-02  -8.6 1.41e-02  -4.0 1.00e+00 3.05e-05h 16
  13  4.3019194e+02 2.42e-06 4.32e-05  -8.6 1.41e-02  -4.5 1.00e+00 1.00e+00s 22
  14  4.3019194e+02 3.05e-12 3.45e-04  -8.6 1.66e-05    -  1.00e+00 0.00e+00S 22
  15  4.3019194e+02 2.85e-14 1.81e-10  -8.6 1.83e-08  -5.0 1.00e+00 1.00e+00h  1

Number of Iterations....: 15

                                   (scaled)                 (unscaled)
Objective...............:   2.7451564684464381e+02    4.3019194276647835e+02
Dual infeasibility......:   1.8148423762684358e-10    2.8440293900758524e-10
Constraint violation....:   2.8463342793827451e-14    2.8463342793827451e-14
Complementarity.........:   2.5059035596800630e-09    3.9269875255390520e-09
Overall NLP error.......:   2.5059035596800630e-09    3.9269875255390520e-09


Number of objective function evaluations             = 69
Number of objective gradient evaluations             = 16
Number of equality constraint evaluations            = 69
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.166
Total CPU secs in NLP function evaluations           =      0.008

EXIT: Optimal Solution Found.

Covariance Matrix via Reduced Hessian¶

When the objective function is the sum of squared errors (SSE) for homogeneous data, defined as $\text{SSE} = \sum_{i = 1}^{n} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i}; \boldsymbol{\theta})\right)^\text{T} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i}; \boldsymbol{\theta})\right)$ , the covariance matrix can be computed as:

\boldsymbol{V}_{\boldsymbol{\theta}} = 2 \sigma^2 \left(\frac{\partial^2 \text{SSE}} {\partial \boldsymbol{\theta}^2}\right)^{-1}_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}}

(1)

Similarly, when the objective function is the weighted SSE (WSSE) for heterogeneous data, defined as $\text{WSSE} = \frac{1}{2} \sum_{i = 1}^{n} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta})\right)^\text{T} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta})\right)$ , the covariance matrix is:

\boldsymbol{V}_{\boldsymbol{\theta}} = \left(\frac{\partial^2 \text{WSSE}} {\partial \boldsymbol{\theta}^2}\right)^{-1}_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}}

(2)

Where $\boldsymbol{V}_{\boldsymbol{\theta}}$ is the covariance matrix of the estimated parameters $\hat{\boldsymbol{\theta}} \in \mathbb{R}^p$ , $\boldsymbol{y}_{i} \in \mathbb{R}^m$ are observations of the measured variables, $\boldsymbol{f}$ is the model function, $\boldsymbol{x}_{i} \in \mathbb{R}^{q}$ are the input variables, $n$ is the number of experiments, $\boldsymbol{\Sigma}_{\boldsymbol{y}}$ is the measurement error covariance matrix, and $\sigma^2$ is the variance of the measurement error.

This method computes the inverse of the Hessian by scaling the objective function (SSE or WSSE) with a constant probability factor, $\frac{1}{n}$ , and constructing the model optimality conditions.

# compute the covariance using the reduced hessian method in ParmEst
cov_red = pest.cov_est(method="reduced_hessian")

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...:    15309
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     7203

Total number of variables............................:     3610
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     1804
                     variables with only upper bounds:        0
Total number of equality constraints.................:     3606
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  4.3019194e+02 2.85e-14 1.24e-09  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  4.3019265e+02 1.97e-08 9.42e-02  -1.0 6.99e-04    -  9.91e-01 1.00e+00f  1
   2  4.3019197e+02 6.26e-10 1.02e-03  -1.7 4.62e-04    -  1.00e+00 1.00e+00h  1
   3  4.3019194e+02 2.94e-11 6.47e-05  -3.8 1.17e-04    -  1.00e+00 1.00e+00h  1
   4  4.3019194e+02 8.52e-14 3.63e-09  -5.7 1.67e-06    -  1.00e+00 1.00e+00h  1
   5  4.3019194e+02 7.11e-14 2.50e-10  -8.6 1.73e-06    -  1.00e+00 1.00e+00H  1

Number of Iterations....: 5

                                   (scaled)                 (unscaled)
Objective...............:   4.3019194276645783e+02    4.3019194276645783e+02
Dual infeasibility......:   2.5001809993250180e-10    2.5001809993250180e-10
Constraint violation....:   7.1109784727241276e-14    7.1109784727241276e-14
Complementarity.........:   2.5067148368443800e-09    2.5067148368443800e-09
Overall NLP error.......:   2.5067148368443800e-09    2.5067148368443800e-09


Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 6
Number of equality constraint evaluations            = 7
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 6
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 5
Total CPU secs in IPOPT (w/o function evaluations)   =      0.095
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

Covariance Matrix via Finite Difference¶

In this method, the covariance matrix, $\boldsymbol{V}_{\boldsymbol{\theta}}$ , is computed by differentiating the Hessian, $\frac{\partial^2 \text{SSE}}{\partial \boldsymbol{\theta}^2}$ or $\frac{\partial^2 \text{WSSE}}{\partial \boldsymbol{\theta}^2}$ , and applying Gauss-Newton approximation which results in:

\boldsymbol{V}_{\boldsymbol{\theta}} = \left(\sum_{i = 1}^n \boldsymbol{G}_{i}^{\text{T}} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} \boldsymbol{G}_{i} \right)^{-1}

(3)

where

\boldsymbol{G}_{i} = \frac{\partial \boldsymbol{f}(\boldsymbol{x}_i;\boldsymbol{\theta})} {\partial \boldsymbol{\theta}} \approx \frac{\boldsymbol{f}(\boldsymbol{x}_i;\theta_k + \Delta \theta_k) \vert_{\hat{\boldsymbol{\theta}}} - \boldsymbol{f}(\boldsymbol{x}_i;\theta_k - \Delta \theta_k) \vert_{\hat{\boldsymbol{\theta}}}}{2 \Delta \theta_k} \quad \forall \quad \theta_k \, \in \, [\theta_1,\cdots, \theta_p]

(4)

# compute the covariance using the finite difference method in ParmEst
cov_finite = pest.cov_est(method="finite_difference")

Fetching long content....

Covariance Matrix via Automatic Differentiation¶

Similar to the finite difference method, the covariance matrix can be computed as:

\boldsymbol{V}_{\boldsymbol{\theta}} = \left( \sum_{i = 1}^n \boldsymbol{G}_{\text{kaug},\, i}^{\text{T}} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} \boldsymbol{G}_{\text{kaug},\, i} \right)^{-1}

(5)

However, this method uses implicit differentiation and the optimality or Karush–Kuhn–Tucker (KKT) conditions of the model to compute the Jacobian matrix, $\boldsymbol{G}_{\text{kaug},\, i}$ , for experiment $i$ .

\boldsymbol{G}_{\text{kaug},\,i} = \frac{\partial \boldsymbol{f}(\boldsymbol{x}_i,\boldsymbol{\theta})} {\partial \boldsymbol{\theta}} + \frac{\partial \boldsymbol{f}(\boldsymbol{x}_i,\boldsymbol{\theta})} {\partial \boldsymbol{x}_i}\frac{\partial \boldsymbol{x}_i}{\partial \boldsymbol{\theta}}

(6)

# compute the covariance using the automatic differentiation method in ParmEst
cov_auto = pest.cov_est(method="automatic_differentiation_kaug")

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...:    10800
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:     3602
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     1800
                     variables with only upper bounds:        0
Total number of equality constraints.................:     3602
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  0.0000000e+00 3.95e-01 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 2.84e-14 9.12e-02  -1.0 1.01e+01    -  8.47e-01 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   2.8421709430404007e-14    2.8421709430404007e-14
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.8421709430404007e-14    2.8421709430404007e-14


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.007
Total CPU secs in NLP function evaluations           =      0.000

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...:    10800
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      900

Total number of variables............................:     3602
                     variables with only lower bounds:        0
                variables with lower and upper bounds:     1800
                     variables with only upper bounds:        0
Total number of equality constraints.................:     3602
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  3.7689956e+05 3.95e-01 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  4.3019194e+02 2.84e-14 9.12e-02  -1.0 1.01e+01    -  8.47e-01 1.00e+00f  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   2.7451564684464461e+02    4.3019194276647960e+02
Dual infeasibility......:   2.2231105845094135e-12    3.4838242276079054e-12
Constraint violation....:   2.8421709430404007e-14    2.8421709430404007e-14
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.8421709430404007e-14    3.4838242276079054e-12


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.007
Total CPU secs in NLP function evaluations           =      0.003

EXIT: Optimal Solution Found.
[K_AUG] 0.1.0, Part of the IDAES PSE framework
Please visit https://idaes.org/
: dsdp_mode=
dsdp_mode=
W[K_AUG]...	[K_AUG_ASL]No n_rhs declared.
W[K_AUG]...	[K_AUG_ASL]Target log10mu:= -11.
W[K_AUG]...	[K_AUG_ASL]This problem has no degrees of freedom
Pass the option square_override for the desired behaviour
W[K_AUG_ASL]...	[K_AUG_ASL]No ipopt_zL_in suffix declared, setting zL = 0.
W[K_AUG_ASL]...	[K_AUG_ASL]No ipopt_zU_in suffix declared, setting zU = 0.
W[K_AUG]...	[K_AUG_ASL]No f_timestamp suffix declared, Fallback to default writing mode.
I[K_AUG]...	[K_AUG_ASL] Filename for dot_sens dsdp_in_.in
W[K_AUG]...	[K_AUG_ASL]dsdp for linear C(x) + I*p = 0 override.
W[K_AUG]...	[MA57_DRIVER]

		Could not fix the accuracy of the problem.
		Try re-writing the problem or give a different point or change "max_refinement_steps"
		Warning: results might be incorrect.
		Current residual ratio 9.62313e-07; Max residual ratio 1e-10.

I[K_AUG]...	[K_AUG_ASL]File read successfully.
I[K_AUG]...	[K_AUG_ASL]Number of Right hand sides: 0
I[K_AUG]...	[K_AUG_ASL]Number of variables       : 3606
I[K_AUG]...	[K_AUG_ASL]Number of constraints     : 3606
I[K_AUG]...	[K_AUG_ASL]Number of valid n_dof     : 0
I[K_AUG]...	[ADJUST_MU]	Warning no relevant info from the problem can predict logmu
I[K_AUG]...	[FIND_INEQUALITIES]summary: eq: 3606, leq: 0, geq: 0 
I[K_AUG]...	[K_AUG_ASL]Nonzeroes in the sparse Jacobian 15305
I[K_AUG]...	[GET_HESS_ASL]Objective found
I[K_AUG]...	[GET_HESS_ASL]Nonzeroes in the sparse hessian 7203
I[K_AUG]...	[GET_HESS_ASL]Minimization problem detected
I[K_AUG]...	[GET_HESS_ASL]Current objective 430.191943
I[K_AUG]...	[GET_HESS_ASL]Missing nz in the Hessian of the Lag: 2706
I[K_AUG]...	[K_AUG_ASL]Barrier term added.
I[K_AUG]...	[K_AUG_ASL]MC19 scaling...
I[K_AUG]...	[ASSM_RHS_DCDP]According to the suffixes declared len p is 4 
I[MA57]...	[]***
I[MA57]...	[ma57bd_]ma57ad_: Status Ok.
I[MA57]...	[MA57_FACTOR]
I[K_AUG]...	[MA57_FACTOR]n_neig = 3606
I[K_AUG]...	[INERTIA_STRATEGY]Inertia check successful.
I[MA57]...	[Factorize Success.]***
I[MA57]...	[MA57_RATIO]
I[K_AUG]...	[MA57_SOLVE]: Ratio of norm of scaled residuals (reported); 9.623133e-07 
I[K_AUG]...	[MA57_SOLVE]: The norm of residuals is larger than max ratio(computed)
I[K_AUG]...	[MA57_SOLVE]: Attempting to reduce residuals

W[K_AUG]...	[MA57_SOLVE]Attempting to make dc > 0. (Jacobian regularization)
I[K_AUG]...	[INERTIA_STRATEGY]Inertia check successful.
I[MA57]...	[MA57_FACTOR]
W[MA57]...	[MA57_FACTOR]always_pert_jacobian is on before fact
I[K_AUG]...	[MA57_FACTOR]n_neig = 3606
I[K_AUG]...	[INERTIA_STRATEGY]Inertia check successful.
I[MA57]...	[MA57_RATIO]
I[K_AUG]...	[MA57_SOLVE]Accuracy improvement tries 0, Trial pivtol 1e-06, Residual ratio 9.623132789624708e-07.
I[K_AUG]...	[MA57_SOLVE]Accuracy is not improving.

W[K_AUG]...	[MA57_DRIVER]Inertia check OK neig=3606, (neig == m).
I[K_AUG]...	[K_AUG_ASL]Linear solver done. 
I[K_AUG]...	[K_AUG_ASL]Timings..Ev&Assem 0.106, Fact 0.227, Overall 0.237.
WARNING: The covariance matrix is not positive semi-definite.

Compare the Covariance Matrix Results¶

extract_trace_covariance(cov_red, "Reduced Hessian")

extract_trace_covariance(cov_finite, "Finite Difference")

extract_trace_covariance(cov_auto, "Automatic Differentiation")


The covariance matrix from the Reduced Hessian method is:
                    Ua            Ub       inv_CpH       inv_CpS
Ua       1.866396e-10  3.970963e-11  1.763796e-09 -8.017280e-08
Ub      -4.894069e-12  1.203907e+04  8.061843e+03 -2.361608e+06
inv_CpH  1.733479e-09  8.061843e+03  5.398533e+03 -1.581427e+06
inv_CpS -7.153954e-08 -2.361608e+06 -1.581427e+06  4.632578e+08
The trace of the covariance matrix from the Reduced Hessian method is: 4.633e+08

The covariance matrix from the Finite Difference method is:
                    Ua            Ub       inv_CpH       inv_CpS
Ua       1.141433e-09  5.223740e-02  3.497822e-02 -1.025104e+01
Ub       5.223740e-02  2.855139e+06  1.911805e+06 -5.602909e+08
inv_CpH  3.497822e-02  1.911805e+06  1.280147e+06 -3.751714e+08
inv_CpS -1.025104e+01 -5.602909e+08 -3.751714e+08  1.099512e+11
The trace of the covariance matrix from the Finite Difference method is: 1.100e+11

WARNING: The covariance matrix from Automatic Differentiation is not positive semi-definite.
                    Ua            Ub       inv_CpH       inv_CpS
Ua       1.857019e-10 -6.253419e-08 -4.015440e-08  1.219855e-05
Ub      -6.259810e-08  1.784463e+06  1.194877e+06 -3.501819e+08
inv_CpH -4.019351e-08  1.194877e+06  8.000902e+05 -2.344819e+08
inv_CpS  1.221150e-05 -3.501819e+08 -2.344819e+08  6.871948e+10
The trace of the covariance matrix from the Automatic Differentiation method is: 6.872e+10

Takeaways from the Covariance Matrix Results¶

Significant differences in the covariance matrix indicates non-estimable model parameters
When the model structure is well posed such that all the model parameters are estimable, the covariance matrix obtained from the three methods should be similar. Significant differences in the covariance matrix indicates that some of the model parameters are not estimable.
Large variance in the estimated model parameters indicates poorly estimable parameters
Parameters like Ub, inv_CpH, and especially inv_CpS show very large variances across all methods, suggesting:
- High correlation between parameters
- Limited information in the data to estimate the parameters independently
- Lack of unique parameter estimates under the current model structure