Gaussian mixture model fitting method for uncertainty quantification by conditioning to production data

Abstract

For most history matching problems, the posterior probability density function (PDF) may have multiple local maxima, and it is extremely challenging to quantify uncertainty of model parameters and production forecasts by conditioning to production data. In this paper, a novel method is proposed to improve the accuracy of Gaussian mixture model (GMM) approximation of the complex posterior PDF by adding more Gaussian components. Simulation results of all reservoir models generated during the history matching process, e.g., using the distributed Gauss-Newton (DGN) optimizer, are used as training data points for GMM fitting. The distance between the GMM approximation and the actual posterior PDF is estimated by summing up the errors calculated at all training data points. The distance is an analytical function of unknown GMM parameters such as covariance matrix and weighting factor for each Gaussian component. These unknown GMM parameters are determined by minimizing the distance function. A GMM is accepted if the distance is reasonably small. Otherwise, new Gaussian components will be added iteratively to further reduce the distance until convergence. Finally, high-quality conditional realizations are generated by sampling from each Gaussian component in the mixture, with the appropriate relative probability. The proposed method is first validated using nonlinear toy problems and then applied to a history-matching example. GMM generates better samples with a computational cost comparable to or less than other methods we tested. GMM samples yield production forecasts that match production data reasonably well in the history-matching period and are consistent with production data observed in the blind test period.

Appendices

Appendix A: Finite mixture models

A finite mixture model [25] with K components, p (x |Θ), is defined as

$$ p_{MM}\left( x\left| \mathrm{{\Theta}} \right. \right)=\sum\nolimits_{k\mathrm{=1}}^{K} {f^{(k)}p^{(k)}(x\mathrm{,}\theta^{(k)})} . $$

(1)

In Eq. 1, x is an n-dimensional random vector, p^(k)(x,𝜃^(k)) is the probability density function (PDF) with parameters 𝜃^(k) for the k th mixture component with 1 ≤k≤K, and f^(k) is the positive-valued mixture weight for the k th mixture component, representing the probability to generate a sample from this component. By definition, we have ${\sum }_{k\mathrm {=1}}^{K} f^{(k)} \mathrm {=1}$. The complete set of parameters to define a mixture model with K components is Θ = {f⁽¹⁾,f⁽²⁾,…,f^(K);𝜃⁽¹⁾,𝜃⁽²⁾,…,𝜃^(K)}.

If each of the K mixture components is a Gaussian with mean μ^(k) and covariance matrix C^(k), i.e.,

$$ \begin{array}{@{}rcl@{}} p^{(k)}\left( x, \theta^{(k)}\right)&=&\left( 2\pi \right)^{-n \left/2\right.}{\text{Det}}^{-1/2}\left( C^{(k)}\right)\\ &\times& \exp\left\{ -\frac{1}{2}\left[ x - \mu^{(k)} \right]^{T}\left[ C^{(k)} \right]^{-1}\left[ x - \mu^{(k)} \right] \right\}, \end{array} $$

(2)

then the finite mixture model defined in Eq. 1 becomes a Gaussian mixture model (GMM). In Eq. 2, 𝜃^(k) = {μ^(k),C^(k)}.

Given a set of N_T samples (or sets of training data points): $T=\left \{ x^{\mathrm {(1)}}\mathrm {,}x^{\mathrm {(2)}}\mathrm {,\mathellipsis ,}x^{(N_{T})} \right \}$, which are assumed to be independent and identically distributed (i.i.d.). For given parameters Θ, the conditional probability of component k for a given sample x⁽ⁱ⁾ is called “membership weight” and according to the Bayes theorem, it can be computed as

$$ w^{(i\mathrm{,}k)}=\frac{f^{(k)}p^{(k)}\left( x^{(i)}\mathrm{,}\theta ^{(k)}\right)}{\sum\nolimits_{j\mathrm{=1}}^{K} {f^{(j)}p^{(j)}\left( x^{(i)}\mathrm{,}\theta ^{(j)}\right)}} , \text{ for } \mathrm{1\le} k\mathrm{\le} K \text{ and } \mathrm{1\le} i\mathrm{\le} N_{T}. $$

(3)

