Error modeling for surrogates of dynamical systems using machine learning

TL;DR

This paper introduces a machine learning framework to model and correct errors in surrogate models of parameterized dynamical systems, improving accuracy and enabling statistical error analysis.

Contribution

The paper presents a novel machine learning approach that automatically models surrogate errors without feature selection, using local regression models for improved accuracy.

Findings

Enhanced accuracy in surrogate model predictions.

Effective statistical modeling of surrogate errors.

Applicable to nonlinear oil-water flow simulations.

Abstract

A machine-learning-based framework for modeling the error introduced by surrogate models of parameterized dynamical systems is proposed. The framework entails the use of high-dimensional regression techniques (e.g., random forests, LASSO) to map a large set of inexpensively computed `error indicators' (i.e., features) produced by the surrogate model at a given time instance to a prediction of the surrogate-model error in a quantity of interest (QoI). This eliminates the need for the user to hand-select a small number of informative features. The methodology requires a training set of parameter instances at which the time-dependent surrogate-model error is computed by simulating both the high-fidelity and surrogate models. Using these training data, the method first determines regression-model locality (via classification or clustering), and subsequently constructs a `local' regression model to predict the time-instantaneous error within each identified region of feature space. We consider two uses for the resulting error model: (1) as a correction to the surrogate-model QoI prediction at each time instance, and (2) as a way to statistically model arbitrary functions of the time-dependent surrogate-model error (e.g., time-integrated errors). We apply the proposed framework to model errors in reduced-order models of nonlinear oil--water subsurface flow simulations. The reduced-order models used in this work entail application of trajectory piecewise linearization with proper orthogonal decomposition. When the first use of the method is considered, numerical experiments demonstrate consistent improvement in accuracy in the time-instantaneous QoI prediction relative to the original surrogate model, across a large number of test cases. When the second use is considered, results show that the proposed method provides accurate statistical predictions of the time- and well-averaged errors.

Figures

Tables

Table 1. Table 1 : A subset of the features 𝒇 n superscript 𝒇 𝑛 {\boldsymbol{{f}}^{{n}}} used in EMML

No.	Feature	No.	Feature
1.	${\widetilde{𝐏}}_d{𝐳}^n$	2.	${\widetilde{𝐏}}_d{\acute{𝐳}}^i$
3.	${\widetilde{𝐏}}_d{𝐳}^{n-1}$	4.	${\widetilde{𝐏}}_d{\acute{𝐳}}^{i-1}$
5.	$\left({\text{e}}_{2d}^T𝚽\right)\left({𝐳}^n-{𝐳}^{n-1}\right)/\Delta t^n$	6.	$\left({\text{e}}_{2d}^T𝚽\right)\left({\acute{𝐳}}^i-{\acute{𝐳}}^{i-1}\right)/\Delta t^i$
7.	$u_d^n$	8.	${\acute{u}}_d^i$
9.	${𝐏}_d{\acute{𝐉}}^i{𝐏}_d^T$	10.	${𝐏}_d{\acute{𝐁}}^i{𝐏}_d^T$
11.	${𝐏}_d{\acute{𝐂}}^i{𝐏}_d^T$	12.	$\frac{⟨{𝐳}^n,{\acute{𝐳}}^i⟩}{{\parallel {𝐳}^n\parallel }_2{\parallel {\acute{𝐳}}^i\parallel }_2}$
13.	$\Delta t^n$	14.	$\Delta {\acute{t}}^i$
15.	${\mathrm{PVI}}^n$	16.	${\mathrm{PVI}}^i$
17.	$t^n$	18.	${\acute{t}}^i$
19.	${\widetilde{𝐏}}_{2d}{𝐳}^n-{\widetilde{𝐏}}_{2d}{\acute{𝐳}}^i$	20.	${\widetilde{𝐏}}_{2d}{𝐳}^{n-1}-{\widetilde{𝐏}}_{2d}{\acute{𝐳}}^{i-1}$

Table 2. Table 2 : Classification features 𝒇 c n superscript subscript 𝒇 c 𝑛 {{\boldsymbol{{f}}_{\rm{c}}^{{n}}}} corresponding to a production well in grid block d ∈ 𝒟 P 𝑑 subscript 𝒟 𝑃 d\in\mathcal{D}_{P}

No.	Feature	No.	Feature
1.	${\left(\sum _{k=1}^n{\left(u_d^k-{\acute{u}}_d^k\right)}^2\right)}^{1/2}$	2.	$\frac{1}{n}\sum _{k=1}^n\left({\left(u_d^k\right)}_{\mathrm{RL}}-{\acute{u}}_d^k\right)$
3.	${\left(\sum _{k=1}^n{\left(\left({\text{e}}_{2d-1}^T-{\text{e}}_{2d^{\prime }-1}^T\right)𝚽{\text{z}}^k\right)}^2\right)}^{1/2}$, $d^{\prime }\in {𝒟}_I$	4.	${\left(\sum _{k=1}^n{\left(\left({\text{e}}_{2d-1}^T𝚽\right)\left({\text{z}}^k-{\acute{\text{z}}}^k\right)\right)}^2\right)}^{1/2}$
5.	$\frac{1}{n}\sum _{k=1}^n\left({\text{e}}_{2d-1}^T-{\text{e}}_{2d^{\prime }-1}^T\right)𝚽{\text{z}}^k$, $d^{\prime }\in {𝒟}_I$	6.	$\frac{1}{n}\sum _{k=1}^n\left({\text{e}}_{2d-1}^T𝚽\right)\left({\text{z}}^k-{\acute{\text{z}}}^k\right)$

Table 3. Table 3 : EMML for QoI correction: additional test cases. Median value of the time-integrated errors in POD–TPWL and corrected solutions over cases in EMML test set 𝝁 ∈ 𝒮 EMML 𝝁 subscript 𝒮 EMML \boldsymbol{\mu}\in\mathcal{S}_{\text{EMML}}

Method	$E_i(o,{𝒟}_P)$	$E_i(w,{𝒟}_P)$	$E_i(w,{𝒟}_I)$
POD–TPWL ($i=\mathrm{RL}$)	4.5%	6.3%	6.6%
EMML correction ($i=$corr)	2.8%	4.4%	3.7%

Table 4. Table 4 : EMML for QoI correction: alternative EMML parameters. Median value of the time-integrated errors in POD–TPWL and corrected solutions over cases in EMML test set 𝝁 ∈ 𝒮 EMML 𝝁 subscript 𝒮 EMML \boldsymbol{\mu}\in\mathcal{S}_{\text{EMML}} with different EMML parameters

Method					$E_i(o,{𝒟}_P)$	$E_i(w,{𝒟}_P)$	$E_i(w,{𝒟}_I)$
					(%)	(%)	(%)
POD–TPWL ($i=\mathrm{RL}$)					4.49	6.34	6.62
EMML correction ($i=$corr)
Approach	$\tau$	$N_{\text{train}}$	Locality	Regression
1	1	30	classification	RF	2.78	4.39	3.67
1	0	30	classification	RF	2.79	4.38	3.68
1	1	15	classification	RF	3.58	5.63	4.12
1	1	30	classification	LS	3.33	5.67	11.32
1	1	30	clustering	RF	3.23	5.22	3.68
1	1	30	none	RF	3.26	5.18	3.67
2	1	30	classification	RF	2.78	5.62	11.06

Equations

$$\begin{array}{rl}{g}^n& :\mathbf{\mu }\mapsto {\mathbf{x}}^n,n=1,\dots ,N_t,\\ & :{\mathbb{R}}^{N_{\mu }}\to {\mathbb{R}}^{N_x},\end{array}$$

$$\begin{array}{rl}q& :\mathbf{P}\text{x}\mapsto q(\mathbf{P}\text{x}),\\ & :{\mathbb{R}}^{N_P}\to \mathbb{R}.\end{array}$$

$$\begin{array}{rl}{g}_{\mathrm{surr}}^n& :\mathbf{\mu }\mapsto {\mathbf{z}}^n,n=1,\dots ,N_t,\\ & :{\mathbb{R}}^{N_{\mu }}\to {\mathbb{R}}^{\ell },\end{array}$$

$$\begin{array}{rl}{q}_{\mathrm{surr}}& :\text{z}\mapsto q(\tilde{\mathbf{P}}\text{z}),\\ & :{\mathbb{R}}^{\ell }\to \mathbb{R},\end{array}$$

$$\begin{array}{rl}{\mathbf{f}}^n& :\mathbf{\mu }\mapsto {\mathbf{f}}^n(\mathbf{\mu }),n=1,\dots ,N_t,\\ & :{\mathbb{R}}^{N_{\mu }}\to {\mathbb{R}}^{1\times N_f}.\end{array}$$

$${\delta }_q^n(\mathbf{\mu }):=q^n(\mathbf{\mu })-q_{\mathrm{surr}}^n(\mathbf{\mu }),n=1,\dots ,N_t,$$

$${\delta }_q^n=r_q({\mathbf{f}}^n)+\epsilon ,n=1,\dots ,N_t.$$

$${\hat{\delta }}_q^n={\hat{r}}_q({\mathbf{f}}^n),n=1,\dots ,N_t.$$

$${\bar{\delta }}_q^n:=\frac{{\delta }_q^n}{q(\mathbf{P}{\mathbf{x}}^n)},n=1,\dots ,N_t,$$

$${\bar{\delta }}_q^n={\bar{r}}_q({\mathbf{f}}^n)+\epsilon ,n=1,\dots ,N_t,$$

$${\hat{\bar{\delta }}}_q^n={\hat{\bar{r}}}_q({\mathbf{f}}^n),n=1,\dots ,N_t.$$

$$q(\mathbf{P}{\mathbf{x}}^n)=\frac{q(\tilde{\mathbf{P}}{\mathbf{z}}^n)}{1-{\bar{\delta }}_q^n},n=1,\dots ,N_t,$$

$${\delta }_q^n=q(\tilde{\mathbf{P}}{\mathbf{z}}^n)\left(\frac{{\bar{\delta }}_q^n}{1-{\bar{\delta }}_q^n}\right),n=1,\dots ,N_t.$$

$${\hat{\delta }}_q^n=q(\tilde{\mathbf{P}}{\mathbf{z}}^n)\left(\frac{{\hat{\bar{\delta }}}_q^n}{1-{\hat{\bar{\delta }}}_q^n}\right),n=1,\dots ,N_t.$$

$$[{\mathbf{\delta }}_{ss}^n{]}_i$$

$$[{\mathbf{\delta }}_{ss}^n{]}_i=r_{ss,i}({\mathbf{f}}^n)+\epsilon ,n=1,\dots ,N_t,i=1,\dots ,N_P,$$

$$[{\hat{\mathbf{\delta }}}_{ss}^n{]}_i={\hat{r}}_{ss,i}({\mathbf{f}}^n),n=1,\dots ,N_t,i=1,\dots ,N_P.$$

$${\delta }_q^n$$

$${\hat{\delta }}_q^n$$

$$\begin{array}{rl}[{\bar{\mathbf{\delta }}}_{ss}^n{]}_i& :=\frac{[{\mathbf{\delta }}_{ss}^n{]}_i}{[\mathbf{P}{\mathbf{x}}^n{]}_i},n=1,\dots ,N_t,i=1,\dots ,N_P,\end{array}$$

