Representing Model Discrepancy in Bound-to-Bound Data Collaboration

TL;DR

This paper enhances the B2BDC framework by explicitly modeling and incorporating model discrepancy through linear basis functions, improving uncertainty quantification and dataset consistency in deterministic models.

Contribution

It introduces a method to explicitly represent model discrepancy within B2BDC using basis functions, extending the framework's capability for uncertainty quantification.

Findings

Incorporates model discrepancy into B2BDC framework.

Provides formulas for modified predictions including discrepancy.

Generalizes dataset consistency to account for discrepancy.

Abstract

We extended the existing methodology in Bound-to-Bound Data Collaboration (B2BDC), an optimization-based deterministic uncertainty quantification (UQ) framework, to explicitly take into account model discrepancy. The discrepancy was represented as a linear combination of finite basis functions and the feasible set was constructed according to a collection of modified model-data constraints. Formulas for making predictions were also modified to include the model discrepancy function. Prior information about the model discrepancy can be added to the framework as additional constraints. Dataset consistency, a central feature of B2BDC, was generalized based on the extended framework.

Tables

Table 1. Table 1: Results of the dataset consistency analysis

$ϵ$	$n$
	0	1	2	3	4
0.05	inconsistent	inconsistent	inconsistent	$k^{\ast }\in \widetilde{ℱ}$	$k^{\ast }\in \widetilde{ℱ}$
0.10	inconsistent	$k^{\ast }\notin \widetilde{ℱ}$	$k^{\ast }\in \widetilde{ℱ}$	$k^{\ast }\in \widetilde{ℱ}$	$k^{\ast }\in \widetilde{ℱ}$

Table 2. Table 2: Projection of the joint feasible set, computed with the inadequate model defined in Equation   3.4 and different model discrepancy functions defined in Equation   3.6 , onto coordinate directions of model parameter k 𝑘 k and discrepancy-function coefficients { c i } i = 1 n superscript subscript subscript 𝑐 𝑖 𝑖 1 𝑛 \{c_{i}\}_{i=1}^{n} , as well as computed volume ratio of the joint feasible set to the box made of the projected intervals. The symbol ∅ \varnothing represents an empty posterior uncertainty due to dataset inconsistency.

$n$	$k$	$c_0$	$c_1$	$c_2$	$c_3$	Volume ratio
$ϵ=0.05$
	[0.24, 0.25]∗
0	$\varnothing$
1	$\varnothing$	$\varnothing$
2	$\varnothing$	$\varnothing$	$\varnothing$
3	[0.20, 0.30]	[-0.02, 0.03]	[-0.08, 0.06]	[-0.05, 0.00]		$5.2\times {10}^{-3}$
4	[0.20, 0.30]	[-0.02, 0.03]	[-0.12, 0.06]	[-0.13, 0.14]	[-0.05, 0.03]	$2.3\times {10}^{-4}$
$ϵ=0.10$
	[0.24, 0.26]∗
0	$\varnothing$
1	[0.20, 0.21]	[0.00, 0.01]				$5.2\times {10}^{-1}$
2	[0.20, 0.30]	[-0.01, 0.04]	[-0.11, 0.01]			$4.9\times {10}^{-2}$
3	[0.20, 0.30]	[-0.02, 0.04]	[-0.13, 0.06]	[-0.05, 0.03]		$1.2\times {10}^{-2}$
4	[0.20, 0.30]	[-0.04, 0.05]	[-0.19, 0.19]	[-0.31, 0.22]	[-0.07, 0.09]	$2.0\times {10}^{-4}$

Table 3. Table 3: Design conditions for the training data.

Design index	$s_1$	$T$ (K)	$s_2$	$P$ (atm)	$s_3$	$\varphi$
1	1	1200	1	10	1	1.25
2	1	1200	1	10	-1	0.75
3	1	1200	-1	1	1	1.25
4	1	1200	-1	1	-1	0.75
5	-1	1600	1	10	1	1.25
6	-1	1600	1	10	-1	0.75
7	-1	1600	-1	1	1	1.25
8	-1	1600	-1	1	-1	0.75
9	0	1370	0	3.2	0	1
10	1.215	1170	0	3.2	0	1
11	-1.215	1660	0	3.2	0	1
12	0	1370	1.215	12.8	0	1
13	0	1370	-1.215	0.78	0	1
14	0	1370	0	3.2	1.215	1.3
15	0	1370	0	3.2	-1.215	0.7

Table 4. Table 4: Tested model discrepancy functions.

Model discrepancy	Number of basis function $n$
No $\delta$	0
$\delta =c_0$	1
$\delta =c_0+{\sum }_{i=1}^3{c}_i{s}_i$	4
$\delta =c_0+{\sum }_{i=1}^3{c}_i{s}_i+{\sum }_{i,j=1;i\le j}^3{c}_{i,j}{s}_i{s}_j$	10

Table 5. Table 5: Results of dataset consistency and the distance between true model parameter value λ ∗ superscript 𝜆 \lambda^{*} and the feasible set.

$n$	Dataset consistency	$d_{{\lambda }^{\ast }}$
	$ϵ=0.01$
0	Inconsistent	—
1	Inconsistent	—
4	Consistent	0.167
10	Consistent	0.047
	$ϵ=0.005$
0	Inconsistent	—
1	Inconsistent	—
4	Inconsistent	—
10	Consistent	0.072

Table 6. Table 6: Scenario parameter values for model prediction.

Case index	Prediction	$s_1$	$T$ (K)	$s_2$	$P$ (atm)	$s_3$	$\varphi$
1	Interpolation	-0.6	1500	0.4	5	0	1
2	Extrapolation	-1.67	1800	0.4	5	0	1
3	Extrapolation	-1	1600	1.16	12	0	1
4	Extrapolation	-1	1600	0.4	5	1.6	1.4
5	Extrapolation	-1.67	1800	1.16	12	1.6	1.4

Table 7. Table 7: Projection of the joint feasible set, computed with linear δ 𝛿 \delta at ϵ = 0.01 italic-ϵ 0.01 \epsilon=0.01 , onto coordinate directions of each discrepancy-function coefficient.

