index.rst

Uncertainty Propagation Toolbox

The uncertainty_propagation module quantifies and propagates parametric uncertainties through optimization and simulation problem based on IDAES models. The module has two core features:

1. Given a parameter estimation model and data, calculate the least squares best fit parameter values and estimate their covariance.
2. Given a process model and the covariance for its parameters, estimate the variance of the optimal solution and the objective function.

Consider the optimization problem:

$$\begin{aligned} \begin{align*} \mathrm{minimize} ~~ & f(x,p) \\\ \mathrm{s.t.} \quad ~ & c(x,p) = 0 \\\ & x_{lb} \leq x \leq x_{ub} \end{align*} \end{aligned}$$

Here x ∈ ℝ^n × 1 are the decision variables, p ∈ ℝ^m × 1 are the parameters, f(x, p): ℝ^n × 1  × ℝ^m × 1 → ℝ is the objective function, c(x, p) = {c₁(x, p), …, c_k(x, p)} : ℝ^n × 1  × ℝ^m × 1 → ℝ^k × 1 are the constraints, and x_lb and x_ub are the lower and upper bounds, respectively.

Let x^* represent the optimal solution given parameters p^*. In many process systems engineering problems, p^* is estimated from data and has some uncertainty represented with covariance matrix Σ_p. This toolbox estimates the uncertainty in the optimal solution x^* and objective function value f(x^*, p^*) induced by uncertainty in p^*.

Based on a first-order error propagation formula, the variance of the objective function is:

$$\text{Var}[f(x^*,p^*)] = \left(\frac{\partial f}{\partial p} + \frac{\partial f}{\partial x}\frac{\partial x}{\partial p} \right) \Sigma_p \left(\frac{\partial f}{\partial p} + \frac{\partial f}{\partial x}\frac{\partial x}{\partial p} \right)^T$$

Likewise, the variance in constraint k is:

$$\text{Var}[c_k(x^*,p^*)] = \left(\frac{\partial c_k}{\partial p} + \frac{\partial c_k}{\partial x}\frac{\partial x}{\partial p} \right) \Sigma_p \left(\frac{\partial c_k}{\partial p} + \frac{\partial c_k}{\partial x}\frac{\partial x}{\partial p} \right)^T$$

Note that Var[c_k(x^*, p^*)] ≈ 0 because the constraints remain feasible for a small perturbation in p, i.e., c(x, p) = 0.

All gradients are calculated with k_aug [1]. More specifically, $\frac{\partial f}{\partial p}, \frac{\partial f}{\partial x}, \frac{\partial c_1}{\partial p}, \frac{\partial c_1}{\partial x}, \ldots, \frac{\partial c_k}{\partial p}, \frac{\partial c_k}{\partial x}$ evaluated at (x^*, p) are computed via automatic differentiation whereas $\frac{\partial x}{\partial p}$ are computed via nonlinear programming sensitivity theory.

The covariance matrix Σ_p is either user supplied or obtained via regression (with Pyomo.ParmEst).

Dependencies

k_aug [1] is required to use uncertainty_propagation module. The k_aug solver executable is easily installed via the idaes get-extensions command.

Basic Usage

This toolbox has two core functions:

1. propagate_uncertainty: Given an IDAES (Pyomo) process model with parameters p and covariance Σ_p, estimate Var[f(x^*, p)].
2. quantify_propagate_uncertainty: Given an IDAES (Pyomo) regression model and data, first estimate parameters p and covariance Σ_p. Then given a second IDAES (Pyomo) process model, estimate Var[f(x^*, p)].

The following example shows the usage of quantify_propagate_uncertainty for a reaction kinetic problem from Rooney and Biegler [2]. Consider the mathematical model:

y_i(θ₁, θ₂, t_i) = θ₁(1 − e^{− θ₂t_i}), i = 1, ..., n

Here y_i is the concentration of the chemical product at time t_i and p = (θ₁, θ₂) are the fitted model parameters.

# Required imports
>>> from idaes.apps.uncertainty_propagation.uncertainties import quantify_propagate_uncertainty
>>> import pyomo.environ as *
>>> import pandas as pd

# Define rooney_biegler_model
>>> def rooney_biegler_model(data):
        model = ConcreteModel()
        model.asymptote = Var(initialize = 15)
        model.rate_constant = Var(initialize = 0.5)

        def response_rule(m, h):
            expr = m.asymptote * (1 - exp(-m.rate_constant * h))
            return expr
        model.response_function = Expression(data.hour, rule = response_rule)

        return model

# Define rooney_biegler_model_opt
>>> def rooney_biegler_model_opt(data):
        model = ConcreteModel()
        model.asymptote = Var(initialize = 15)
        model.rate_constant = Var(initialize = 0.5)
        model.obj = Objective(expr=model.asymptote*(1-exp(-model.rate_constant*10))
                              , sense=minimize)
        return model

# Define data
>>> data = pd.DataFrame(data=[[1,8.3],[2,10.3],[3,19.0],[4,16.0],[5,15.6],[7,19.8]],
                        columns=['hour', 'y'])

# Define variable_name
>>> variable_name = ['asymptote', 'rate_constant']

# Define SSE_obj_function
>>> def SSE_obj_function(model, data):
        expr = sum((data.y[i] - model.response_function[data.hour[i]])**2 for i in data.index)
        return expr

# Run quantify_propagate_uncertainty
>>> results = quantify_propagate_uncertainty(rooney_biegler_model, rooney_biegler_model_opt,  
                                             data, variable_name, obj_function)  

The Python function rooney_biegler_model generates a Pyomo regression model using the input Pandas DataFrame data. This model contains the attributes asymptote and rate_constant which are the parameters p to be estimated by minimizing the sum of squared errors (SSE). The list variable_name contains strings with these attribute names.

Similarly, the Python function rooney_biegler_model_opt returns a concrete Pyomo model which is the process optimization problem. This specific example has no degrees of freedom once p is specified; the objective function computes the product concentration at time t = 10.

The function quantify_propagate_uncertainty returns the object results which contains several important attributes:

• results.theta_names contains the names of parameters p
• results.theta contains the estimate values for parameters p
• results.gradient_f contains the gradient $\frac{\partial f}{\partial x},~ \frac{\partial f}{\partial p}$
• results.gradient_c contains the Jacobians $\frac{\partial c}{\partial x},~ \frac{\partial c}{\partial p}$
• results.dsdp contains the Jacobians $\frac{\partial x}{\partial p},~\frac{\partial p}{\partial p}$`
• results.propagation_f contains the estimate variance of the objective function

Important Notes:

1. The uncertain parameters p should be declared as Var in Pyomo.
2. The uncertain parameters p should not be fixed in Pyomo. Instead, set the upper and lower bounds equal. If they are fixed, k_aug will give an error message that the optimization problem has too few degrees of freedom.

Available Functions

idaes.apps.uncertainty_propagation.uncertainties

quantify_propagate_uncertainty

propagate_uncertainty

clean_variable_name

Examples

Two examples are provided to illustrate the detailed usage of uncertainty_propagation. In each case, a Jupyter notebook with explanations as well as an equivalent Python script is provided.

References

[1] David Thierry (2020). k_aug, https://github.com/dthierry/k_aug

[2] Rooney, W. C. and Biegler, L. T. (2001). Design for model parameter uncertainty using nonlinear confidence regions. AIChE Journal, 47(8), 1794-1804.