$${\mathcal{T}}_{\text{EMML}}:=\{{\mathbf{\mu }}_{\text{train}}^1,\dots ,{\mathbf{\mu }}_{\text{train}}^{N_{\text{train}}}\}\subseteq {\mathbb{R}}^{N_{\mu }},$$

$$\{\left({\delta }_q^n(\mathbf{\mu }),{\mathbf{\delta }}_{ss}^n(\mathbf{\mu }),{\mathbf{f}}^n(\mathbf{\mu })\right){\}}_{\mathbf{\mu }\in {\mathcal{T}}_{\text{EMML}},n=1,\dots ,N_t}.$$

$$q_{\mathrm{corr}}^n(\mathbf{\mu }):=q_{\mathrm{surr}}^n(\mathbf{\mu })+{\hat{\delta }}_q^n(\mathbf{\mu }),$$

$$\{\left(h({\delta }_q^1(\mathbf{\mu }),\dots ,{\delta }_q^{N_t}(\mathbf{\mu })),h({\hat{\delta }}_q^1(\mathbf{\mu }),\dots ,{\hat{\delta }}_q^{N_t}(\mathbf{\mu }))\right){\}}_{\mathbf{\mu }\in {\mathcal{T}}_{\text{ROMES}}},$$

$$\frac{\partial }{\partial t}\left(\varphi {\rho }_j{S}_j\right)+\nabla \cdot \left({\rho }_j{\text{v}}_j\right)+{\rho }_j{\tilde{q}}_j=0,$$

$${\text{v}}_j=-{\lambda }_j\text{k}\nabla p,$$

$$\text{x}={\left[p_1{S}_1\cdots p_{N_c}{S}_{N_c}\right]}^T\in {\mathbb{R}}^{2N_c},$$

$$(q_j{)}_d=(T_{\mathrm{well}}{)}_d{\left({\lambda }_j\right)}_d\left(p_d-u_d\right),j=o,w\text{for}d\in {\mathcal{D}}_P,j=w\text{for}d\in {\mathcal{D}}_I.$$

$$\mathbf{u}={\left[u_{d_1}\cdots u_{d_{N_u}}\right]}^T\in {\mathbb{R}}^{N_u}.$$

$$\boldsymbol{\mu}=\left[\begin{array}[]{c}\textbf{u}^{1}\\ \vdots\\ \textbf{u}^{N_{t}}\end{array}\right]\in\mathbb{R}^{N_{\mu}},$$

$${\mathbf{g}}^n:=\text{g}\left({\mathbf{x}}^n,{\mathbf{x}}^{n-1},{\mathbf{u}}^n\right)=\text{0},n=1,\dots ,N_t,$$

$${\mathbf{g}}_{\text{L}}^n$$

Full text

1 Introduction

Computational simulation is being increasingly employed for real-time and many-query problems such as design, optimal control, uncertainty quantification, and inverse modeling. However, these applications require the rapid and repeated simulation of a (parameterized) computational model, which often corresponds to a discretization of partial differential equations (PDEs) that can be nonlinear, time dependent, and multiscale in nature. Accurate predictive models can therefore incur substantial computational costs. While recent advances in parallel computing have reduced simulation times for high-fidelity models, the rapid, repeated simulation of such models remains a bottleneck in many applications.

This computational challenge has motivated the development of a wide range of surrogate-modeling methods. Surrogate models—which can be categorized as data fits, lower-fidelity models, or reduced-order models (ROMs)—are approximations of high-fidelity models (HFMs) that aim to provide large computational savings while preserving accuracy. Unfortunately, these models often introduce non-negligible errors due to the assumptions and approximations employed in their construction, and these errors can have deleterious effects on the resulting analysis. Thus, to more rigorously employ surrogate models, it is critical to quantify the errors they introduce. Without reliable error quantification —which can be accomplished via statistical approaches, rigorous error bounds, error indicators, or error models—the accuracy of surrogate-model predictions is unknown, and the trustworthiness of the resulting analysis may be questionable.

For this reason, researchers have developed a variety of methods for quantifying surrogate-model errors. Data-fit surrogate models construct a deterministic function (e.g., polynomial fit [1], artificial neural network [2]), or stochastic process (e.g., Gaussian process/kriging [3, 4]) that explicitly approximates the mapping from the model input parameters to the model output quantity/quantities of interest (QoI). For this class of surrogate model, statistical approaches such as the $R^{2}$ value [5], cross validation [5, 6], confidence intervals, and prediction intervals [5, 6] can be applied to quantify surrogate error. When the data fit corresponds to a stochastic-process model, the prediction variance can be applied to quantify the uncertainty in the prediction directly [7, 8]. Although data-fit surrogates are nonintrusive to implement (they require only ‘black-box’ queries of the HFM), their predictions are not physics based, which can lead to inaccurate predictions, especially for high-dimensional input-parameter spaces. For this reason, many applications demand more sophisticated physics-based surrogates such as lower-fidelity or reduced-order models.

Lower-fidelity models are physics-based surrogates that apply simplifications to the original HFM, such as coarser discretizations, lower-order approximations or linearization, or they neglect some physics. In this case, one approach for quantifying the error would be to explicitly model the error via a ‘data-fit’ mapping between input parameters and lower-fidelity-model error. These approaches—which have been pioneered in the field of multifidelity design optimization—typically enforce ‘global’ zeroth-order consistency between the corrected low-fidelity-model QoI and the HFM QoI at training points [9, 10, 11, 12, 13, 14, 15], or ‘local’ first- or second-order consistency at trust-region centers [16, 17]. Unfortunately, in the context of dynamical systems (which we consider in this work), such corrections are only applicable to scalar-valued QoI (e.g., time-averaged quantities). Quantifying the error in time-dependent quantities has been largely ignored, with the exception of recent work that interpolates time-dependent error models in the input-parameter space [18]. In addition, the error may exhibit a complex or oscillatory dependence on the input parameters, which adds further challenges in high-dimensional input-parameter spaces and can cause the approach to fail [13, 19]. Alternatively, adjoint-based error estimation (i.e., dual-weighted-residual error indicators) can also be applied to approximate QoI errors in the context of coarse finite-element [20, 21, 22, 23], finite-volume [24, 25, 26], and discontinuous Galerkin [27, 28] discretizations. Unfortunately, for dynamical systems, dual-weighted residuals require an additional (time-local or time-global) dual solve, which can incur a non-negligible additional cost.

Projection-based ROMs apply projection to reduce the dimensionality of the equations governing the HFM. Typically, the low-dimensional bases are derived empirically by evaluating the HFM at training points or by performing other analyses, e.g., solving Lyapunov equations or computing a Krylov subspace. ROM error is typically estimated by deriving rigorous a priori and a posteriori error bounds for the state, QoI, or transfer function; such bounds have been derived for the reduced-basis method [29, 30, 31], system-theoretic approaches (e.g., balanced truncation, rational interpolation) [32], and proper-orthogonal-decomposition (POD) Galerkin [33, 34] and least-squares Petrov–Galerkin (LSPG) [35] methods; see, e.g., [36] for a review. However—for dynamical systems—such error bounds typically grow exponentially in time, causing the bound to significantly overpredict the error [37], which can limit the practical utility of these bounds. Note that the approaches developed for quantifying the low-fidelity-model error could also be adopted to quantify ROM errors [13, 38], as ROMs can be interpreted as lower-dimensional models with empirically derived basis functions.

To address the above issues in the context of quantifying ROM errors, Drohmann and Carlberg [19] devised the reduced-order-model error surrogates (ROMES) approach , which can be considered an error-modeling approach, as it constructs a Gaussian process that maps error indicators (e.g., error bounds, residual norms, dual-weighted residuals) produced inexpensively by the ROM to a distribution over the ROM QoI error. This study demonstrated that—even in the presence of high-dimensional input-parameter spaces—the ROM produces a small number of inexpensively computable error indicators that can be employed to derive accurate, low-variance predictions of the ROM error. Follow-on work also investigated the use of statistical modeling and regression methods to model ROM errors from indicators [39]. While promising, the ROMES approach requires the user to hand-select the error indicators that are most strongly correlated with the ROM error; i.e., feature selection is left to the user. This task can be challenging in general applications, as the user may not have strong a priori knowledge of which (small number of) features inform the error. Further, the ROMES method was demonstrated only on a stationary (i.e., steady-state) problem, and its extension to dynamical systems is nontrivial. Finally, application of the ROMES method to other physics-based surrogates (e.g., coarsened or upscaled models) is not obvious, as different error indicators will likely be required for such surrogates.

In this paper, we propose a machine-learning-based framework for modeling the error introduced by physics-based surrogate models of dynamical systems. The approach applies statistical techniques for high-dimensional regression (e.g., random forests) to map a large set of inexpensively computed quantities or features generated by the surrogate model to a prediction of the time-instantaneous surrogate-model QoI error. This method is referred to as error modeling via machine learning (EMML). In contrast to the ROMES method, the proposed EMML approach enables a large number of potential error indicators to be included in the candidate feature set and thus does not require the user to manually select a small number of features that inform the error, which can be challenging as mentioned above. Thus, feature selection is included in the process of constructing the error model, and extension to multiple types of physics-based surrogates is straightforward—we assume only that the surrogate produces a large set of features that can be mined for potential error indicators.

While the proposed framework can be applied to any surrogate model in principle, in this study we apply the method to model errors introduced by a TPWL ROM [40, 41, 42, 43] applied with a proper orthogonal decomposition (POD) basis and LSPG projection [44, 35], which we refer to hereon as POD–TPWL. We employ this ROM rather than other approaches (e.g., (D)EIM with POD–Galerkin, GNAT) because TPWL is less intrusive: it simply requires extracting linear operators (i.e., Jacobians of the residual with respect to the current state, previous state, and controls) from the HFM simulation code during the offline training stage, and it is entirely independent of the HFM simulation code during the online stage. This particular ROM was also employed in previous subsurface-flow studies involving oil–water and oil–gas compositional problems [45, 46, 47]. At each time step during test simulations, POD–TPWL performs linearization around the nearest training solution; LSPG projection is then applied—with a low-dimensional POD basis—to reduce the dimensionality of the linearized system. While the approach has been shown to yield $\mathcal{O}(10^{2})-\mathcal{O}(10^{3})$ speedups [42, 43, 45], POD–TPWL incurs non-negligible errors due to the approximations it introduces, namely (1) linearization error, and (2) ‘out-of-plane error’ [34] arising from employing a low-dimensional POD trial subspace, (3) ‘in-plane error’ arising from projection, and (4) error propagated from the previous time step. See [46] for further discussion of POD–TPWL error. We aim to apply the proposed EMML framework to model the QoI error resulting from these approximations.

Finally, we note that while machine learning and its application across various disciplines have been extensively studied, the use of machine learning within the domain of physics-based modeling and simulation is relatively new, although it is quite promising [48, 49]. In the context of improving surrogate-model predictions, researchers have recently investigated the use of machine learning techniques to identify (via classification) the spatial locations for high low-fidelity-model error [50]. Machine learning has also been used to quantify the inherent error with kriging data-fit surrogates [51] and to derive improved closure models in the context of computational fluid mechanics [52, 53, 54, 55, 56]. Machine learning was also applied to derive the source term for the transport of an intermittency variable while transitioning from laminar to turbulent flow [57]. Similarly, regression techniques (e.g., LASSO [58]) have been used to calibrate and infer uncertainties in viscosity-model coefficients for transonic flow applications [59]. We note that although machine learning has been used in a post-processing step to identify the physical regions of high surrogate-model error [50], to our knowledge, the direct approximation of error through construction of a regression model has not been pursued.