Coefficient	Posterior uncertainty bounds
$c_0$	[-0.139, 0.112]
$c_1$	[-0.018, 0.005]
$c_2$	[-0.015, -0.006]
$c_3$	[-0.042, -0.013]

Table 8. Table 8: Detailed and reduced H 2 -O 2 reaction sets and associated parameters of the rate coefficients, A ​ T n ​ e − E / R ​ T 𝐴 superscript 𝑇 𝑛 superscript 𝑒 𝐸 𝑅 𝑇 AT^{n}e^{-E/RT} , in the units of cm 3 , mol, s, cal, K (from [ 52 ] ).

Reduced		Reactions	$A$	$n$	$E$
$✓$	1	H + O2 = O + OH	$2.65\times {10}^{16}$	-0.6707	17041
$✓$	2	O + H2 = H + OH	$3.87\times {10}^4$	2.7	6260
$✓$	3	OH + H2 = H + H2O	$2.16\times {10}^8$	1.51	3430
	4	OH + OH = O + H2O	$3.57\times {10}^4$	2.4	-2110
	5111Collision efficiency: Ar = 0.63.	H + H + M = H2 + M	$1.00\times {10}^{18}$	-1.0	0
		H + H + H2 = H2 + H2	$9.00\times {10}^{16}$	-0.6	0
		H + H + H2O = H2 + H2O	$6.00\times {10}^{19}$	-1.25	0
	6222Collision efficiencies: H2 = 2.4, H2O = 15.4, Ar = 0.83.	O + O + M = O2 + M	$1.20\times {10}^{17}$	-1.0	0
	7333Collision efficiencies: H2 = 2, H2O = 12, Ar = 0.7.	O + H + M = OH + M	$4.71\times {10}^{18}$	-1.0	0
	8444Collision efficiencies: H2 = 0.73, H2O = 3.65, Ar = 0.38.	H + OH + M = H2O + M	$2.20\times {10}^{22}$	-2.0	0
$✓$	9555Collision efficiencies: H2O = 12, Ar = 0.53; Troe parameters: $a$ = $0.5$, $T^{\ast \ast \ast }$ = ${10}^{-30}$, $T^{\ast }$ = ${10}^{30}$, $T^{\ast \ast }$ = ${10}^{100}$.	H + O2 + M = HO2 + M	$5.75\times {10}^{19}$	-1.4	0
		H + O2 = HO2	$4.65\times {10}^{12}$	0.44	0
	10	H + HO2 = O + H2O	$3.97\times {10}^{12}$	0.0	671
$✓$	11	H + HO2 = H2 + O2	$2.99\times {10}^6$	2.12	-1172
	12	H + HO2 = OH + OH	$8.40\times {10}^{13}$	0.0	635
	13	O + HO2 = OH + O2	$2.00\times {10}^{13}$	0.0	0
	14	OH + HO2 = H2O + O2	$2.89\times {10}^{13}$	0.0	-497
	15	HO2 + HO2 = H2O2 + O2	$1.30\times {10}^{11}$	0.0	-1630
		HO2 + HO2 = H2O2 + O2	$4.20\times {10}^{14}$	0.0	12000
	16666Collision efficiencies: H2 = 2, H2O = 6, Ar = 0.67; Troe parameters: $a$ = $1.0$, $T^{\ast \ast \ast }$ = ${10}^{-30}$, $T^{\ast }$ = ${10}^{30}$, $T^{\ast \ast }$ = ${10}^{30}$.	OH + OH + M = H2O2 + M	$1.46\times {10}^{11}$	0.868	-8548
		OH + OH = H2O2	$8.71\times {10}^9$	0.869	-2191
	17	H + H2O2 = H2O + OH	$1.00\times {10}^{13}$	0.0	3600
	18	H + H2O2 = HO2 + H2	$1.21\times {10}^7$	2.0	5200
	19	O + H2O2 = HO2 + OH	$9.63\times {10}^6$	2.0	4000
	20	OH + H2O2 = H2O + HO2	$1.74\times {10}^{12}$	0.0	318
		OH + H2O2 = H2O + HO2	$7.59\times {10}^{13}$	0.0	7272
	213	O + OH + M = HO2 + M	$8.00\times {10}^{15}$	0.0	0

Equations

$$y=\mathcal{R}(s)+ϵ,$$

$$\mathcal{R}(s)=M(x^{\ast },s)+\delta (s),$$

$$\mathcal{F}=\{x|x\in \mathcal{H},L_e\le M(x,s_e)\le U_e,e=1,2,\dots ,N\}.$$

$$C_{\text{D}}:=$$

$$L_e+\gamma (\frac{U_e-L_e}{2})\le M_e(x)\le U_e-\gamma (\frac{U_e-L_e}{2}),$$

$$e=1,2,\dots ,N.$$

$$L_p=\underset{x\in \mathcal{F}}{\min }M(x,s_p),$$

$$U_p=\underset{x\in \mathcal{F}}{\max }M(x,s_p).$$

$$L_e\le M(x,s_e)+\delta (s_e)\le U_e,e=1,2,\dots ,N.$$

$$\delta (s)=\sum _{i=1}^n{c}_i{\Phi }_i(s),$$

$$L_e\le M(x,s_e)+\sum _{i=1}^n{c}_i{\Phi }_i(s_e)\le U_e,e=1,2,\dots ,N.$$

$${\mathcal{F}}_{\delta }=\{(x,c)|x\in \mathcal{H},c\in {\mathcal{H}}_c,L_e\le M(x,s_e)+\sum _{i=1}^n{c}_i{\Phi }_i(s_e)\le U_e,e=1,2,\dots ,N\},$$

$$\tilde{\mathcal{F}}=\{x|\exists \hat{c}:(x,\hat{c})\in {\mathcal{F}}_{\delta }\},$$

$$L_p=\underset{(x,c)\in {\mathcal{F}}_{\delta }}{\min }M(x,s_p)+\sum _{i=1}^n{c}_i{\Phi }_i(s_p),$$

$$U_p=\underset{(x,c)\in {\mathcal{F}}_{\delta }}{\max }M(x,s_p)+\sum _{i=1}^n{c}_i{\Phi }_i(s_p).$$

$${\delta }^1(s)=\sum _{i=1}^n{c}_i{\Phi }_i(s)$$