In this expression, the f^(k) is the (“prior”) probability for component k and the likelihood function, i.e., the probability (density) of the sample x⁽ⁱ⁾ from this given component is p^(k)(x⁽ⁱ⁾,𝜃^(k)). The denominator in Eq. 3 is the probability density for sample x⁽ⁱ⁾, and hence, the probability density to generate the N_T i.i.d. samples (for given parameters Θ) is

$$ p\left( x^{\left( 1 \right)},x^{\left( 2 \right)},\mathellipsis, x^{\left( N_{T} \right)}\left|{\Theta}\right.\right)=\prod\nolimits_{i=1}^{N_{T}}\left[\sum\nolimits_{j=1}^{K} {f^{(j)}p^{(j)}\left( x^{(i)},\theta^{(j)}\right)}\right]. $$

(4)

The maximum likelihood estimate of the parameters Θ can now be found using an iterative expectation and maximization (EM) algorithm.

Appendix B: Analytical derivatives of GMM with respect to unknown GMM parameters

To reduce the computational cost for GMM fitting, we need to apply a gradient-based optimization algorithm to minimize the error function, which requires calculation of the analytical gradient of the error function defined in Eq. ?? and/or calculation of the sensitivity matrix of the normalized shift defined in Eq. ??. By the chain rule, we have

$$ \mathrm{\nabla} G\left( {\Theta} \right)\mathrm{=-2}\sum\nolimits_{j\mathrm{=1}}^{N_{T}} {E\left( x^{(j)}\mathrm{;}{\Theta} \right)\mathrm{\nabla} p_{\text{GMM}}^{\left( j \right)}\left( {\Theta} \right)} . $$

(5)

The GMM approximation at a training data point x^(j) is

$$ p_{\text{GMM}}^{\left( j \right)}\left( {\Theta} \right)=p_{\text{GMM}}\left( x^{(j)}\mathrm{,}{\Theta} \right)=\sum\nolimits_{l\mathrm{=1}}^{K} \left[ \beta_{G}^{\left( l \right)}p_{G}^{\left( l \right)}(x^{\left( j \right)}\mathrm{,}H^{(l)}) \right] , $$

(6)

$$ p_{G}^{\left( l \right)}\left( x^{\left( j \right)}\mathrm{,}H^{\left( l \right)} \right)=e^{-\frac{1}{2}\left\{ \left[ x^{\left( j \right)}-\mu^{(l)} \right]^{T}H^{(l)}\left[ x^{\left( j \right)}-\mu^{(l)} \right] \right\}}. $$

(7)

The derivatives of $p_{\text {GMM}}^{\left (j \right )}\left ({\Theta } \right )$ with respect to w^(l), $\beta _{G}^{\left (l \right )}$, and v^(l,k) (for l= 1, 2,…,K and k= 1, 2,…,m) are

$$ \frac{\partial p_{\text{GMM}}^{\left( j \right)}\left( {\Theta} \right)}{\partial \beta_{G}^{\left( l \right)}}=p_{G}^{\left( l \right)}(x^{\left( j \right)},H^{(l)}). $$

(8)

$$ \frac{\partial p_{\text{GMM}}^{\left( j \right)}\left( {\Theta} \right)}{\partial w^{(l)}}\mathrm{=-}\frac{1}{2}\beta_{G}^{\left( l \right)}p_{G}^{\left( l \right)}(x^{\left( j \right)},H^{(l)})\left[ x^{\left( j \right)}\! -\!\mu^{(l)} \right]^{T}H_{G0}^{(l)}\left[ x^{\left( j \right)}-\mu^{(l)} \right]. $$

(9)

$$ \mathrm{\nabla}_{v^{(l\mathrm{,}k)}}p_{\text{GMM}}^{\left( j \right)}\left( {\Theta} \right)\mathrm{=-}\beta_{G}^{\left( l \right)}p_{G}^{\left( l \right)}(x^{\left( j \right)},H^{(l)})\left[ x^{\left( j \right)}\! -\!\mu^{(l)} \right]^{T}v^{(l\mathrm{,}k)}\left[ x^{\left( j \right)}\! -\! \mu^{(l)} \right]. $$

(10)