The paper proceeds as follows. In Section 2 we describe the general EMML framework. Following the definition of the error in the QoI, we introduce four methods to map this error to a set of inexpensively computed features using high-dimensional regression techniques. The particular HFM (subsurface flow) and surrogate model (POD–TPWL) considered as an application are discussed in Section 3. In Section 4, we present numerical results for modeling the POD–TPWL error in flow quantities driven by time-varying control variables over a large number of test cases. Different algorithmic treatments are also considered. A summary and suggestions for future work are provided in Section 5. Appendices A and B present descriptions of random-forest and LASSO regression.

2 General Problem Statement

In this section, we describe the overall EMML framework for error modeling. We begin by introducing both the high-fidelity model (HFM) and the surrogate model in a general setting.

2.1 Dynamical-system high-fidelity model

Given input parameters $\boldsymbol{\mu}\in\mathbb{R}^{N_{\mu}}$, we assume that the HFM generates a time history of states ${{\textbf{x}}}^{n}\in\mathbb{R}^{N_{x}}$, $n=1,\ldots,N_{t}$, and a scalar-valued output QoI $q$ that depends on the state, i.e.,

[TABLE]

Here, $\mathbf{P}\in\{0,1\}^{N_{P}\times N_{x}}$ with $N_{P}\leq N_{x}$ is a sampling matrix comprising selected rows of the identity matrix. This operator extracts the elements of the state vector required to compute the QoI $q$.

2.2 Dynamical-system surrogate model

We assume that the inexpensive surrogate model generates a time history of surrogate-model states ${\textbf{z}}^{n}\in\mathbb{R}^{\ell}$, $n=1,\ldots,N_{t}$ (ideally with $\ell\ll N_{x}$), given the input parameters $\boldsymbol{\mu}\in\mathbb{R}^{N_{\mu}}$. The QoI can be computed from the surrogate model states z using the function $q_{\rm{surr}}$. Our critical assumption is that the surrogate model also produces auxillary data in the form of a time history of ‘features’ ${\boldsymbol{{f}}^{{n}}}\in\mathbb{R}^{1\times{N_{f}}}$, $n=1,\ldots,N_{t}$, given these parameters $\boldsymbol{\mu}$. The surrogate model is described as follows:

[TABLE]

Here, $\widetilde{\mathbf{P}}\in\mathbb{R}^{N_{P}\times\ell}$ denotes a prolongation operator that transforms the surrogate model state into the elements of the high-fidelity state required to compute the QoI. The decision of what to include in the set of features is motivated by the underlying form of $g_{\rm{surr}}^{{n}}$, as discussed in Section 3.3. For notational simplicity, from hereon we use $\boldsymbol{{f}}^{{n}}$ in place of $\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu})$.

2.3 Error modeling

Our objective is to predict the error in the QoI at each time step $n$, which we define as

[TABLE]

where $\delta_{q}^{{n}}\in\mathbb{R}$, $q^{{n}}(\boldsymbol{\mu})\vcentcolon=q\circ\mathbf{P}g^{{n}}(\boldsymbol{\mu})$, denotes the QoI computed using the HFM at time instance $n$, and $q_{\rm{surr}}^{{n}}(\boldsymbol{\mu})\vcentcolon=q_{\rm{surr}}\circ g_{\rm{surr}}^{{n}}(\boldsymbol{\mu})$ denotes the QoI computed using the surrogate model at time instance $n$. We propose to approximate this QoI error as $\hat{\delta}_{q}^{{n}}\approx\delta_{q}^{{n}},\ n=1,\ldots,N_{t}$, by constructing error surrogates using high-dimensional regression methods developed in the context of machine learning.

In particular, we propose four regression-based approaches that construct a mapping from the surrogate-model features to a QoI-error prediction $\hat{\delta}_{q}^{{n}}$. While it is possible to construct an error prediction $\hat{\delta}_{q}^{{n}}$ as a function of only the input parameters $\boldsymbol{\mu}$, this approach can fail due to the oscillatory behavior of certain surrogate-model errors in the input space [13, 19], as discussed in the Introduction. For notational simplicity, hereon $\delta_{q}^{{n}}=\delta_{q}^{{n}}(\boldsymbol{\mu})$.

2.3.1 Method 1: QoI error.

The first method models the nondeterministic mapping $\boldsymbol{{f}}^{{n}}\mapsto\delta_{q}^{{n}}$ as a sum of a deterministic function $r_{q}:\mathbb{R}^{1\times{N_{f}}}\rightarrow\mathbb{R}$, and nondeterministic noise $\varepsilon$ as follows:

[TABLE]

Here, $\varepsilon$ is a zero-mean random variable that accounts for potentially unknown features, deficiencies in the model form of $r_{q}$, and the error introduced due to sampling variability. Thus, $\varepsilon$ denotes irreducible error induced by the error model Equation (7), and it may in principle depend on the features (i.e., $\varepsilon=\varepsilon(\boldsymbol{{f}}^{{n}})$); this heteroscedasticity can occur, for example, when the mean and variance of the predicted error are larger for larger-magnitude features. For simplicity, we neglect this dependence in the current study.

We note that employing the error model in Equation (7) allows us to (1) approximate the form of $r_{q}$ using data, which in turn enables us to express the error in the QoI as a function of the features only, and (2) account for the possibility that the same feature vector may yield different values of the QoI error, i.e., $\delta_{q}^{n}\neq\delta_{q}^{m}$ with $n\neq m$ but $\boldsymbol{{f}}^{n}=\boldsymbol{{f}}^{m}$.

Next, we construct a model $\hat{r}_{q}:\mathbb{R}^{1\times{N_{f}}}\rightarrow\mathbb{R}$ of the function $r_{q}$, such that $\hat{r}_{q}(\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu}))\approx r_{q}(\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu})),n=1,\ldots,N_{t}$, $\forall\boldsymbol{\mu}\in\mathbb{R}^{N_{\mu}}$. This model allows the error to be approximated as a function of the features as

[TABLE]

Note that we consider modeling $\delta_{q}^{{n}}$ as a prediction problem rather than as a time-series-analysis problem. This is because, in the problem under consideration, we perform numerical experiments for different input-parameter instances over the same time interval. Thus, we include time in the feature set, as described in Section 3.3. Implicitly, we assume that samples are independent and identically distributed (i.i.d.), with each sample corresponding to the quantity at a given time step.

2.3.2 Method 2: relative QoI error.

In many cases, the QoI errors $\delta_{q}^{{n}}(\boldsymbol{\mu})$ can exhibit a wide range of observed values. This can make the machine-learning task more challenging, as the associated regression model must be predictive across this entire range of values. To address this, we can instead apply regression to the relative QoI errors—which typically exhibit a narrower range of values—and subsequently approximate the QoI errors in a postprocessing step.

We define the relative QoI error at time step $n$ by

[TABLE]

and—following Method 1—express the mapping $\boldsymbol{{f}}^{{n}}\mapsto\bar{\delta}_{q}^{{n}}$ as

[TABLE]

where $\bar{r}_{q}:\mathbb{R}^{1\times{N_{f}}}\rightarrow\mathbb{R}$ is an unknown deterministic function. We model $\bar{r}_{q}$ by constructing an approximation $\hat{\bar{r}}_{q}:\mathbb{R}^{1\times{N_{f}}}\rightarrow\mathbb{R}$ such that $\hat{\bar{r}}_{q}(\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu}))\approx r_{q}(\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu}))$, $n=1,\ldots,N_{t}$, $\forall\boldsymbol{\mu}\in\mathbb{R}^{N_{\mu}}$. This allows the approximated relative error to be expressed as

[TABLE]

From Equation (6), it follows that

[TABLE]

which allows the QoI error to be related to the relative QoI error by

[TABLE]

Therefore, we can model the QoI error as $\hat{\delta}_{q}^{{n}}$ from the relative QoI regression model $\hat{\bar{\delta}}_{q}^{{n}}$ in a postprocessing step as

[TABLE]

2.3.3 Method 3: state error.

As an alternative to modeling the QoI error, we can model the error in the relevant state(s), and then use this quantity to approximate the QoI error. This method is advantageous when the QoI error exhibits a more complex behavior than the state error, as may be the case for highly nonlinear QoI. We denote the error in the sampled state at time step ${n}$ by

[TABLE]

where $\boldsymbol{\delta}_{ss}^{{n}}\in\mathbb{R}^{N_{P}}$ and $[\cdot]_{i}$ denotes the $i$th element of the vector-valued argument.

The next steps follow closely the derivation of Method 1. As before, we model the mappings $\boldsymbol{{f}}^{{n}}\mapsto[\boldsymbol{\delta}_{ss}^{{n}}]_{i}$ as

[TABLE]

where $r_{ss,i}:\mathbb{R}^{1\times{N_{f}}}\rightarrow\mathbb{R}$, $i=1,\ldots,N_{P}$, denote unknown deterministic functions that allow the state-variable errors to be computed as a function of the features. Analogous to Equation (7), we construct regression models $\hat{r}_{ss,i}:\mathbb{R}^{1\times{N_{f}}}\rightarrow\mathbb{R}$, such that $\hat{r}_{ss,i}(\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu}))\approx r_{ss,i}(\boldsymbol{{f}}^{{n}}(\boldsymbol{\mu}))$, $n=1,\ldots,N_{t}$, $i=1,\ldots,N_{P}$, $\forall\boldsymbol{\mu}\in\mathbb{R}^{N_{\mu}}$. This model allows the approximated state error to be expressed as

[TABLE]

Note that instead of pursuing multi-response multivariate regression, we execute $N_{P}$ independent multivariate regressions, i.e., we construct a unique and independent mapping $\hat{r}_{ss,i}$ for each of the $N_{P}$ sampled states. Because the QoI error can be expressed as

[TABLE]

we can model the QoI error from the modeled state error in a postprocessing step as

[TABLE]

2.3.4 Method 4: relative state error.

Finally, if the errors in a particular sampled state exhibit a wide range of observed values, we can construct a regression model for the relative state error, which we define as

[TABLE]

with $\bar{\boldsymbol{\delta}}_{ss}^{{n}}\in\mathbb{R}^{N_{P}}$, and subsequently use this model to predict the QoI error. Analogous to Methods 1–3, we construct a regression model, which maps the features ${\boldsymbol{{f}}^{{n}}}$ to the relative error in the sampled state $[\bar{\boldsymbol{\delta}}_{ss}^{{n}}]_{i}$, thus enabling the computation of the approximated relative error $[\hat{\bar{\boldsymbol{\delta}}}_{g}^{{n}}]_{i}$. Analogous to Method 2, we relate the relative error in the sampled state to the actual error. After using Method 3 or Method 4 or some combination thereof to compute the value of $[\hat{\boldsymbol{\delta}}_{ss}^{{n}}]_{i},~{}\forall i=1,\ldots,N_{P}$, we apply Equation (19) to determine the QoI error model $\hat{\delta}_{q}$.

2.4 Training data