$${\delta }^2(s)={\delta }^1(s)+\sum _{j=1}^{n^{\prime }}{c}_j^{\prime }{\Phi }_j^{\prime }(s).$$

$$\sum _{i=1}^n{c}_i{\varphi }_i(s_j)>0,j=1,2,\dots .$$

$$L_{\delta }=\underset{(x,c)\in {\mathcal{F}}_{\delta }}{\min }\sum _{i=1}^n{c}_i{\Phi }_i(s_p),$$

$$U_{\delta }=\underset{(x,c)\in {\mathcal{F}}_{\delta }}{\max }\sum _{i=1}^n{c}_i{\Phi }_i(s_p).$$

$$F=-kz,$$

$$\frac{d^{2}z}{dt}=-kz-b\frac{dz}{dt},$$

$${z|}_{t=0}=z_0=-1.5,$$

$${\frac{dz}{dt}|}_{t=0}=v_0=1,$$

$$z^{\ast }(k,t)=e^{-bt/2}\left[\frac{v_0+0.5b{z}_0}{\sqrt{k-b^2/4}}\sin \left(\sqrt{k-b^2/4}t\right)+z_0\cos \left(\sqrt{k-b^2/4}t\right)\right].$$

$$z(k,t)=\frac{v_0}{\sqrt{k}}\sin \left(\sqrt{k}t\right)+z_0\cos (\sqrt{k}t).$$

$$z_e=z^{\ast }(k^{\ast },t_e)+ϵu_e,$$

$$u_e\sim \mathcal{U}(-1,1),e=1,2,\dots ,20.$$

$$\delta (t)=\sum _{i=1}^n{c}_{i-1}{t}^{i-1},n=1,2,3,4,$$

$$\frac{1}{N+1}(\left|\delta (t_p)\right|+\sum _{i=1}^N\left|\delta (t_i)\right|)\le \alpha ,$$

$$s_1=$$

$$T_{\text{center}}=1370\text{K},T_{\text{low}}=1200\text{K},T_{\text{high}}=1600\text{K},$$

$$s_2=$$

$$P_{\text{center}}=3.2\text{atm},P_{\text{low}}=1\text{atm},T_{\text{high}}=10\text{atm},$$

$$s_3=$$

Full text

1 Introduction

During the past few decades, computational capabilities and data availability have seen substantial growth in many scientific and engineering fields. The growing demand for predictive models with quantifiable uncertainty has developed into an active research area, uncertainty quantification (UQ) [32]. Two principal objectives of UQ are inference of model parameters, also known as the inverse or calibration problem utilizing a set of known data (the training set), and model prediction outside such a set. Theories and methods have been developed from both statistical and deterministic perspectives. In the statistical perspective, specifically under a Bayesian inference framework [19], the prior distribution of model parameters is updated by the likelihood of observations through Bayes’ theorem. The produced posterior distribution is utilized for model parameter inference and predictions. In the deterministic perspective, for example the Bound-to-Bound Data Collaboration (B2BDC) method [10, 16, 17, 41] employed in this paper, posterior ranges on model parameters and predictions are obtained through inequality-constrained optimization.

The two perspectives, while complementing each others in some aspects, address essentially the same problems and produce comparable outcomes. Both methodologies emphasize the critical role of identifying the source of uncertainty [15]. In addition to Bayesian calibration methods, we also note that B2BDC shares conceptual similarities with Bayesian history matching [7] and methods under the set membership framework [6, 18, 30, 40]. Bayesian history matching retains a probabilistic interpretation of the data and defines a non-implausible region in the parameter space that contains all acceptable parameter vectors. Although both methods seek a region containing valid parameter vectors, Bayesian history matching uses stochastic emulators and improves the quality of the emulators through iterative updates of the non-implausible region, whereas B2BDC employs polynomial response surfaces and focuses primarily on evaluation of uncertainty in predictions [16, 17]. Similar to the setup in B2BDC, the set-membership framework in robust control describes uncertainty in prior information and data by constraints. Methods developed under this framework also formulate estimation/prediction of quantities of interest (QOIs) as regions defined by inequality constraints or solutions to optimization problems [30, 48]. However, these methods differ from B2BDC in that they often pursue a simply shaped approximation of the complex regions formed therein, e.g., a bounding ellipsoid [12, 40] or a minimum-volume bounding parallelotope [6]. B2BDC, on the other hand, uses polynomial surrogate models and handles the resulting nonconvexity through convex relaxation, leading to global guarantees on optimality [41]. A prior B2BDC analysis showed that approximation of the feasible set by a bounding ellipsoid or polytope may lead to overly conservative prediction results [38].

When analysis suggests disagreement among models and data, there are three possible causes: the model is correct and the data are flawed, the data are correct and the model is flawed, or both are flawed. In the present study, we focus on the second case, where we have more confidence in the data than the model. The first case has been a subject of our past studies [10, 23], and the third case will be left as a challenge for future work. In the work of Kennedy and O’Hagan [27], experimental observations are assumed to be noisy measurement of the underlying true process which represents reality,

[TABLE]

where $\epsilon$ is the measurement noise, $s$ are the scenario parameters, and $\mathcal{R}(s)$ represents reality. The scenario parameters are controllable properties known from the experimental setup, like the initial temperature and pressure, and can vary from experiment to experiment. In any scientific endeavor, knowledge of the true process is an idealization; a model, considered as tentatively entertained [4], may have a systematic error in prediction. Kennedy and O’Hagan [27] suggested to describe the uncertainty in the model form as an additive term, $\delta$, referred to as model inadequacy, to the model output,

[TABLE]

where $x^{*}$ are the underlying true calibration parameters. The model parameters are uncertain parameters intrinsic to the model, $M(\cdot)$, and share common values across all experiments.

