Abstract
We compare numerically the performance of reversible and non-reversible Markov Chain Monte Carlo algorithms for high-dimensional oil reservoir problems; because of the nature of the problem at hand, the target measures from which we sample are supported on bounded domains. We compare two strategies to deal with bounded domains, namely reflecting proposals off the boundary and rejecting them when they fall outside of the domain. We observe that for complex high-dimensional problems, reflection mechanisms outperform rejection approaches and that the advantage of introducing non-reversibility in the Markov Chain employed for sampling is more and more visible as the dimension of the parameter space increases.
Article PDF
Similar content being viewed by others
References
Eclipse – Industry Reference Reservoir Simulator, Reference Manual. Version 2015.1
Bouchard-Côte, A., Doucet, A., Vollmer, S.J.: The bouncy particle sampler: a non-reversible rejection-free Markov Chain Monte Carlo method submitted (2015)
Duncan, A.B., Lelievre, T., Pavliotis, G.A.: Variance reduction using nonreversible Langevin samplers. J. Stat. Phys. 163(3), 457–491 (2016)
Bierkens, J., Bouchard-Cote, A., Doucet, A., Duncan, A.B., Fearnhead, P., Lienart, T., Roberts, G, Vollmer, S.J.: Piecewise deterministic Markov processes for scalable Monte Carlo on restricted domains. arxiv preprint (2018)
Lelievre, T., Rousset, M., Stoltze, G.: Free Energy Computations: a Mathematical Perspective. Imperial College Press, London (2010)
Ma, Y.-A., Fox, E.B., Chen, T., Wu, L.: Irreversible samplers from jump and continuous Markov processes. Stat Comput, 1–26 (2018)
Ottobre, M.: Markov chain Monte Carlo and irreversibility. Reports on Math Phys (2016)
Ottobre, M., Pillai, N., Pinski, F., Stuart, A.M.: A function space HMC algorithm with second order Langevin diffusion limit Bernoulli (2016)
Ottobre, M., Pillai, N., Spiliopoulos, K: Optimal scaling of the MALA algorithm with irreversible proposals for Gaussian targets. arXiv:1702.01777
Beskos, A., Pinski, F., Sanz-Serna, J.M., Stuart, A.M.: Hybrid Monte Carlo on Hilbert spaces. Stoch. Proc Appl. (2011)
Fletcher, R.: Practical Methods of Optimization, 2nd edn. Wiley-Interscience Publication (2000)
Horowitz, A.M.: A generalized guided Monte Carlo algorithm. Phys. Lett. B 268(2), 247–252 (1991)
Lolon, E.P., Archer, R.A., Ilk, D., Blasingame, T.A.: New semi-analytical solutions for multilayer reservoirs. SPE 114946
Poe, B.D. Jr., Atwood, W.K., Kohring, J., Brook, K.: Characterization of multilayer reservoir properties using production logs. SPE 101721
Popa, C., Popa, A., Cover, A.: Zonal allocation and increased production opportunities using data mining in Kern River . SPE 90266
Neal, R.M.: MCMC using Hamiltonian dynamics handbook of Markov chain Monte Carlo (2010)
Rey-Bellet, L., Spiliopoulos, K.: Irreversible Langevin samplers and variance reduction: A large deviations approach. Nonlinearity 28(7), 2081–2103 (2015)
Rey-Bellet, L., Spiliopoulos, K.: Variance reduction for irreversible Langevin samplers and diffusion on graphs. Electron. Commun. Probab., 20 (2015)
Rodrigues, J.R.P.: Calculating derivatives for automatic history matching. Computational Geosciences (2006)
Sanz-Serna, J.M.: Markov Chain Monte Carlo and Numerical Differential Equations. Current Challenges in Stability Issues for Numerical Differential Equations, pp 39–88. Springer, Cham (2014)
Villani, C.: Hypocoercivity. Mem. Amer. Math. Soc., 202 (950) (2009)
Yudin, E., Lubnin, A.: Simulation of Multilayer Wells Operating. SPE 149924
Zangl, G., Hermann, R.: Waterflood pattern optimization using genetic algorithms with multi-tank material balance. SPE 90259
Funding
The work of I. Fursov and G. J. Lord was supported by the EPSRC EQUIP grant (EP/K034154/1). P. Dobson was supported by the Maxwell Institute Graduate School in Analysis and its Applications (MIGSAA), a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix 1
This Appendix gathers some basic results about the SOL-HMC-bounce algorithm, presented in Section 2. Throughout we use the notation introduced in Section 2.
Proposition A.1
The SOL-HMC-bounce algorithm with reflections preserves the target measure.
Proof
It is easy to see that the operator \(\mathcal {O}^{\varepsilon }\) preserves the target measure \(\tilde {\pi }\). Indeed \(\mathcal {O}^{\varepsilon }\) leaves the x-variable untouched so, because \(\tilde {\pi }\) is the product of π(x) and a standard Gaussian in the p variable, looking at the definition (2.12)–(2.13) of \(\mathcal {O}^{\varepsilon }\), all one needs to show is that if p is drawn from a standard Gaussian then \(\hat {p}:=pe^{-\varepsilon } + i\xi \) is also a Gaussian random variable—here ξ is a standard Gaussian independent of p. This is readily seen as, by definition, \(\hat {p}\) has expectation 0 and variance 1, since e− 2ε + i2 = 1. Therefore, if (x,p) are drawn from \(\tilde {\pi }\), then \(\mathcal {O}^{\varepsilon }(x,p)=(x,\hat {p})\) is also distributed according to \(\tilde {\pi }\).
Let \(\chi =\chi ^{\delta }_{S,bounce}\) denote the integrator described in the SOL-HMC-bounce algorithm. Since \(\mathcal {O}^{\varepsilon }\) preserves the target measure \(\tilde {\pi }\), it remains to show that the combination of the integrator χ and the accept–reject mechanism preserves the target measure.
It is well known, for instance see [20, Theorem 9], that if the integrator \(\chi ^{\delta }_{S,bounce}\) is reversible under momentum flip (that is, \(\chi ^{\delta }_{S,bounce}\circ S= S\circ (\chi ^{\delta }_{S,bounce})^{-1}\) where S(x,p) = (x,−p)) and volume preserving, then the composition of \(\chi ^{\delta }_{S,bounce}\) and of the accept–reject move satisfies the detailed balance equation. In particular, this step also preserves the target measure \(\tilde {\pi }\).
Therefore, it is sufficient to show that \(\chi ^{\delta }_{S,bounce} = {\Theta }^{\delta /2}\circ R_{\text {bounce}} \circ {\Theta }^{\delta /2}\) is reversible under momentum flip and volume preserving. Note that both Θδ and Rδ are flows corresponding to a Hamiltonian system so they must be reversible and volume preserving, see [20, Section 8.2.2 and 8.2.3]. The composition of these operators also has these two properties and including reflection preserves these two properties, therefore Rbounce is volume preserving and reversible, and hence so is \(\chi ^{\delta }_{S,bounce}\). □
Proposition A.2
The SOL-HMC-bounce algorithm defined in Section 2 is non-reversible.
Proof Proof of Proposition A.2
For simplicity, we will only consider the case when N = 1,C = 1 and V (x) = 0. That is, we consider the target measure to be the “truncation of a standard two dimensional Gaussian”, namely
where Za is a normalising constant. In this case, Θδ/2 is the identity map, and Rδ can be written as
With these observations, if at time k the chain is (xk,pk), then we can write the proposed move \((\tilde {x}^{k+1},\tilde {p}^{k+1})\) in the k + 1-th step of SOL-HMC-bounce as
where ξ is drawn from a standard normal distribution. In this case, the acceptance probability is given by
Now we wish to calculate the transition kernel, K((x,p), (y,q)), for this Markov chain, i.e. find the probability density corresponding to the move from (x,p) to (y,q).
Observe that Rδ is a rotation about the origin and hence preserves radial distance, that is if \((\tilde {x},\tilde {p}):=R^{\delta }(x,p)\), then
Flipping momentum sign, i.e. applying reflections S, also preserves radial distance; therefore, the operator \(R_{\text {bounce}}^{\delta }\) preserves radial distance. In particular, if \(\tilde {x}^{2}+\tilde {p}^{2}< a^{2}\) (or equivalently x2 + p2 < a2), then Rδ(x,p) must remain in the strip \([-a,a]\times \mathbb {R}\), so in this situation \(R_{\text {bounce}}^{\delta }(x,p)=R^{\delta }(x,p)\).
Suppose that y2 + q2 ≤ a2. Fix some \(x\in [-a,a],p\in \mathbb {R}\), then let \((\hat {x},\hat {p})=\mathcal {O}^{\varepsilon }(x,p)=(x,pe^{-\varepsilon }+i\xi )\), where ξ is a standard normal random variable. In which case we have that \(\hat {p}\) is normally distributed with mean pe−ε and variance i2. Set \((y,q)=R^{\delta }(\hat {x},\hat {p})\), then
Therefore, y is normally distributed with mean \(x\cos \limits (\delta )+pe^{-\varepsilon }\sin \limits (\delta )\) and variance \(i^{2}\sin \limits (\delta )^{2}\). Once y has been determined, we may solve for q and find
In which case the transition kernel is given by
where α is the acceptance probability and is given by
Now the algorithm is reversible if and only if the detailed balance condition holds, that is
To see that this does not hold consider the point (x,p) = (0, 0) and let (y,q) be some point in the ball of radius a. Then by (A.1), we must have \(y=q\tan (\delta )\), and the left hand side of (A.2) becomes
On the other hand, if we suppose 0 < δ < π/4, then to move from \((q\tan (\delta ),q)\) to (0,0) is not possible unless q = 0, since (A.1) in this case becomes \( q\tan (\delta )=0 \). Therefore, for any q≠ 0, the right hand side of (A.2) must be zero, in particular we have that the algorithm is not reversible. □
Appendix 2. Description of reservoir model and simulator
The simulator we use is an in-house single-phase simulator working on an unstructured grid with finite volumes spatial discretisation and backward Euler time discretisation, calculating the dynamics of pressures and fluid flows in the subsurface porous media.
To obtain the observed pressure data, a fine grid three-phase model was run in the first place, using Schlumberger Eclipse black oil simulator [1]. The resulting output Eclipse pressures were perturbed by the uncorrelated Gaussian noise, with standard deviation σBHP = 20 bar for the well BHP data, and σb = 3 bar for the reservoir (block) pressure data. Altogether 380 measurement points were considered (365 for the BHP, 15 for the reservoir pressure), taken with time step of 6 months. The data errors’ covariance matrix Cd is diagonal, with the entries equal to either σBHP2 or σb2.
In the forward simulation mode, the reservoir properties are fixed, the producing and injecting wells (indexed by w) are controlled by the volumetric flow rates qw and the output modelled data are the time-dependent pressures at the blocks Pℓ and the bottom hole pressures at the wells \(P_{w}^{\text {BHP}}\). The equations describing the fluid flow are as follows. First, the volumetric flow rate Qℓj between the pair of connected blocks ℓ,j is proportional to the pressure difference between them, which can be regarded as Darcy’s law:
where ρ is the known liquid density, hℓj is the known depth difference between the block centers and g is the acceleration due to gravity.
The inflow qwℓ into the perforation of well w in block ℓ is proportional to the difference of the bottom hole pressure (BHP) and the block pressure:
where \(h_{w}^{\ell }\) is the depth difference between the block center and the BHP reference depth. The total inflow into well w is obtained by summing up contributions related to this well; that is,
Finally, the volumetric inflows and outflows for block ℓ are balanced, with the excessive/deficient fluid volume leading to the block pressure change via the following compressibility equation:
where t denotes time, and the first (second, respectively) summation on the right hand side is taken over all the blocks j connected to the block ℓ (all the wells w perforated in block ℓ, respectively). The compressibility cℓ of the block is supposed to be known. The simulated reservoir time spans 12 years.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Dobson, P., Fursov, I., Lord, G. et al. Reversible and non-reversible Markov chain Monte Carlo algorithms for reservoir simulation problems. Comput Geosci 24, 1301–1313 (2020). https://doi.org/10.1007/s10596-020-09947-4
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10596-020-09947-4