Each of the four methods described above entails the use of a regression model to predict the output (response)—which corresponds to the actual or relative error in the QoI or in the sampled state—given the inputs (features). Constructing any such regression model relies on training data. We denote points in the input-parameter space used to collect these data as

[TABLE]

where $\boldsymbol{\mu}^{i}\in\mathbb{R}^{N_{\mu}}$, $i=1,\ldots,N_{\text{train}}$ denotes the $i$th EMML training instance of the input parameters, and $N_{\text{train}}$ denotes the number of training points.

Next, we simulate both the HFM and the surrogate model for input-parameter instances in the training set $\mathcal{T}_{\text{EMML}}$. This produces the EMML training data, which comprise errors $\delta_{q}^{{n}},\boldsymbol{\delta}_{ss}^{{n}}$, and features ${\boldsymbol{{f}}^{{n}}}$ over all time steps and training simulations, i.e.,

[TABLE]

As mentioned previously, we assume that the associated samples are i.i.d. In the following, we denote a general training error by $\delta^{n}(\boldsymbol{\mu})$, $\boldsymbol{\mu}\in\mathcal{T}_{\text{EMML}},\ n=1,\ldots,N_{t}$.

2.5 Regression-model locality

As an alternative to constructing a single global regression function, $r_{ss,i},\bar{r}_{ss,i},r_{q}$, or $\bar{r}_{q}$, that is valid over the entire feature space, we can instead construct multiple ‘local’ regression models that are tailored to particular physical regimes or feature-space regions with the intent of improving prediction accuracy. To realize this, we partition the training data into subsets corresponding to different feature-space regions and construct separate regression functions for each subset. We consider two methods for determining regression-model locality: classification and clustering.

Classification is a supervised machine learning technique [6] that predicts the label (or category) to which an observation belongs. In this work it entails constructing a statistical model from EMML training data containing samples whose category membership is known, along with categorization criteria for those samples. In the current context, we propose applying classification using ‘classification features’ $\boldsymbol{{f}}_{\rm{c}}\in\mathbb{R}^{1\times N_{f_{c}}}$ that may in general be different from the EMML features $\boldsymbol{{f}}$. We employ these features to identify the subsets of the EMML training data associated with different physical regimes of the problem, for which different regression models are appropriate. Then, given a new observation for which we require error prediction, we first identify its category using the classification model, and subsequently apply the associated local regression model for error prediction.

Clustering is an unsupervised machine learning method [6] that can be applied to partition the training data according to the (e.g., Euclidean) distance in feature space between the training samples. We propose partitioning feature space (or a lower-dimensional space embedded in feature space computed, e.g., via principal component analysis) according to the Voronoi diagram produced by the cluster means. In this work, we employ $k$-means clustering where the number of clusters $k$ is determined by identifying the elbow of the curve reporting the sum of squared errors as a function of the number of clusters [6]. Given a new observation, we first identify its cluster from the Voronoi diagram, and then apply the local regression model.

2.6 High-dimensional regression methods

While, in principle, standard regression methods (e.g., linear regression) could be used to construct error surrogates $r_{ss,i},\bar{r}_{ss,i},r_{q}$, and $\bar{r}_{q}$, such approaches may be ineffective when the number of candidate features ${N_{f}}$ produced by the surrogate model is large. This occurs, for example, when a projection-based reduced-order model is employed as the surrogate. This ineffectiveness arises due to (1) the lack of available guidelines for feature-subset selection, (2) the time-consuming and challenging nature of a priori identification of the relevant subset of features, and (3) the fact that the response–feature relationship may depend on nonlinear interactions between a large number of features. We therefore propose applying high-dimensional regression methods that incorporate automatic feature selection.

A wide range of methods—such as tree-based methods (gradient boosting, random forests), support vector machines, $K$-nearest neighbors, elastic-net, and artificial neural network—fits within this category. While the specific choice of regression technique depends on the problem at hand, we pursue two specific methods in this work: random-forest regression (RF) and LASSO (least absolute shrinkage and selection operator [58]) regression (LS). For completeness, Appendices A and B provide brief summaries of these two techniques.

2.7 Application of error models

We propose two practical ways to use a QoI-error prediction $\hat{\delta}_{q}^{{n}}$: (1) as a correction to the surrogate-model QoI prediction at a given time instance, or (2) as an error indicator to be used within the Gaussian-process-based ROMES framework [19] for statistically modeling arbitrary functions of the time-dependent surrogate-model error. These approaches are now described in turn.

2.7.1 QoI correction

The most obvious way in which to employ the QoI-error prediction $\hat{\delta}_{q}^{{n}}$ is simply to apply it as a correction to the time-instantaneous QoI computed using the surrogate model, i.e., employ

[TABLE]

as a corrected QoI at time instance ${n}$. Of course, our expectation is that the corrected QoI error has smaller magnitude than the surrogate QoI error, i.e., $|q^{{n}}(\boldsymbol{\mu})-q_{\rm{corr}}^{{n}}(\boldsymbol{\mu})|=|\delta_{q}^{{n}}(\boldsymbol{\mu})-\hat{\delta}_{q}^{{n}}(\boldsymbol{\mu})|<|q^{{n}}(\boldsymbol{\mu})-q_{\rm{surr}}^{{n}}(\boldsymbol{\mu})|=|\delta_{q}^{{n}}(\boldsymbol{\mu})|$.

2.7.2 QoI error modeling

Alternatively, we can adopt the perspective of the ROMES method [19], which aims to construct a statistical model of the surrogate-model error via Gaussian-process regression. The key insight of the method is that one particular type of surrogate—reduced-order models—produce inexpensively computable error indicators such as error bounds, residual norms, and dual-weighted residuals that correlate strongly with the surrogate-model error. The method exploits such error indicators by constructing a Gaussian process that maps the chosen error indicator to a (Gaussian) distribution over the true surrogate-model error.

In the present context, given an arbitrary function of the time-dependent surrogate-model error that we would like to predict, ${h}(\hat{\delta}_{q}^{1}(\boldsymbol{\mu}),\ldots,\hat{\delta}_{q}^{N_{t}}(\boldsymbol{\mu}))$, we propose employing the same function applied to the EMML-predicted surrogate QoI errors, ${h}(\delta_{q}^{1}(\boldsymbol{\mu}),\ldots,\delta_{q}^{N_{t}}(\boldsymbol{\mu}))$, as an error indicator in the ROMES framework. We expect this to perform well if the EMML QoI-error predictions $\hat{\delta}_{q}^{{n}}$, $n=1,\ldots,N_{t}$, are accurate representations of the true QoI errors $\delta_{q}^{{n}}$, $n=1,\ldots,N_{t}$.

Specifically, defining a probability space $(\Omega,\mathcal{F},P)$, we approximate the deterministic mapping $\boldsymbol{\mu}\mapsto{h}(\delta_{q}^{1}(\boldsymbol{\mu}),\ldots,\delta_{q}^{N_{t}}(\boldsymbol{\mu}))$ by a stochastic mapping ${h}(\hat{\delta}_{q}^{1}(\boldsymbol{\mu}),\ldots,\hat{\delta}_{q}^{N_{t}}(\boldsymbol{\mu}))\mapsto\hat{h}$, with $\hat{h}:\Omega\rightarrow\mathbb{R}$ a scalar-valued Gaussian random variable that can be considered a statistical model of the true error function ${h}(\delta_{q}^{1}(\boldsymbol{\mu}),\ldots,\delta_{q}^{N_{t}}(\boldsymbol{\mu}))$. The stochastic mapping is constructed via Gaussian-process regression using ROMES training data

[TABLE]

where $\mathcal{T}_{\text{ROMES}}\subset\mathbb{R}^{N_{\mu}}$ denotes the ROMES training set, which should be distinct from the EMML training set $\mathcal{T}_{\text{EMML}}$.

3 Application to subsurface flow

In this section, we first present the governing equations for a two-phase oil–water subsurface flow system, followed by the POD–TPWL reduced-order (surrogate) model used in this work. We then describe the specialization of the EMML components (error-modeling approach, feature design, training/test data, determining locality for the regression model) employed for this application. Please refer to [42] for a description of the oil–water flow equations and associated finite-volume discretization, and for a detailed development of POD–TPWL for such systems. The use of LSPG projection with POD–TPWL is described in [45, 46, 47].

3.1 High-fidelity model: two-phase oil–water system

The HFM for the two-phase oil–water problem entails statements of conservation of mass for the oil and water components combined with Darcy’s law for each phase. Assuming immiscibility (which means that components only exist in their corresponding phase), and neglecting capillary pressure and gravitational effects, the equations for phase $j$ can be written as

[TABLE]

where $j$=$o$ designates the oil phase and $j$=$w$ the water phase, $t$ is time, $\phi$ is porosity (void fraction within the rock), $\rho_{j}$ denotes phase density, $S_{j}$ is phase saturation (i.e., phase fraction), $\textbf{v}_{j}$ is the Darcy phase velocity, ${\tilde{q}}_{j}$ denotes the well phase flow rate per unit volume (${\tilde{q}}_{j}>0$ for production/sink wells and ${\tilde{q}}_{j}<0$ for injection/source wells), k denotes the permeability tensor (taken to be isotropic in the examples here), $\lambda_{j}=k_{rj}(S_{j})/\mu_{j}$ designates the phase mobility, with $k_{rj}$ the relative permeability to phase $j$ and $\mu_{j}$ the phase viscosity, and $p$ denotes pressure (note that $p_{o}=p_{w}=p$ because capillary pressure is neglected). We additionally have the saturation constraint $S_{w}+S_{o}=1$. For subsurface flow problems, we are often interested in predicting the phase flow rates for all production and injection wells.

The two-phase system described by Equation (22) is discretized using a finite-volume method with pressure $p$ and water saturation $S\vcentcolon=S_{w}$ in each grid block as the primary unknowns. Thus, $N_{P}=2$. Then, the (time-dependent) states can be represented as

[TABLE]

where $N_{{c}}$ denotes the number of grid blocks; thus, $N_{x}=2N_{{c}}$ for this application. At each time step, we consider the set of QoI to be the well phase flow rates $q_{j}$, at a subset of $N_{u}\ll N_{{c}}$ grid blocks referred to as the well blocks, which we represent by the indices $\mathcal{D}\vcentcolon=\{{d_{1}},\ldots,{d_{N_{u}}}\}\subset\{1,\ldots,N_{{c}}\}$ with $\mathcal{D}=\mathcal{D}_{P}\cup\mathcal{D}_{I}$ and $\mathcal{D}_{P}\cap\mathcal{D}_{I}=\emptyset$, where $\mathcal{D}_{P}\subseteq\mathcal{D}$ denotes the set of producer wells and $\mathcal{D}_{I}\subseteq\mathcal{D}$ denotes the set of injector wells.

We compute these flow-rate QoI using the standard Peaceman well model [60]:

[TABLE]

Here, subscript $d$ indicates the value of a variable at grid block $d$, $T_{\rm{well}}\in\mathbb{R}_{+}$ denotes the well index, which depends on the wellbore radius and well-block permeability and geometry ($T_{\rm{well}}$ is essentially the transmissibility linking the well to the well block), and $u\in\mathbb{R}_{+}$ denotes the prescribed wellbore pressure, also referred to as the bottom-hole pressure (BHP). Equation (24) is written for a discrete finite-volume model. Thus, $(q_{j})_{d}={\tilde{q}}_{j}V_{d}$, where $V_{d}$ is the volume of grid block $d$ $\left((q_{j})_{d}\text{ is of units volume/time}\right)$.