This approach of compensating for model discrepancy has received substantial interest and following (see, e.g., [2, 5, 24, 26, 28, 35, 36, 46, 49]), some using a Gaussian process (GP) [39, 45] to represent $\delta(s)$ [5, 24, 35, 36, 46, 49] and others a functional decomposition [26, 28], while referring to model inadequacy as model discrepancy, model bias, model form uncertainty, model error, and model form error. Efforts have also been made to overcome the difficulty in identifying model discrepancy and model parameters individually, and to improve prediction performance at conditions different from the training data. For example, Brynjarsdóttir and O’Hagan [5] put constraints on the GP realization of model discrepancy at specific conditions derived from domain knowledge. Plumlee [35] argued that the prior distribution of model discrepancy should be orthogonal to the gradient of the model under certain assumptions. Wang et al. [49] estimated the model discrepancy and model parameters separately. Joseph and Melkote [26] constructed a statistical model of discrepancy in a sequential manner, limiting its contribution to the prediction. In Bayesian history matching, model discrepancy is taken into account by including a structured term in its statistical model of the relation between data and simulator’s output [20].

Our objective here is to resolve disagreement among models and data in the deterministic setting of B2BDC using the perspective of Equation 1.2. The optimization-based framework of B2BDC represents uncertainty by sets and has been successfully applied in several domains, including combustion science [10, 13, 17, 25, 37, 42] and engineering [34], atmospheric chemistry [43], quantum chemistry [8], and system biology [9, 11, 51]. In the present work, we expand the B2BDC formalism by adding a deterministic model-form discrepancy function to the constraints derived using model and data, conceptually following Kennedy and O’Hagan [27]. We start with a brief recount of B2BDC in Section 2.1, followed by reformulating the feasible set and prediction problems with model discrepancy in Section 2.2. Application of the proposed methodology is presented in Section 3 for two examples, a simple mass-spring-damper and a more realistic combustion system. Further interpretation of model discrepancy as a general consistency measure is discussed in Section 4. We conclude with summarizing comments in Section 5.

2 Methodology

Bound-to-Bound Data Collaboration (B2BDC) is a deterministic optimization-based framework for systematically combining models and experimental data with quantified uncertainties [10, 16, 17, 41]. In this framework, the uncertainties in experimental data are specified by intervals $[L_{e},U_{e}]$, where $L_{e}$ and $U_{e}$ are the lower and upper bounds assessed for the $e$-th quantity of interest (QOI). A prior uncertainty region of model parameters derived from domain knowledge is also given and denoted by $\mathcal{H}$. The collection of all provided information is referred to as a dataset.

2.1 B2BDC without model discrepancy

The experimental data are utilized to carve out a smaller region in $\mathcal{H}$, referred to as the feasible set and denoted by $\mathcal{F}$,

[TABLE]

The constraints $L_{e}\leq M(x,s_{e})\leq U_{e}$ are referred to as model-data constraints since they are derived by connecting model outputs with experimental bounds. The feasible set constitutes the posterior region of the model parameters which satisfy all model-data constraints. The dataset is referred to as consistent if its feasible set is non-empty and inconsistent otherwise. The feasibility is determined by calculating a numerical measure termed the scalar consistency measure (SCM) [10]. The quantity is defined by

[TABLE]

The dataset is consistent if $C_{\rm{D}}\geq 0$ as its non-negative value proves the existence of a parameter vector that satisfies all the constraints. For a consistent dataset, the predicted interval $[L_{p},\,U_{p}]$ for an unmeasured QOI can be calculated by solving

[TABLE]

and

[TABLE]

where $s_{p}$ are the corresponding scenario parameters.

An inconsistent dataset implies that the models, experimental data, and prior information are fundamentally incompatible with each other, making prediction essentially meaningless. In prior work, inconsistency was resolved by identifying likely offending experimental data. This was accomplished by computing sensitivities of the consistency measure [10] and/or minimal relaxations of the bounds to recover consistency [23]. In the present work, we focus on regaining dataset consistency and the ability to make predictions through scenario-dependent model corrections.

In practice, the optimization problems in Equations 2.2, 2.3 and 2.4 are often challenging to solve numerically because the underlying model is not convex and expensive to evaluate. The B2BDC approach is to bracket the solution to the optimization problems with two quantities, referred to as the inner bound and outer bound, whose computation is more tractable [11, 41]. Computation of the outer bound, in particular, motivates the use of a quadratic/polynomial or rational quadratic surrogate model for the underlying model at each scenario condition, i.e., $M(x,s_{e})$. The outer bound is then obtained by solving a semidefinite program [41]. The inner bound is a local solution to the original optimization problem found by nonlinear constrained optimization solvers. In our B2BDC toolbox [29], the inner bound is calculated with MATLAB function fmincon and the outer bound using SDP solver SeDuMi [47] or CVX package [21, 22]. Considering that many QOIs depend mainly on a small fraction of the model parameters, a phenomenon termed effective sparsity [3], a subset of the model parameters that have a significant effect on $M(x,s_{e})$, referred to as active parameters, is selected by a screening sensitivity analysis [14]. The coefficients of the surrogate model are computed by fitting simulations of $M(x,s_{e})$ at selected points, referred to as design points, in the active parameter space. Development of the active parameters results in a reduced computational cost of building and evaluating the surrogate models. In a typical B2BDC application, generation of surrogate models usually consumes the most CPU time. In addition to the savings provided by selecting active parameters for each QOI, the cost can be further reduced by running simulations at design points in parallel. Previous applications [42, 44] showed the computation of inner and outer bounds took at most a few minutes on a desktop computer for datasets consisted of 55 and 102 model parameters along with 167 and 77 QOIs, respectively.

2.2 B2BDC with model discrepancy

In practice, one may encounter an inconsistent dataset where the experimental data are more reliable than the model form (e.g., as in [33]). In such a case, resolution of the inconsistency may be sought by compensating the difference between the model and data through introduction of a discrepancy function. We adopt the definition of model discrepancy proposed by Kennedy and O’Hagan in Equation 1.2 and on that basis introduce a scenario-dependent discrepancy function $\delta(s)$ into the model-data constraints,

[TABLE]

We assume that the discrepancy takes the form of a linear combination of $n$ basis functions, $\{\Phi_{i}\}_{i=1}^{n}$,

[TABLE]

where $\{c_{i}\}_{i=1}^{n}$ are unknown coefficients, and $n=0$ refers hereafter to a zero discrepancy function (i.e., $\delta=0$). Similar forms of the discrepancy function have also been used by others [26]. Substitution of Equation 2.6 into Equation 2.5 results in,

[TABLE]