In the systems considered here, only water is injected. The output-function $q$ introduced in Equation (2) is defined by Equation (24), where the sampling matrix $\mathbf{P}$ simply extracts the pressure and saturation from the appropriate well block, i.e., the output corresponding to $(q_{j})_{d}$, $j=o,w$, employs a sampling matrix $\mathbf{P}_{d}=\left[{\textbf{e}_{2d-1}}\ {\textbf{e}_{2d}}\right]^{T}$, where ${\textbf{e}_{i}}$ denotes the $i$th canonical unit vector. Thus, the QoI depend on both pressure and saturation, i.e., $N_{P}=2$ for each QoI. Note that the treatment described here is for cases where a particular well penetrates only a single grid block. Our procedures can be generalized, however, for multiblock well penetrations and for cases where rates (rather than BHPs) are specified.

In this work, we employ the time-varying well BHPs as the control variable. We denote the (time-dependent) control vector as

[TABLE]

The time-varying BHP profiles constitute the input parameters, i.e.,

[TABLE]

thus, $N_{\mu}=N_{u}N_{t}$. Alternatively, well flow rates $(q_{j})_{d}$, $d\in\mathcal{D}$, or a combination of rates and BHPs, could be prescribed as the control variables.

Following  [42, 43, 45], the discretized set of nonlinear algebraic equations (obtained using fully implicit discretization111Note that this assumes that a linear multistep method with $k=1$ step is applied for time integration.) describing the HFM is represented as

[TABLE]

where $\textbf{g}:({\bar{\textbf{x}}^{1}},{\bar{\textbf{x}}^{2}},{\bar{\bf{u}}})\mapsto\textbf{g}({\bar{\textbf{x}}^{1}},{\bar{\textbf{x}}^{2}},{\bar{\bf{u}}})$ and $\textbf{g}:\mathbb{R}^{2N_{{c}}}\times\mathbb{R}^{2N_{{c}}}\times\mathbb{R}^{N_{u}}\rightarrow\mathbb{R}^{2N_{{c}}}$ designates the residual vector we seek to drive to zero, and superscript $n$ denotes the value of a variable at time step $n$. The state operator $g^{{n}}$ defined in Equation (1) is implicitly defined by the sequential solution to Equation (27).

3.2 Surrogate model: POD–TPWL

We now briefly describe the POD–TPWL formulation for oil–water systems, which will be our surrogate model for this application. For further details, the reader is referred to [42, 43, 45, 46, 47]. Given a set of ‘test’ controls ${\bf{u}}^{n}$, $n=1,\ldots,N_{t}$, the POD–TPWL model linearizes the residual around a previously saved ‘training’ simulation solution. Then, at time step $n$ in the test simulation, rather than solving the system of nonlinear algebraic equations (27) using, e.g., Newton’s method, we instead solve the system of linear algebraic equations

[TABLE]

where we have used the fact that ${\mathbf{g}}({\bf{{\acute{\textbf{x}}}}}^{i},{\bf{{\acute{\textbf{x}}}}}^{i-\rm{1}},{\bf{\acute{u}}}^{i})=\textbf{0}$ and defined

[TABLE]

Here, $\acute{\textbf{x}}$ denotes the (saved) training simulation state, $\acute{\mathbf{u}}$ denotes the training controls, and $i$ denotes the time step associated with the training state about which we linearize. Note that an accent $\acute{\cdot}$ indicates that the quantity has been saved during training simulations. While the criterion for determining the ‘closest’ training configuration $({\bf{{\acute{\textbf{x}}}}}^{i},{\bf{{\acute{\textbf{x}}}}}^{i-\rm{1}},{\bf{\acute{u}}}^{i})$ about which to linearize is application dependent, here we use pore volume injected (PVI) to determine the appropriate training solution. PVI quantifies the fraction of the system pore space that has been filled by injected fluid (water in our case), and as such corresponds to a dimensionless time. Thus, we seek to linearize around a solution that has progressed to the same PVI as the current test solution. See [45] for further discussion and details on the computation of PVI.

To reduce the computational cost associated with solving Equation (28), we approximate the state x in a low-dimensional affine subspace, using POD, as ${\textbf{x}}\approx{\mathbf{\Phi}}{\textbf{z}}+\bar{\textbf{x}}$, where $\mathbf{\Phi}\in\mathbb{R}^{2N_{{c}}\times\ell}$ denotes a POD basis, $\textbf{z}\in\mathbb{R}^{\ell}$ designates the reduced state, and $\bar{\textbf{x}}\in\mathbb{R}^{2N_{{c}}}$ indicates a reference state, which is often taken to be the mean of the snapshots. Replacing x with ${\mathbf{\Phi}}{\textbf{z}}+\bar{\textbf{x}}$ in Equation (28) yields

[TABLE]

where ${\bf{\acute{z}}}\vcentcolon=\mathbf{\Phi}^{T}(\acute{\textbf{x}}-\bar{\textbf{x}})$.

Because Equation (29) is overdetermined ($2N_{{c}}$ equations and $\ell<2N_{{c}}$ unknowns), it may not have a solution. Thus, we reduce the number of equations to $\ell$ by forcing the residual in Equation (29) to be orthogonal to the range of a test basis ${\mathbf{\Psi}}\in\mathbb{R}^{2N_{{c}}\times\ell}$. In line with previous studies on the application of POD–TPWL for subsurface flow models [45, 46, 47], we employ the least-squares Petrov–Galerkin (LSPG) test basis [44, 61, 35], i.e., ${\bf{\mathbf{\Psi}}}^{i}={\bf{\acute{J}}}^{i}{\mathbf{\Phi}}$. Premultiplying Equation (29) by $({\bf{\mathbf{\Psi}}}^{i})^{T}$, the linear system of equations in the low-dimensional space is now expressed as

[TABLE]

where

[TABLE]

and the subscript $\rm{RL}$ indicates that this is the POD–TPWL representation. Thus, the surrogate-state operator $g_{\rm{surr}}^{{n}}$ defined in Equation (3) is implicitly defined by the sequential solution to Equation (30) for this application. Further, this implies $\widetilde{\mathbf{P}}=\mathbf{P}\mathbf{\Phi}$. We also define the prolongation operator associated with the output in well block $d$ as $\widetilde{\mathbf{P}}_{d}\vcentcolon=\mathbf{P}_{d}\mathbf{\Phi}$ for this system. During online POD–TPWL computations, well-block saturations computed as $\widetilde{\mathbf{P}}_{d}{\bf{z}}^{n}$ can fall outside of the physical range. In this case, any $S<0$ is mapped to $S=0$, and $S>1$ is mapped to $S=1$.

POD–TPWL requires some number of training runs to be performed during an offline stage. These training simulations involve solving Equation (27) for prescribed time-varying BHPs $\boldsymbol{\mu}\in\mathcal{T}_{\text{TPWL}}\subseteq\mathbb{R}^{N_{\mu}}$, where $\mathcal{T}_{\text{TPWL}}$ denotes the TPWL training points. The state snapshots generated from the training simulations are saved and used to construct the POD basis $\mathbf{\Phi}$ by performing SVD on the (centered) snapshots. The constructed POD basis is then used to perform offline processing, which involves computing and saving quantities such as ${\bf{\acute{J}}}^{i}_{r},{\bf{\acute{B}}}^{i}_{r}$ and ${\bf{\acute{C}}}^{i}_{r}$. Although several training runs are used to construct the POD basis, only one of these, referred to as the primary training run, is used for linearization (i.e., ${\bf{\acute{J}}}^{i}$, ${\bf{\acute{B}}}^{i}$ and ${\bf{\acute{C}}}^{i}$ all derive from the primary training run).

3.3 Error modeling

To model QoI errors we adopt two approaches. Approach 1 is a hybrid treatment wherein we apply Method 4 (Section 2.3.4) to model relative errors in the well-block pressure and Method 3 (Section 2.3.3) to model errors in the well-block saturation. We pursue this approach because the error in the well-block pressure can exhibit a wide range of values, which makes modeling the relative error an easier task. On the other hand, the error in the well-block saturation spans a narrow range of absolute values because $0\leq S\leq 1$; thus, directly modeling the state error is appropriate. Approach 2 simply applies Method 2 (Section 2.3.2) to model the relative QoI error directly.

3.4 Feature design

The definition of the QoI in well block $d$ (24), the HFM governing equations (27), and the surrogate-model governing equations (30) suggest that the corresponding sampled-state error $\boldsymbol{\delta}_{ss}^{{n}}$ and QoI error $\delta_{q}^{{n}}$ will likely depend on data such as the states/controls about which POD–TPWL is linearized $(\widetilde{\mathbf{P}}_{d}{\bf{\acute{z}}}^{i},\widetilde{\mathbf{P}}_{d}{\bf{{\acute{z}}}}^{i-\rm{1}},{\bf{\acute{u}}}^{i})$, the current well-block state $\widetilde{\mathbf{P}}_{d}{\bf{z}}^{n}$, the previous well-block state $\widetilde{\mathbf{P}}_{d}{\bf{z}}^{n\rm{-1}}$, and operators associated with the linearized system such as $\mathbf{P}_{d}{\bf{\acute{J}}}^{i}\mathbf{P}_{d}^{T}\in\mathbb{R}^{2\times 2}$, $\mathbf{P}_{d}{\bf{\acute{B}}}^{i}\mathbf{P}_{d}^{T}\in\mathbb{R}^{2\times 2}$, and $\mathbf{P}_{d}{\bf{\acute{C}}}^{i}\mathbf{P}_{d}^{T}\in\mathbb{R}^{2\times 2}$, where the $2\times 2$ matrix is converted to a $1\times 4$ row vector. Similarly, these errors may also depend on other quantities such as the PVI in the primary training (${\rm{PVI}}^{i}$) and the test case (${\rm{PVI}}^{{n}}$), as well as data such as time-step sizes $\Delta t^{{n}}$ and $\Delta\acute{t}^{i}$, and time instances $t^{{n}}$ and $\acute{t}^{i}$. Some of the features used in this work are shown in Table 1. We note that most of these features are already computed during the course of the ROM simulation (e.g., the control input $u^{{n}}_{d}$), while others (e.g., $\dfrac{\langle{\bf{z}}^{n},{\bf{{\acute{z}}}}^{i}\rangle}{{\lVert{\bf{z}}^{n}\rVert}_{2}~{}{\lVert{\bf{{\acute{z}}}}^{i}\rVert}_{2}}$) can be computed via inexpensive computations, with cost scaling with the ROM dimension $\ell$. Further, while Hessian information is useful for informing the linearization error of TPWL, we do not include it in the feature set due to the high cost of its computation [47].

Note that solving Equation (30) for time step ${n}$ requires information from time step ${n}-1$; however, we could also include data from multiple previous time steps in the set of features. To include features generated over a ‘memory’ of $\tau$ previous time steps, we define a new feature vector ${{\boldsymbol{{f}}_{\rm{mem}}^{{n}}}}(\tau)\vcentcolon=\left[{\boldsymbol{{f}}^{{n}}}\ \boldsymbol{{f}}^{{n}-1}\ \cdots\ \boldsymbol{{f}}^{{n}-\tau}\right]\in\mathbb{R}^{1\times\left(1+\tau\right){N_{f}}}$. Note that some features in ${{\boldsymbol{{f}}_{\rm{mem}}^{{n}}}}(\tau)$ will be strongly correlated (some will be identical); we remove such highly correlated features in a preprocessing step by computing the feature–feature Pearson correlation coefficients for all pairs. Further details on correlation-criteria-based feature selection can be found in [62].