The linear form in Equation 2.6 is motivated by the fact that the existing tools in B2BDC are directly applicable with Equation 2.7 in the extended parameter space $(x,\,c)$: while $\{\Phi_{i}\}_{i=1}^{n}$ can be any set of nonlinear functions, the modified model-data constraints are linear in $\{c_{i}\}_{i=1}^{n}$. Representing model discrepancy by a linear combination of basis functions has been used by others [26]. However, our method is different from those in that it does not set as an objective to fit a particular model discrepancy function. Hence, orthogonality among basis functions, usually imposed to improve regression performance, is not required.

We define the joint feasible set in the extended parameter space of $x$ and $c$ by combining the prior uncertainty and modified model-data constraints,

[TABLE]

where $\mathcal{H}_{c}$ represents the prior uncertainty region of the discrepancy-function coefficients. The feasibility of $\mathcal{F}_{\delta}$ can be calculated by modifying the model-data constraints accordingly in the SCM formula in Equation 2.2. The projection of $\mathcal{F}_{\delta}$ on the model parameter space is,

[TABLE]

which represents the set of feasible model parameters after including the discrepancy function. When the joint feasible set is not empty, prediction at an unmeasured scenario $s_{p}$ can be obtained by solving the modified versions of Equations 2.3 and 2.4:

[TABLE]

and

[TABLE]

The proposed framework treats the data as the sum of model output and discrepancy function as implied by Equation 1.2. The joint feasible set, defined over the extended parameter space $(x,\,c)$, therefore represents only combinations of the model and discrepancy function that are consistent with the data. It is not possible to elucidate the model or the discrepancy function separately without further assumptions and/or additional information.

A challenge with the above approach is the choice of basis functions. A practitioner may simply choose a set with the least number of basis functions that resolves dataset inconsistency, following for instance the Akaike information criterion [1] of penalizing larger number of parameters, and thus prevent overfitting. If extra insight is provided, various forms of discrepancy function can be investigated before making the final decision based on considerations in addition to the requirement that dataset consistency is recovered.

The developed framework with model discrepancy expressed using Equation 2.6 has a general feature that, for a given dataset and a prediction QOI, the prediction interval becomes systematically wider if additional basis functions are included. To understand this, suppose two sets of basis functions are used in an analysis, with the second being a superset of the first. Let vector $c$ represent the coefficients for the shared basis functions $\{\Phi_{i}\}_{i=1}^{n}$ and $c^{\prime}$ the coefficient vector for the additional basis functions $\{\Phi^{\prime}_{j}\}_{j=1}^{n^{\prime}}$, i.e.,

[TABLE]

The corresponding joint feasible sets formed by Equation 2.8 are denoted by $\mathcal{F}_{\delta}^{1}$ and $\mathcal{F}_{\delta}^{2}$. Any feasible point $(x,\,c)\in\mathcal{F}_{\delta}^{1}$ is also feasible for $\mathcal{F}_{\delta}^{2}$ by setting $c^{\prime}$ to zero, i.e., $(x,\,c,\,c^{\prime})|_{c^{\prime}=0}\in\mathcal{F}_{\delta}^{2}$. Therefore, the posterior uncertainty interval of the QOI, predicted on $\mathcal{F}_{\delta}^{1}$, is always contained by that predicted on $\mathcal{F}_{\delta}^{2}$. The increased uncertainty in the prediction interval can depend on the prediction QOI, the dataset, and the selected basis functions, as will be demonstrated in Section 3.

Previous work (e.g., [5, 35]) has demonstrated the value of including prior knowledge of the model discrepancy function when applying statistical UQ methods. In B2BDC, this can be accomplished by incorporating additional constraints. For example, sign constraints on the discrepancy function, or its derivatives, can be enforced at specified scenario conditions by introducing linear inequalities in $c$. An example of forcing model discrepancy function to be positive at selected scenarios is

[TABLE]

The effect of such constraints is automatically propagated to predictions through augmenting the feasibility constraint in Equations 2.10 and 2.11. Another example of constraining the magnitude of model discrepancy function is given in Section 3.1.4.

The posterior uncertainty of the model discrepancy function at any specified scenario $s_{p}$ can be calculated by solving the prediction problems in Equations 2.10 and 2.11 with the objective replaced by $\delta(s_{p})$, i.e.,

[TABLE]

and

[TABLE]

Repeating this computation at various conditions in the scenario parameter space can identify regions where uncertainty in the model discrepancy function is large. In a similar manner, the posterior uncertainty can be calculated for each discrepancy-function coefficient $c_{i}$.

3 Numerical examples

Application of the extended B2BDC framework is demonstrated with two examples, an illustrative mass-spring-damper system and a realistic hydrogen combustion system. In each example, we started with a postulated “true” model to represent the underlying true process. An inadequate model for the analysis was then created by omitting some parts of the true model and a true calibration parameter value was selected. The developed framework was applied to the inadequate model and the following results are reported and discussed for different choices of basis functions.

Dataset consistency 2. 2.

If the dataset is consistent, a) the predicted intervals at interpolated and extrapolated conditions; b) whether the true process values are contained in the predicted intervals, and to a secondary point, whether the true calibration parameter $x^{*}$ is in the feasible set.

The computation was conducted for two levels of experimental uncertainty to provide a more comprehensive characterization of the developed method. To clarify the nomenclature, $\delta^{*}$ and $\delta$ are used to represent the true model discrepancy defined in Equation 1.2 and the linear combination in Equation 2.6, respectively. We also differentiate between interpolated and extrapolated predictions, where the former refers to $s_{p}$ lying within the training domain and the latter outside.

All the example scripts along with the general B2BDC software [29]—i.e., everything required to reproduce the results reported in this paper—can be found using the GitHub link https://github.com/B2BDC/.

3.1 One-dimensional mass-spring-damper system

The force, $F$, needed to extend or compress a spring by a small distance, $z$, is expressed using Hooke’s law

[TABLE]

where $k$ is a constant characteristic of the spring, its stiffness. We now consider a simple system: a ball attached to a spring, whose other end is fixed at a wall, sketched in Figure 1.

The ball has a mass $m=1$ and is placed initially at $z_{0}=-1.5$ with an initial velocity $v_{0}=1$. In addition to the force exerted by the spring, the motion of the ball is also affected by a damping force proportional to the ball’s velocity. Thus, the evolution of the ball’s displacement is described by

[TABLE]

where $b$ is the constant coefficient of the damping force and its value was set to 0.05. For a given $k$, displacement evolution of system described by Equation 3.2—the “true” model in this example—is the solution to a second order, constant coefficient, ordinary differential equation and has analytic form

[TABLE]

The “inadequate” model was constructed by neglecting the damping force (i.e., $b=0$), which results in the solution

[TABLE]

In both the true and inadequate models, the stiffness $k$ is the model parameter and the time $t$ is the scenario parameter. The true stiffness of the spring — the true calibration parameter value — was selected to be $k^{*}=0.25$ with the prior uncertainty interval $\mathcal{H}=[0.2,\,0.3]$. The real displacement was evaluated with $z^{*}(k^{*},t)$. The displacements computed by the two models with $k=k^{*}$ and their difference, the model discrepancy defined in Equation 1.2, are demonstrated in Figure 2 for $t\in[0,\,4]$.

The QOIs for this example were chosen to be the displacements of the ball at specified times $t_{e}$. The dataset was composed of twenty of these QOIs in the scenario region $t\in[0,\,3]$. For each QOI, an “observed” value was generated by adding uniform noise with a prescribed maximum magnitude $\epsilon$ to the true process value,

[TABLE]

QOI uncertainty bounds were generated by setting $L_{e}=z_{e}-\epsilon$ and $U_{e}=z_{e}+\epsilon$. The present analysis was performed with $\epsilon$ values of 0.05 and 0.1. Three prediction QOIs were generated for $t$ values of 1.5, 3.2, and 4. The first prediction case occurs at a scenario within the training-set domain of $[0,\,3]$ and is an interpolated prediction. The second and third cases occur at scenarios outside the training-set domain and are extrapolated predictions.

3.1.1 Dataset consistency and QOI prediction

We first considered the ideal situation where the true model, given by Equation 3.3, and the formulas in Section 2.1 were used in the B2BDC calculations. The prediction results are displayed in Figure 3. With this setup, the dataset is consistent with $k^{*}$ being feasible and the predicted intervals contain the true process values for both tested $\epsilon$’s. The length of the predicted intervals at each prediction scenario is shorter for a smaller value of $\epsilon$, indicating that more accurate measurements produce more accurate predictions.

We then moved to a more realistic situation where an inadequate model, given by Equation 3.4, was examined through the modified B2BDC framework, described in Section 2.2. Four different model discrepancy functions,

[TABLE]

were tested in addition to the case where $\delta=0$. The discrepancy function is a polynomial in $t$ of degree $n-1$.

The outcome of the dataset consistency analysis is summarized in Table 1. Examination of these results shows that the dataset is inconsistent for both values of $\epsilon$ when $n=0$, i.e., a zero model discrepancy function was used. For $\epsilon=0.05$, a quadratic $\delta$ is required to obtain dataset consistency. In this case, $k^{*}$ is also found to be feasible. For $\epsilon=0.1$, a constant $\delta$ is enough to achieve consistency. However, $k^{*}$ becomes feasible only after using a linear $\delta$.

The predicted QOI intervals are displayed in Figure 4 for $t=1.5$, $3.2$, and $4$. As expected, the prediction intervals with a higher order $\delta$ are wider for both $\epsilon$ values. In the cases where $\delta$ produced a consistent dataset for both $\epsilon$ values, a shorter prediction interval is observed with the smaller $\epsilon$. For $\epsilon=0.1$, the QOI interval predicted using a constant $\delta$ does not contain the true value at all time instances. With a linear $\delta$, the predicted interval contains the true value at all three time instances. The predicted interval contains the true value for all tested times and for both values of $\epsilon$ with a quadratic and cubic $\delta$.

3.1.2 Posterior bounds of model parameter and discrepancy-function coefficients

We now examine the posterior uncertainty bounds of model parameter $k$ and discrepancy-function coefficients $\{c_{i}\}_{i=1}^{n}$ obtained for different polynomial orders of $\delta$. These bounds are the 1-dimensional projection of the joint feasible set $\mathcal{F}_{\delta}$ onto coordinate directions. The volume ratio of the joint feasible set and the multidimensional box, whose sides are the posterior projections of the parameters, was calculated as the fraction of $10^{6}$ samples uniformly distributed in the box, that lay in $\mathcal{F}_{\delta}$. The results are presented in Table 2. For comparison, the posterior interval of $k$ obtained with the true model is also listed. The computed volume ratio results show that the joint feasible set becomes progressively smaller relative to the box as dimension increases.

The B2BDC analysis with the true model resulted in a significantly narrower posterior uncertainty interval for model parameter $k$ as compared to its prior; the interval in this case, by design, contains the true calibration parameter value. With the inadequate model and a constant model discrepancy function ($n=1$) at $\epsilon=0.1$, an even narrower posterior interval was obtained; however, the true value, $k^{*}$, is completely missed. With a higher order $\delta$, the posterior interval covers the same range as the prior. This outcome can be explained by considering two factors that affected the posterior interval of $k$.

Firstly, the inadequate model has a different functional dependency on the model parameter $k$, resulting in a problem specific change of the posterior bounds: feasible $k$ values for the true model can become infeasible for the inadequate model and vice versa. In the current example, this can be visually shown by comparing the displacement predicted using the true and inadequate models and its dependency on model parameter $k$, as demonstrated in Figure 5. Plotted in this figure are $z^{*}(k,t)$ and $z(k,t)$ computed for different $k$’s drawn from its prior interval along with the experimental bounds for the case where $\epsilon=0.1$. For given $t$, larger $k$ values produced larger displacements for both models. The resulting displacement bands (shown in cyan) cover similar vertical regions at smaller $t$ values but the band for the inadequate model gradually shifts upward with increasing magnitude at larger $t$ values comparing to that for the true model. For the last two observations, shown in the right inset plot, only a small portion of the band satisfies the QOI bounds. Note that this portion corresponds to smaller $k$ values. However, predictions with these smaller $k$ values invalidated at least one other QOI bound, motivating the use of $\delta$ to resolve inconsistency.