3.5 EMML training and test data

Following [47], we generate a set of $N_{\text{BHP}}$ BHP controls by adding unique random perturbations to the ‘primary’ training BHP control $\acute{\boldsymbol{\mu}}\equiv\left[(\acute{\textbf{u}}^{1})^{T}\ \cdots\ (\acute{\textbf{u}}^{N_{t}})^{T}\right]^{T}\in\mathcal{T}_{\text{TPWL}}$ used in the primary training run. For a given set of BHP controls, we define the perturbation in producer BHPs as

[TABLE]

The perturbation in injection BHPs $\Delta u^{I}(\boldsymbol{\mu})$ is defined by instead performing summation over the set of injector wells $\mathcal{D}_{I}$.

We partition the $N_{\text{BHP}}$ BHP schedules into $N_{\text{train}}<N_{\text{BHP}}$ clusters according to their representation in the (two-dimensional) space defined by $\Delta u^{P}$ and $\Delta u^{I}$. We select the BHP schedules closest to the cluster centers as ‘representative’ schedules for which we simulate both the POD–TPWL and the high-fidelity models. This subset of controls $\mathcal{T}_{\text{EMML}}\vcentcolon=\{\boldsymbol{\mu}_{\text{train}}^{1},\ldots,\boldsymbol{\mu}_{\text{train}}^{N_{\text{train}}}\}\subseteq\mathbb{R}^{N_{\mu}}$ constitutes the EMML training data. The remaining set of $N_{\text{BHP}}-N_{\text{train}}$ test schedules comprises the EMML test data $\mathcal{S}_{\text{EMML}}\vcentcolon=\{\boldsymbol{\mu}_{\text{test}}^{1},\ldots,\boldsymbol{\mu}_{\text{test}}^{N_{\text{BHP}}-N_{\text{train}}}\}\subseteq\mathbb{R}^{N_{\mu}}$. We note that the test set is not used in the construction of the EMML model; it is simply used to assess EMML performance after the model has been constructed.

3.6 Regression-model locality for POD–TPWL

As described in Section 2.5, we determine regression-model locality using classification and clustering to construct tailored local regression models for error prediction. Recall that QoI, described in Equation (24), correspond to the oil and water flow rates at the production wells $\mathcal{D}_{P}$ and water flow rates at the injection wells $\mathcal{D}_{I}$. We construct a local regression model only for production-well QoI $(q_{j})_{d}$, $j=o,w$, $d\in\mathcal{D}_{P}$; in fact, we determine regression-model locality for each production well independently, which is valid for both the oil- and water-flow rate QoI at that well. We do not construct local regression models for injection-well QoI $(q_{w})_{d}$, $d\in\mathcal{D}_{I}$, as a global regression model performs well in this case due to the relatively simple behavior of the associated QoI errors.

3.6.1 Classification

For production-well QoI, we partition the EMML training data into four categories. This partitioning is based on well-block saturation $S_{d}$, $d\in\mathcal{D}_{P}$, as shown in Figure 1, where the blue curve represents the POD–TPWL prediction and the black curve the corresponding HFM prediction.

The four categories correspond to different stages of the system and are referred to as $A,B^{+},B^{-},$ and $C$. All samples with $(S_{d})_{\rm{RL}}\leq\epsilon_{A}$ and $S_{d}\leq\epsilon_{A}$ are assigned to category $A$, where $(S_{d})_{\rm{RL}}$ denotes the well-block saturation predicted by POD–TPWL and, as before, $S_{d}$ denotes well-block saturation from the high-fidelity simulation. Thus, all the samples in category $A$ have close agreement between the POD–TPWL and HFM solutions; this category corresponds to the state before water ‘breakthrough’ occurs at a particular production well. Samples with $\epsilon_{A}<(S_{d})_{\rm{RL}}\leq S_{d}\leq\epsilon_{C}$ are assigned to category $B^{+}$, while samples with $\epsilon_{A}\leq S_{d}<(S_{d})_{\rm{RL}}\leq\epsilon_{C}$ are assigned to category $B^{-}$. Finally, samples with $S_{d}>\epsilon_{C}$ are assigned to category $C$. This category corresponds to significant water production. The actual values used in this work for $\epsilon_{A}$ and $\epsilon_{C}$ are given in Section 4.

As mentioned in Section 2.5, we perform classification using classification features $\boldsymbol{{f}}_{\rm{c}}\in\mathbb{R}^{1\times N_{f_{c}}}$. These features quantify (1) the perturbation in the prescribed control variables $u^{k}_{d},\ k=1,\ldots,n,~{}d\in\mathcal{D}_{P}$ relative to the primary training BHP controls $\acute{u}^{k}_{d},\ k=1,\ldots,n,~{}d\in\mathcal{D}_{P}$, and (2) the differences in well-block pressure for producer–injector pairs, which in turn may impact the velocity field as indicated by Equation (22b). For a production well located in grid block $d\in\mathcal{D}_{P}$, classification features include quantities such as the difference between the test BHP schedule and the primary training run BHP controls $\acute{\boldsymbol{\mu}}$, i.e., \Big{(}\sum\limits_{k=1}^{n}\left({u}^{k}_{d}-{\acute{u}}^{k}_{d}\right)^{2}\Big{)}^{\nicefrac{{1}}{{2}}}, the average well-block pressure difference between all producer-injector pairs, $\frac{1}{n}\sum\limits_{k=1}^{n}({\textbf{e}}_{2d-1}^{T}-{\textbf{e}}_{2d^{\prime}-1}^{T})\mathbf{\Phi}\textbf{z}^{k}$, $d^{\prime}\in\mathcal{D}_{I}$, and the average well-block pressure difference between the test case and the primary training simulation, represented by $\frac{1}{n}\sum\limits_{k=1}^{n}({\textbf{e}}_{2d-1}^{T}\mathbf{\Phi})(\textbf{z}^{k}-\acute{\textbf{z}}^{k})$. Table 2 reports some of the classification features employed in the current application.

4 Numerical results

In this section, we present numerical results for the application described in Section 3. The specific problem involves flow simulation in a synthetic two-dimensional horizontal reservoir. The reservoir model contains $50\times 50$ grid blocks such that $N_{{c}}=2500$ and $N_{x}=5000$. It contains three production wells $|\mathcal{D}_{P}|=3$, which we label as $P_{1}$, $P_{2}$, and $P_{3}$, and three injection wells $|\mathcal{D}_{I}|=3$, which we label as $I_{1}$, $I_{2}$, and $I_{3}$. The six wells ($N_{u}=6$) are shown in Figure 2. The permeability field is isotropic, i.e., $\textbf{k}=\mathrm{diag}(k)$, and the porosity is set to $\phi=0.2$. The relative permeability functions are prescribed to be $k_{rw}(S)=S^{2}$ and $k_{ro}(S)=(1-S)^{2}$. We apply a backward Euler time integrator with adaptive time-step selection.

Three training simulations, $\lvert\mathcal{T}_{\text{TPWL}}\rvert=3$, are performed to construct the POD–TPWL model (the three runs provide a sufficient number of snapshots for the POD basis), from which $\ell=150$ POD basis vectors are extracted. Of these, 90 correspond to the saturation state variables and 60 to the pressure state variables. Figure 3 depicts the BHP controls $\acute{\boldsymbol{\mu}}\in\mathcal{T}_{\text{TPWL}}$ applied in the primary training simulation (recall that this is the run used for linearization). These time-varying BHPs, as well as those considered in the test runs, are meant to be representative of the BHP schedules that can arise during oil production optimization computations. In such optimizations, the goal is to determine the time-varying BHPs that maximize an economic metric, or the cumulative oil recovered from the reservoir.

We consider $N_{\text{BHP}}=200$ sets of BHP controls to construct the EMML training $\mathcal{T}_{\text{EMML}}$ and test sets $\mathcal{S}_{\text{EMML}}$. As described in Section 3.5, each of these sets is characterized by time-varying BHPs ${\bf{u}}^{n}$, $n=1,\ldots,N_{t}$, obtained by adding a unique (time-varying) random perturbation to the primary training BHP profiles. The time-varying BHPs for a particular case (Case 1) are shown in Figure 4. The frequency of change in the primary training BHP schedule is every 200 days (Figure 3), while the frequency of change in the BHP schedule for Case 1 is every 175 days (Figure 4).

We consider a total of ${N_{f}}=84$ features, which include those listed in Table 1. After neglecting highly correlated features, the total number of retained features is reduced to around 40 (note that each QoI may retain a different subset of features). If we consider a memory of $\tau={1}$, we obtain ${N_{f}}=84\times 2=168$ features; this is reduced to about 64 after neglecting highly correlated features.

Based on extensive numerical experiments, we observed that the highest EMML accuracy was obtained on the test set $\mathcal{S}_{\text{EMML}}$ using Approach 1 (i.e., error-modeling Method 4 in Section 2.3.4 for the well-block pressure, and Method 3 in Section 2.3.3 for the well-block saturation), a memory of $\tau={1}$, $N_{\text{train}}={30}$ EMML training points, classification for determining regression-model locality for production-well QoIs, and random forests (RF) for regression. While performing classification, we set $\epsilon_{A}=0.05$ and $\epsilon_{C}=0.6$. These values are somewhat heuristic, but they appropriately identify the basic behaviors (solution stages) we wish to capture through classification. We compute the hyper-parameters for random forests by minimizing the out-of-bag error. We first report the numerical results corresponding to these best-case parameters. Then, in Section 4.3, we quantify EMML performance for other choices of algorithmic parameters ; for example, we assess the effect of employing clustering to determine regression-model locality, using only global regression models, and applying LASSO for regression. From Section 2.7 we recall that there are two possible applications of the EMML QoI-error prediction: (1) as a correction to the surrogate-model QoI, or (2) as an error indicator to be used within the ROMES framework. Here we consider the first application, and in Section 4.4 the second.

4.1 EMML for QoI correction: Test Case 1

We first present results for Test Case $1$, represented by $\boldsymbol{\mu}_{\text{test}}^{1}\in\mathcal{S}_{\text{EMML}}$. As will be described in Section 4.2, this case corresponds to the median time-integrated POD–TPWL error in the test set $\mathcal{S}_{\text{EMML}}$. Figure 5 reports results for the pressure for wells $P_{1}$ and $I_{3}$ and the saturation for well $P_{1}$; note that these quantities associate with the sampled state that is corrected as an intermediate step in Methods 3 and 4. Figure 6 presents several QoI; these correspond to the oil and water production rates in well $P_{1}$, and the water injection rate for well $I_{3}$. We focus on wells $P_{1}$ and $I_{3}$, as they are the wells with the highest cumulative liquid production and injection.