The second factor, as discussed in Section 2.2, is that inclusion of a higher order $\delta$ always results in a wider posterior interval. For a constant $\delta$ at $\epsilon=0.1$, the posterior interval widened from the empty set (a zero $\delta$) to an interval with finite length. With the constant $\delta$, feasible $(k,c_{0})$ can be found with $k$ limited to a very small region close to the prior lower bound. The red dashed curve in Figure 5 corresponds to the prediction with one of the feasible $(k,c_{0})$.

The posterior uncertainty intervals of $\{c_{i}\}_{i=1}^{n}$ also become systematically wider with a higher polynomial order of $\delta$, as expected. The enlarged posterior uncertainty intervals associated with individual parameters $k$ and $\{c_{i}\}_{i=1}^{n}$ are related to the phenomenon usually referred to in statistical literature (e.g., [5]) as confounding, manifesting itself in the presence of a strong correlation between model parameter(s) and model discrepancy despite their relatively wider marginal posterior distributions. We demonstrate this from a deterministic perspective by the plots shown in Figure 6, generated for the case of a linear $\delta$ at $\epsilon=0.1$. The plots display the joint feasible set of $k$, $c_{0}$ and $c_{1}$ along with its 2-dimensional projections. The three-dimensional plot clearly shows that the joint feasible set occupies only a small fraction of the enclosing cube. Inspection of the projections indicates that at a fixed $k$ value, the uncertainty in $c_{0}$ and $c_{1}$ is reduced, on average, to 66 and 46%, of their posterior ranges, and at fixed $c_{0}$, the uncertainty in $c_{1}$ is reduced, on average, to 26% of its posterior range.

3.1.3 Posterior uncertainty of model discrepancy

We now examine the upper and lower bounds of $\delta$ predicted at 1000 discrete time points, $t_{i}$, equally spaced in $[0,\,4]$. The bounds were calculated by solving problems in Equations 2.14 and 2.15. This region is divided into the interpolation zone ($t\in[0,\,3]$), where data exists, and the extrapolation zone ($t\in[3,\,4]$) for comparison. The uncertainty bands are shown in Figure 7 for quadratic and cubic $\delta$; they were generated by linearly interpolating adjacent upper and lower bounds.

Inspection of these results shows that the computed uncertainty bounds enclose $\delta^{*}$ in both the interpolation and extrapolation zones for both quadratic and cubic $\delta$. The width of the predicted uncertainty band is effectively constrained within the interpolation zone. The predicted uncertainty band starts to widen toward the end of the interpolation zone and diverges rapidly in the extrapolation zone. The observed divergence is more dramatic for a cubic $\delta$ than a quadratic $\delta$. The uncertainty band for a fixed $\delta$ is overall narrower with a smaller $\epsilon$ in both interpolation and extrapolation zones.

3.1.4 Additional constraints on model discrepancy

As discussed in Section 2.2, constraints derived from domain knowledge about the model discrepancy can be included in the B2BDC calculations. We illustrate this feature with the following example.

Let us assume that although we introduced a discrepancy function, we would still like to rely on the model more than on the introduced correction when making predictions. This idea reflects the general spirit of some existing work in the literature (e.g., [26]). This requirement can be attained by selecting among all feasible values of $\delta$ those that have their magnitude, averaged over data and prediction scenarios, below a prescribed threshold, $\alpha$,

[TABLE]

where $N$ is the number of experimental QOIs. The constraint was added to the joint feasible set constructed using Equation 2.8 and predictions were made with varying values of $\alpha$. The results for $\epsilon=0.1$ and cubic $\delta$ are shown in Figure 8.

As expected, the predicted interval increases for larger $\alpha$, reaching the value obtained without using eq. 3.7 eventually, as this additional constraint becomes inactive.

3.2 A hydrogen combustion model

In this section we apply the formalism described above to a hydrogen combustion model: a homogeneous adiabatic H2-air reaction system at constant volume. The evolution of the system states, i.e., species concentrations and temperature, is simulated numerically by solving a set of ordinary differential equations. The time derivatives of species concentrations and temperature are calculated based on the specified chemical reaction mechanism and the energy equation. Simulations with detailed (21 reactions [52]) and reduced (5 reactions [50]) mechanisms, listed in Appendix A, were considered as the true and inadequate models, respectively. The model parameters, denoted by $\lambda\in\mathbb{R}^{5}$, are logarithm of the multipliers associated with the five rate constants shared by both mechanisms, with their prior uncertainties taken from [52]. The true calibration parameter value was specified as $\lambda^{*}=\mathbf{0}$, where $\mathbf{0}$ is a vector of zeros.

The normalized scenario parameters, $s_{1}$, $s_{2}$ and $s_{3}$, were defined as

[TABLE]

where $T$, $P$ and $\phi$ are the initial temperature, initial pressure and equivalence ratio of the mixture, respectively. In this example, equivalence ratio is the ratio of hydrogen to oxygen concentrations in the initial mixture to that in the stoichiometric mixture. The use of inverse temperature and logarithm of pressure for defining $s_{1}$ and $s_{2}$ are common in the combustion field, e.g., [50].

A dataset was constructed using a second-order orthogonal design [31] over the scenario region $[-1,\,1]^{3}$. The scenario parameter values are listed in Table 3.

For each of the scenario conditions, the corresponding QOI was defined as the time when the hydrogen concentration drops to half of its initial value. This QOI was computed numerically from the simulated hydrogen concentration profile and denoted by $t_{1/2}^{*}(\lambda,T,P,\phi)$ and $t_{1/2}(\lambda,T,P,\phi)$ for the true and inadequate models, respectively. “Measured” QOIs, denoted by $t_{i}$, were generated by adding a relative noise to the true process values, specified as $t_{1/2}^{*}$ evaluated at the true calibration parameter value $\lambda^{*}$,

[TABLE]

The maximum noise magnitude, $\epsilon$, was assigned values of 0.01 and 0.005. As before, the uncertainty bounds were generated by computing $[(1-\epsilon)t_{i},(1+\epsilon)t_{i}]$. The QOI computed with the inadequate model has no analytic solution and a quadratic surrogate model $S_{i}$ was generated for each QOI such that $S_{i}(\lambda)\approx{\rm{ln}}(t_{1/2}(\lambda,T_{i},P_{i},\phi_{i}))$. As in the previous example, we consider a polynomial model discrepancy function (Table 4), but now with the scenario parameters $s_{1}$, $s_{2}$ and $s_{3}$ defined by Equation 3.8.

3.2.1 Dataset consistency and QOI prediction

Dataset consistency was calculated first and the results are given in Table 5.

Inspection of these results shows that with $\epsilon=0.01$, the dataset is inconsistent for both the zero $\delta$ and constant $\delta$ cases, and becomes consistent when linear and quadratic $\delta$ are used. After $\epsilon$ was lowered to 0.005, linear $\delta$ is insufficient to keep the dataset consistent. For cases where the dataset is consistent, the distance between the true calibration parameter value and the feasible set, denoted by $d_{\lambda^{*}}$ and defined as

[TABLE]

was calculated and the results are also reported in Table 5. In all these cases, the true calibration parameter value is not in the feasible set. Its distance from the feasible set is larger when lower order $\delta$ or smaller $\epsilon$ were used.

For cases where the dataset is consistent, model predictions were computed at one interpolated and four extrapolated scenarios, which are specified in Table 6.

The results are depicted in Figure 9.

Again, the lengths of predicted intervals are shorter with linear $\delta$ as compared to quadratic $\delta$. Similarly, smaller values of $\epsilon$ produced shorter predictions. At $\epsilon=0.01$, the predicted interval with a linear $\delta$ contains the true value for cases 1, 2 and 3, but underpredicts the target for cases 4 and 5. With a quadratic $\delta$, the predicted intervals contain the true values for all tested cases at both $\epsilon$ values.

3.2.2 Inference of model discrepancy

The projection of the feasible set on the parameter space of discrepancy function coefficients, $c$’s, describes not one but a set of discrepancy functions that are consistent with the data. The following analysis with the linear $\delta$ and $\epsilon=0.01$ shows an example of inferring model discrepancy from B2BDC calculations. The posterior uncertainty bounds of the discrepancy-function coefficients were calculated and the results are given in Table 7. The volume ratio of the joint feasible set to the multidimensional box, specified similarly as in Section 3.1, is $2.6\times 10^{-8}$ based on $10^{9}$ samples.

The results show that all feasible $c_{2}$ and $c_{3}$ are negative since the calculated posterior upper bounds are negative for these two coefficients. For the linear $\delta$, the coefficients are also the partial derivatives of the discrepancy function with respect to the scenario parameters, i.e.,

[TABLE]

All feasible $\delta$’s are therefore smaller at larger $s_{2}$ or $s_{3}$ values given other scenario parameters fixed.

The predicted interval of $\delta$, i.e., $[L_{\delta},\,U_{\delta}]$ from Equations 2.14 and 2.15, was then calculated in the $s_{2}$-$s_{3}$ ($P$-$\phi$) space at three fixed $s_{1}$ ($T$) values. The computed intervals were examined to determine the sign of feasible $\delta$’s at each specified scenario and the results are shown in Figure 10.

Similar patterns are observed for the three tested temperature values: except for a lower left triangle region where both pressure and equivalence ratio are relatively small, all feasible $\delta$’s are negative. As a result of dataset consistency, predictions made at the grey-region scenarios always add a negative correction to the model output, suggesting that the inadequate model systematically overpredicts the QOI. Combining the results that $c_{2}$ and $c_{3}$ are negative in the feasible set, the overprediction is likely to be stronger at larger $s_{2}$ and $s_{3}$ values, i.e., at higher pressures and equivalence ratios.

4 Discrepancy as a consistency measure

The above examples and discussion primarily focus on the impact of model discrepancy on prediction. The inclusion of discrepancy into B2BDC also provides the opportunity to calculate a more general consistency measure. For a given collection of basis functions $\{\Phi_{i}\}_{i=1}^{n}$, define this consistency measure as

[TABLE]

where the objective $f(\cdot)$ is a function of only the coefficients $\{c_{i}\}_{i=1}^{n}$ and reflects the “complexity” of the discrepancy function. In essence, Equation 4.1 asks the following question: what is the least complex discrepancy function required to recover dataset consistency? Different choices of $f(\cdot)$ and $\Phi_{i}(\cdot)$ produce different consistency measures. For example, defining the complexity and discrepancy functions as

[TABLE]

where ${n}=N$ and the $\{s_{i}\}_{i=1}^{N}$ are the dataset scenarios. Note that with this choice of discrepancy, the $e$-th model-data constraint in Equation 4.1 becomes

[TABLE]

This is exactly a version of the vector consistency measure presented in earlier work [23, Equation (4.4)]. Other choices of $f$, such as the integral over $\mathcal{H}_{c}$ of the squared discrepancy (or its squared derivatives, should each $\Phi_{i}$ be differentiable), can also be handled in the B2BDC framework.

Consistency measures formulated in this fashion gauge the disagreement between models and observations based on the “simplest” (or least complex) discrepancy required to render a dataset consistent. One potential application of this type of consistency measure is for model comparison. For a fixed set of basis functions $\{\Phi_{i}\}_{i=1}^{n}$, multiple models can be compared by evaluating Equation 4.1. This generalized consistency analysis with model discrepancy is currently being investigated and will be discussed in future work.

5 Conclusions

We examined the inclusion of model discrepancy as a linear combination of finite basis functions in B2BDC. The existing B2BDC framework was extended by reformulating the feasible set to include both model parameters and discrepancy-function coefficients; the prediction formulas were adjusted accordingly. Dataset consistency can be effectively recovered with the developed framework by increasing the complexity of the used model discrepancy function. The developed method offers a flexible construction of discrepancy function structure through the selection of basis functions; prior information on model discrepancy can be included naturally in the optimization problems as additional constraints. The confounding between model parameters and model discrepancy function in the posterior uncertainty, presented and discussed in statistical methodologies (e.g., [5]), was demonstrated from the deterministic perspective.

Appendix A Detailed and reduced hydrogen mechanisms

The detailed and reduced mechanisms used in the present work are listed in Table 8. The reduced mechanism is consisted of the 5 reactions, marked with bold font and check marks.