For Test Case $1$, the POD–TPWL prediction (blue curve) for the saturation for well $P_{1}$ has an error that is most noticeable at around $500$ days in Figure 5b. Similarly, the POD–TPWL error in the production rates is evident at around the same time in Figure 6a,b. In Figure 5b, we observe that the EMML-corrected well-block saturation (red curve) demonstrates improved accuracy around the time of water breakthrough ($\sim 410$ days). Similarly, in Figure 6, we see that the EMML-corrected flow rates display better accuracy than the POD–TPWL results. The improvement is most apparent in the breakthrough prediction in Figure 6b, and in oil production rate (Figure 6a) at a time of about 500 days.

4.2 EMML for QoI correction: additional test cases

We now present results for two additional test cases with control vectors $\boldsymbol{\mu}\in\mathcal{S}_{\text{EMML}}$, which correspond to different POD–TPWL prediction errors. We then assess EMML performance for an ensemble containing the entire EMML test set $\mathcal{S}_{\text{EMML}}$ ($|\mathcal{S}_{\text{EMML}}|=170$ cases).

BHP schedules for Test Cases 2 and 3—represented by control vectors $\boldsymbol{\mu}_{\text{test}}^{2}$ and $\boldsymbol{\mu}_{\text{test}}^{3}$, respectively—are shown in Figures 7 and 9. Test Case 2, for which $\Delta u^{P}(\boldsymbol{\mu}_{\text{test}}^{2})=0.54$ and $\Delta u^{I}(\boldsymbol{\mu}_{\text{test}}^{2})=0.10$, corresponds to a smaller perturbation in the BHPs relative to the primary training run BHPs compared to that in Test Case 1 ($\Delta u^{P}(\boldsymbol{\mu}_{\text{test}}^{1})=0.68,~{}\Delta u^{I}(\boldsymbol{\mu}_{\text{test}}^{1})=0.23$). It also corresponds to lower POD–TPWL error compared to Test Case 1. The results for production and injection rates for Test Case 2 are shown in Figure 8. The POD–TPWL error is again most noticeable at around 500 days for the oil and water production rates at well $P_{1}$. The correction is clearly evident in Figure 8a at around 500 days and in Figure 8b at around 500 and 1250 days. Slight improvement in water injection rate (Figure 8c) is also apparent just before 200 days.

Test Case 3, with $\Delta u^{P}(\boldsymbol{\mu}_{\text{test}}^{3})=0.46$ and $\Delta u^{I}(\boldsymbol{\mu}_{\text{test}}^{3})=0.32$, corresponds to a higher perturbation in the injector BHPs compared to that in Test Case 1, and it leads to a larger POD–TPWL error. The POD–TPWL error is again most evident at around 500 days for both the oil and water production rates as shown in Figure 10a,b. These results are again significantly improved by the proposed EMML-based correction. We note finally that the corrected solutions for the production and injection rates display fluctuations at some times. This is because, when constructing the corrections, we treat each time step as independent, consistent with the i.i.d. assumption.

To quantify EMML performance over the entire test set $\mathcal{S}_{\text{EMML}}$ of $|\mathcal{S}_{\text{EMML}}|=170$ cases, we define the following relative time-integrated error measures for the POD–TPWL and corrected solutions:

[TABLE]

Here, $E_{\rm{RL}}$ denotes the relative average time-integrated error in the POD–TPWL solution and $E_{\rm{corr}}$ designates the average time-integrated error in the EMML-corrected solution. Note that $\mkern 1.5mu\overline{\mkern-1.5mu\mathcal{D}\mkern-1.5mu}\mkern 1.5mu=\mathcal{D}_{P}$ (with $j=o$ or $j=w$) for production wells, while $\mkern 1.5mu\overline{\mkern-1.5mu\mathcal{D}\mkern-1.5mu}\mkern 1.5mu=\mathcal{D}_{I}$ (with $j=w$) for injection wells.

Figure 11 displays the time-integrated errors for the entire set of 200 cases in the EMML training $\mathcal{T}_{\text{EMML}}$ and test $\mathcal{S}_{\text{EMML}}$ sets. The cases are sorted by increasing POD–TPWL error. For each case in the ensemble, the figure reports the time-integrated error in the POD–TPWL prediction (blue) and the EMML-corrected POD–TWPL predictions, which are further distinguished by whether they correspond to cases in the EMML training set $\mathcal{T}_{\text{EMML}}$ (green) or test set $\mathcal{S}_{\text{EMML}}$ (red). Consistent with the results presented in Section 4.1, the time-integrated errors for all test cases are reduced after application of the EMML-driven correction. Test Cases 1, 2 and 3, discussed above, are identified in Figure 11. These three test cases correspond to the $10$th, $50$th and $90$th percentiles in $E_{\rm{RL}}(o,\mathcal{D}_{P})$. Note that the time-integrated error in the EMML training data (green points) is small but nonzero, which indicates that the random forest model $\hat{r}_{ss,i}$ does not perfectly fit the data. This is intentional, as it prevents overfitting, which can potentially lead to large errors in EMML test-case predictions.

Table 3 presents the median errors for the test set $\mathcal{S}_{\text{EMML}}$ results displayed in Figure 11. We observe that by applying the EMML correction, we reduce the three errors, $E_{\rm{RL}}(o,\mathcal{D}_{P})$, $E_{\rm{RL}}(w,\mathcal{D}_{P})$ and $E_{\rm{RL}}(w,\mathcal{D}_{I})$, by about $38\%$ on average. Although we achieve substantial improvements at certain time instances using EMML (as is evident in Figure 6a,b), small errors in time persist. The EMML procedure reduces these errors but it does not completely eliminate them.

We note that one source of error in the EMML predictions is misclassification, which in turn leads to using the local regression model from the incorrect category. The average misclassification error—defined as the ratio of the number of EMML test samples misclassified to the total number of EMML test samples—over cases in the EMML test set $\mathcal{S}_{\text{EMML}}$ is 3% in this set of experiments. The misclassification error is primarily from misclassifying samples whose actual category is $B^{+}$ as $B^{-}$, and vice-versa.

Finally, from the variable importance plots generated by the random forest regression model (see  [64] for more details on variable importance plots), we observe some of the key features in this application to be $\widetilde{\mathbf{P}}_{2d}{\bf{z}}^{n}-\widetilde{\mathbf{P}}_{2d}{\bf{{\acute{z}}}}^{i}$, $\widetilde{\mathbf{P}}_{2d-1}{\bf{z}}^{n}-\widetilde{\mathbf{P}}_{2d-1}{\bf{{\acute{z}}}}^{i}$, $\widetilde{\mathbf{P}}_{2d}{\bf{z}}^{n\rm{-1}}-\widetilde{\mathbf{P}}_{2d}{\bf{{\acute{z}}}}^{i-\rm{1}}$, $\widetilde{\mathbf{P}}_{2d}{\bf{z}}^{n}$, $\widetilde{\mathbf{P}}_{2d}{\bf{{\acute{z}}}}^{i}$, ${({\textbf{e}}_{2d}^{T}\mathbf{\Phi})({\bf{z}}^{n}-{\bf{z}}^{n\rm{-1}})}~{}/~{}{\Delta t^{{n}}}$, $u^{{n}}_{d}$, ${\rm{PVI}}^{{n}}$, and $\dfrac{\langle{\bf{z}}^{n},{\bf{{\acute{z}}}}^{i}\rangle}{{\lVert{\bf{z}}^{n}\rVert}_{2}~{}{\lVert{\bf{{\acute{z}}}}^{i}\rVert}_{2}}$.

4.3 EMML for QoI correction: alternative EMML parameters

For completeness, we now analyze EMML performance using different algorithmic parameters. Table 4 reports the median time-integrated errors $E_{\rm{corr}}(j,\mkern 1.5mu\overline{\mkern-1.5mu\mathcal{D}\mkern-1.5mu}\mkern 1.5mu)$ over the 170 test cases corresponding to $\boldsymbol{\mu}\in\mathcal{S}_{\text{EMML}}$. The EMML parameters used here differ from the best-case parameters employed in Sections 4.1–4.2. We vary the memory $\tau$, the number of high-fidelity simulations $N_{\text{train}}$ used to construct the EMML training data, and the method for determining regression-model locality clustering or classification). For a detailed discussion of these results, we refer the reader to [65].

We observe from Table 4 that the EMML-based corrections obtained using Approach 1 lead to more accurate results than those obtained by Approach 2 (these two approaches are defined in Section 3.3). Additionally, a decrease in the number of high-fidelity simulations ($N_{\text{train}}$) used to build the EMML training dataset leads to a noticeable decrease in accuracy. Reduced accuracy is also observed when LASSO regression is used instead of random-forest regression. In addition, using classification (a supervised machine learning technique) to determine regression-model locality prior to constructing local RF models performs better than the use of clustering (an automated unsupervised learning approach). Even in the absence of employing local regression models (which adds complexity to the EMML method), EMML with a global error model still provides improved accuracy relative to POD–TPWL, though more accurate results are achieved when classification is used to determine locality. We note finally that the impact of memory is very small for this test set.

4.4 EMML as an error indicator for ROMES

We now consider the second application of the EMML error model: as an error indicator for the ROMES method [19]. In particular, we apply the framework proposed in Section 2.7.2 with the following scalar-valued function of the surrogate QoI errors:

[TABLE]

We also denote the average value of ${h}$ associated with QoI errors $\delta_{q}$ for a given phase $j$ and a specified set of wells $\bar{\mathcal{D}}$ as

[TABLE]

Figure 12 displays cross-plots of the true time-integrated error ($\bar{h}(\delta_{q};o,\mathcal{D}_{P})$, $\bar{h}(\delta_{q};w,\mathcal{D}_{P})$, and $\bar{h}(\delta_{q};w,\mathcal{D}_{I})$) versus the EMML- approximated time-integrated error ($\bar{h}(\hat{\delta}_{q};o,\mathcal{D}_{P})$, $\bar{h}(\hat{\delta}_{q};w,\mathcal{D}_{P})$, and $\bar{h}(\hat{\delta}_{q};w,\mathcal{D}_{I})$). Although scatter is apparent, the relationship is essentially linear. However, note that $\bar{h}(\delta_{q};j,\bar{\mathcal{D}})$ is generally greater than $\bar{h}(\hat{\delta}_{q};j,\bar{\mathcal{D}})$, which leads to a systematic bias. Thus, applying the EMML-computed quantities $\bar{h}(\hat{\delta}_{q};o,\mathcal{D}_{P})$, $\bar{h}(\hat{\delta}_{q};w,\mathcal{D}_{P})$, and $\bar{h}(\hat{\delta}_{q};w,\mathcal{D}_{I})$ to predict their ‘true’ counterparts will be biased; this is reflected by the red line in Figure 12, which corresponds to the prediction if the EMML-computed quantity alone is applied for prediction. Note that this bias is not trivial to fix within the regression method itself, as the regression was constructed for predicting time-instantaneous errors, while the observed bias is present for time- and well-averaged errors.

As described in Section 2.7.2, we can address this issue using the ROMES method [19]. This technique applies Gaussian-process (GP) regression to model the (generally unknown) average true error $\bar{h}(\delta_{q};j,\bar{\mathcal{D}})$ using an error indicator, which (in this case) corresponds to the (computable) average error predicted by EMML $\bar{h}(\hat{\delta}_{q};j,\bar{\mathcal{D}})$. To construct the GP, we employ 15 additional high-fidelity and POD–TPWL simulations for parameter instances corresponding to $\mathcal{T}_{\text{ROMES}}\subset\mathcal{S}_{\text{EMML}}$, i.e., $\lvert\mathcal{T}_{\text{ROMES}}\rvert=15$. These simulations provide the ROMES training data $\{(\bar{h}(\delta_{q}(\boldsymbol{\mu});j,\bar{\mathcal{D}}),\bar{h}(\hat{\delta}_{q}(\boldsymbol{\mu});j,\bar{\mathcal{D}}))\}_{\boldsymbol{\mu}\in\mathcal{T}_{\text{ROMES}}}.$ We then construct a GP using the DACE package [66] with a first-order-polynomial mean function, and a Gaussian covariance function.

Figure 13 shows the resulting GP. This figure displays the Gaussian process with a 90% confidence interval (shaded region), the true average error (black), the average error predicted by EMML alone (red points), and the average error predicted by EMML after post-processing with ROMES (purple points). Most importantly, these results show that the (mean) EMML prediction after ROMES postprocessing is significantly more accurate (i.e., closer to the true errors) than the EMML prediction alone. Further, the true time-integrated error for the majority of the test cases lies within the 90% confidence interval predicted by the GP; this demonstrates the importance of a statistical prediction rather than a deterministic prediction, as the confidence interval quantifies the prediction uncertainty. We thus conclude that the method proposed in Section 2.7.2 is effective at modeling time-integrated errors in this application.

4.5 Computational costs

We first discuss the computational costs incurred by POD–TPWL and EMML. All reported timings were obtained on a machine with dual E5520 Intel CPUs (4 cores, 2.26 Ghz) and 24 GB memory using a Matlab implementation of the high-fidelity and surrogate models and an $R$ [67] implementation of EMML. The offline computational costs for POD–TPWL entail (1) executing $|\mathcal{T}_{\text{TPWL}}|=3$ high-fidelity simulations—which can be done in parallel—for $\boldsymbol{\mu}\in\mathcal{T}_{\text{TPWL}}$; the cost of a single high-fidelity simulation is 370 seconds, and (2) assembling POD–TPWL operators via Equation (31), which consumes 23.5 seconds.

The offline computational costs for EMML training entail (1) executing $|\mathcal{T}_{\text{EMML}}|=N_{\text{train}}=30$ high-fidelity and POD–TPWL simulations (in parallel) for $\boldsymbol{\mu}\in\mathcal{T}_{\text{EMML}}$; the POD–TPWL solution takes only 0.5 seconds, which constitutes a speedup of approximately 700 relative to the HFM, (2) constructing a classification model for each of the $|\mathcal{D}_{P}|=3$ producer wells $d\in\mathcal{D}_{P}$; this consumes 1010 seconds per producer well, and (3) constructing a random-forest regression model for each of the 9 QoI defined in Equation (24); this consumes 412 seconds per regression model. The total offline cost of steps 2 and 3—assuming serial computation—is 6735 seconds. We note that this is approximately 18 times costlier than a single high-fidelity simulation. Elapsed timings can be readily reduced through use of parallel processing (each QoI can be treated by a different processor). The ratio of offline costs relative to HFM simulation cost will be smaller for larger-dimensional and more complex HFMs, as the EMML training cost in steps 2 and 3 is independent of the complexity of the HFM. However, the EMML training cost does scale linearly with the number of QoI, assuming serial processing.

In terms of online costs, the POD–TPWL simulation consumes only $0.5$ seconds (as mentioned above), and the online EMML error prediction takes about $3\times 10^{-4}$ seconds for querying the regression model. Thus, online costs are very small relative to offline costs. We note finally that production optimization computations in this setting may require $O(100-1000)$ flow simulations, so an offline cost of $O(10)$ HFM simulations, as required by the EMML-based framework, represents an acceptable overhead.

5 Concluding remarks

In this work we introduced a general method for error modeling using machine learning (EMML). We applied the EMML framework for modeling the error introduced by surrogate models of dynamical systems. The framework employs high-dimensional regression methods from machine learning to map a set of inexpensively computed error indicators (features) to a prediction of the (time-dependent) surrogate-model quantity-of-interest error (response). The method requires first constructing an EMML training dataset by simulating both the surrogate model and the high-fidelity model for some instances of the input parameters. In particular, we proposed:

• •

four different methods for modeling the error (Section 2.3),

• •

two methods (classification and clustering) for determining the notion of locality that is employed to construct local regression models (Section 2.5),

• •

two techniques (random forests and LASSO) for performing regression (Section 2.6), and

• •

two applications of the resulting error models: as a correction to the surrogate-model QoI prediction, and as a way to model functions of the QoI error using the ROMES method (Section 2.7).

We specialized the method to one particular application: subsurface-flow modeling with a POD–TPWL reduced-order model as a surrogate. For this application we proposed specific EMML method ingredients, such as particular choices for error modeling (Section 3.3), feature design (Section 3.4), training and test data (Section 3.5), and classification features to use for determining regression-model locality (Section 3.6).

In the numerical experiments, we observed that the EMML method performed the best using the following algorithmic parameters: error-modeling using Approach 1, memory $\tau={1}$, $N_{\text{train}}={30}$ simulations to train the EMML model, classification to determine regression-model locality, and random-forest regression. When the EMML error models are used as a correction to the surrogate-model prediction, we demonstrated improved accuracy in the output quantities of interest relative to the original POD–TPWL surrogate model for a large number of test cases (Sections 4.1–4.2). When the EMML error models are used to model functions of the QoI error via ROMES, we observed that the EMML prediction—when combined with a ROMES-based Gaussian-process model—produced an accurate prediction with statistical confidence intervals. It is important to note, however, that the EMML offline cost is not negligible (Section 4.5), as it entails (1) constructing the EMML training dataset, which requires executing $N_{\text{train}}$ high-fidelity and surrogate-model simulations, (2) determining regression-model locality for every QoI, and (3) constructing a regression model for every QoI. Thus, this framework is only appropriate for use in many-query problems such as optimization and uncertainty quantification.

The EMML framework provides a general error modeling methodology for surrogate models of dynamical systems; it assumes only that the surrogate produces a large set of features that can be mined for potential error indicators. Thus, it is applicable to a wide variety of surrogates, such as reduced-order models (considered in this work) and those presented in [61, 68, 69, 70, 71], and coarsened models, for example. The application of EMML with upscaled (effectivized) subsurface flow models, within the context of uncertainty quantification, was considered by [65]. Future work should be directed toward modeling error with other types of surrogate models (for a range of applications), and also for modeling errors introduced while performing geological parameterization [72, 73]. Algorithmic improvements within the EMML framework could also be considered. For example, our approach entails univariate regression for each QoI considered; future work could investigate the use of multivariate regression that accounts for interactions between the QoI. It may also be worthwhile to explore other regression methods, such as artificial neural networks with long short-term memory [74]. However, our initial (and ongoing) investigations into the use of deep artificial neural networks for constructing error models demonstrate that the required training time is very large relative to that incurred by high-fidelity simulations, so other approaches should also be considered.

6 Acknowledgments

We thank the industrial affiliates of the Stanford University Smart Fields and Reservoir Simulation Research (SUPRI-B) Consortia for financial support. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

Appendix A Random-forest regression

Random forest is a supervised machine-learning technique based on constructing an ensemble of decision trees; it can be used for both classification and regression. Here, decision trees are constructed by segmenting the feature space along (canonical) directions corresponding to individual features. In Figure 14, which is adapted and modified from [5], we plot the (color-coded) true error—which is the response we aim to predict—as a function of two features corresponding to synthetic training data.

As shown in Figure 14a, segmentation of the feature space involves splitting the domain into regions $R_{1}$ and $R_{2}$. This is achieved by computing the feature index $j\in\{1,\ldots,{N_{f}}\}$ and cutpoint $s\in\mathbb{R}$ corresponding to the segmentation that minimizes the residual sum of squares

[TABLE]

Here, we define feature-space regions $R_{1}(j,s)=\{\boldsymbol{{f}}\ |\ f_{j}<s\}$ and $R_{2}(j,s)=\{\boldsymbol{{f}}\ |\ f_{j}\geq s\}$ and we denote the mean value of the response across the training samples in $R_{k}$ as $\hat{\delta}_{R_{k}}\in\mathbb{R},\ k=1,2$. This segmentation is carried out recursively, as shown in Figure 14b. Because recursive segmentation can be summarized in a tree structure as depicted in Figure 15, this type of approach is known as a decision-tree method. The values assigned to leaves of the tree correspond to the response predicted in the feature-space region associated with those leaves, which is the mean value of the response across training data in that region. Note that segmenting the feature space in this way enables nonlinear interactions among the features to be captured.

Unfortunately, decision-tree models often exhibit high variance and low bias, which can result in overfitting the training data [5]. To resolve this problem, random-forest regression constructs an ensemble of decision trees, and the prediction corresponds to the average prediction across all trees in the ensemble. While constructing a given decision tree, the random-forest technique considers only a randomly selected subset of the feature indices as candidates for performing a split. This randomization serves to decorrelate the trees, thereby reducing the variance incurred by averaging the predictions from different trees, which acts to improve prediction quality.

Additionally, random-forest regression grows trees on a bootstrap-sampled version of the training data. Bootstrapping (a resampling technique involving sampling with replacement) is illustrated in Figure 16a for a data set containing five samples. To construct the first decision tree, the random-forest technique samples, with replacement, five training samples from the data shown in Figure 16a. However, due to sampling with replacement, it is probable that some samples in the resulting data set are repeated, as shown in Figure 16b. Bootstrapping further assists in variance reduction without increasing the bias.

While random forests demonstrate improved accuracy relative to a single decision tree, they lead to a loss of interpretation. We note that the performance of random forests can be improved by various mechanisms, such as pruning or recursively dropping the least-important features. For additional details on random forests, see Ref. [64].

Appendix B LASSO regression

Least absolute shrinkage and selection operator (LASSO) is a regression method that fits a linear model to the feature-to-error mapping, i.e., $\delta^{n}(\boldsymbol{\mu})=\beta_{0}+\sum\limits_{j=1}^{{N_{f}}}f^{n}_{j}(\boldsymbol{\mu})\beta_{j}+\varepsilon$, $n=1,\ldots,N_{t}$ and $\beta_{j}\in\mathbb{R}$, $j=0,\ldots,{N_{f}}$, while reducing the number of nonzero coefficients. It does so by computing coefficients $\beta_{j}$, $j=0,\ldots,{N_{f}}$, that minimize an objective function composed of the sum of squared errors and an $L_{1}$-penalty on the coefficients

[TABLE]

where $\lambda\in\mathbb{R}$ is a penalty parameter, and larger values of $\lambda$ promote sparsity in the computed coefficients.

The value of $\lambda$ is typically chosen by $k$-fold cross validation, which is a process that randomly partitions the training data into $k$ subsamples of equal size. One subsample is withheld while a linear model is constructed using the remaining $k-1$ subsamples. The model is then tested on the withheld sample, and the prediction errors are recorded. To reduce variability, the process is repeated $k$ times, with a different subsample withheld each time. Finally, the prediction errors are averaged, and this average is referred to as the cross-validation error. To determine the optimal value of $\lambda$ in practice, we compute the cross-validation error on only a subset of candidate values for $\lambda$, and select the value yielding the smallest cross-validation error.