Gaussian Process Emulators for Computer Experiments with Inequality Constraints

Physical phenomena are observed in many fields (science and engineering) and are often studied by time-consuming computer codes. These codes are analyzed with statistical models, often called emulators. In many situations, the physical system (computer model output) may be known to satisfy inequality constraints with respect to some or all input variables. The aim is to build a model capable of incorporating both data interpolation and inequality constraints into a Gaussian process emulator. By using a functional decomposition, a finite-dimensional approximation of Gaussian processes such that all conditional simulations satisfy the inequality constraints in the entire domain is proposed. To show the performance of the proposed model, some conditional simulations with inequality constraints such as boundedness, monotonicity or convexity conditions in one and two dimensions are given. A simulation study to investigate the efficiency of the method in terms of prediction is included.     

Closed-loop reservoir management on the Brugge test case

This paper proposes an augmented Lagrangian method for production optimization in which the cost function to be maximized is defined as an augmented Lagrangian function consisting of the net present value (NPV) and all the equality and inequality constraints except the bound constraints. The bound constraints are dealt with using a trust-region gradient projection method. The paper also presents a way to eliminate the need to convert the inequality constraints to equality constraints with slack variables in the augmented Lagrangian function, which greatly reduces the size of the optimization problem when the number of inequality constraints is large. The proposed method is tested in the context of closed-loop reservoir management benchmark problem based on the Brugge reservoir setup by TNO. In the test, we used the ensemble Kalman filter (EnKF) with covariance localization for data assimilation. Production optimization is done on the updated ensemble mean model from EnKF. The production optimization resulted in a substantial increase in the NPV for the expected reservoir life compared to the base case with reactive control.     

Optimal lower bound of passive earth pressure using finite elements and non-linear programming

The present study pertains to the finding of the lower bound solution, formulating it as a non-linear programming problem using the generalized method developed by Lysmer with certain variations to incorporate the non-linear no-yield condition constraints directly in the analysis. The method considers the family of plane stress fields having the property that all stresses vary linearly within each triangular element of some mesh which covers the soil mass under study. For this type of stress field it is possible to express all equilibrium conditions as a set of linear constraints and the no-yield as a set of non-linear constraints. The boundary condition constraints may be of linear equality or inequality type. By expressing some of the design variables in terms of the remaining variables the linear equality constraints are implicity satisfied. Such a technique minimizes the complexity of the problem by eliminating the equality constraints and reduces the dimensionality of the problem, saving much, computational effort. The optimal lower bound is isolated by formulating it as a non-linear programming (NLP) problem subjected to both linear and non-linear inequality constraints. The sequential unconstrained minimization technique using the extended penalty function method as suggested by Kavlie has been used to isolate the optimal lower bound. The method has successfully been applied to the passive earth pressure and bearing capacity problem. Numerical results are obtained and compared with Lysmer's solution to show the effectiveness of the present approach.     

Geochemical and flow modelling as tools in monitoring managed aquifer recharge

Due to a growing world population and the effects of anthropogenic climate change, access to clean water is a growing global concern. Managed aquifer recharge (MAR) is a method that can help society's response to this increasing demand for pure water. In MAR, the groundwater resources are replenished and the quality of the recharged surface water is improved through effects such as the removal of organic matter. This removal occurs through mechanisms such as microbial decomposition, which can be monitored by studying the isotopic composition of dissolved inorganic carbon (DIC). Nevertheless, the monitoring can be difficult when there are other factors, like dissolving calcite, affecting the isotopic composition of DIC.The aims of this study were to establish a method for monitoring the decomposition of organic matter (dissolved organic carbon – DOC) in cases where calcite dissolution adds another component to the DIC pool, and to use this method to monitor the beginning and amount of DOC decomposition on a MAR site at Virttaankangas, southwestern Finland. To achieve this, we calculated the mean residence times of infiltrated water in the aquifer and the fractions of this water reaching observation wells. We conducted geochemical modelling, using PHREEQC, to estimate the amount of DOC decomposition and the mineral reactions affecting the quality of the water.     

Cache-efficient parallel eikonal solver for multicore CPUs

Numerical solution of the eikonal equation is frequently used to compute first-arrival travel times for a given velocity model in seismic applications. Computations for large three-dimensional models become expensive requiring the use of efficient parallel solvers. We present new parallel implementations of the fast sweeping and locking sweeping methods optimized for shared memory systems such as multicore CPUs; we call them block fast sweeping method (BFSM) and block locking sweeping method (BLSM). Proposed methods are based on the domain decomposition approach with a special attention paid to high efficiency of the cache utilization and task execution synchronization. Performance tests on realistic models show high parallel efficiency of 85–95% on modern multicore CPUs and require the same number of iterations to converge as do the serial sweeping methods. We also highlight the importance of properly selecting the stopping criterion in the iterative sweeping methods aiming for a balance between computational time and accuracy of the result required by an application. In particular, we show that in seismic applications one can reach reasonable accuracy of computed travel times while dramatically reducing the number of iterations compared to the case of using the full convergence stopping criterion.     

Strain fields and mechanical response of a highly to medium fissured bentonite clay

This article presents the developments of an ongoing research aimed at modelling the influence of fissuring on the behaviour of clays. In particular, it recalls the main results of an extensive laboratory investigation on a fissured bentonite clay from the south of Italy and presents the data of a new investigation on the evolution with shearing of the strain fields developing within the clay, resulting from Digital Image Correlation (DIC). Element test results are analysed in the framework of continuum mechanics and linked to the clay fissuring features, once characterised using the Fissuring IDentity (F‐ID) chart. This article compares the bentonite behaviour with that of other fissured clays of different F‐IDs, highlighting the common behavioural features. Thereafter, the soil response at the macro level is related to the DIC‐derived strain fields evolving within the clay with loading. For this purpose, DIC was successfully used to investigate the deformation processes active in the fissured clay and the sources of the localisation phenomena. DIC is shown to provide indications of the extent to which highly to medium fissured clays element test results can be of use to model the clay behaviour according to continuum mechanics. Copyright © 2012 John Wiley & Sons, Ltd.     

Interlaboratory study of a method for determining nonvolatile organic carbon in aquifer materials

The organic carbon fraction in aquifer materials exerts a major influence on the subsurface mobilities of organic and organic-associated contaminants. The spatial distribution of total organic carbon (TOC) in aquifer materials must be determined before the transport of hydrophobic organic pollutants in aquifers can be modeled accurately. Previous interlaboratory studies showed that it is difficult to measure TOC concentrations <0.1% in aquifer materials, when total inorganic carbon (TIC) concentrations are >1%. We have tested a new analytical method designed to improve the accuracy and precision of nonvolatile TOC quantitation in geologic materials that also contain carbonate minerals. Four authentic aquifer materials and one NIST standard reference material were selected as test materials for a blind collaborative study. Nonvolatile TOC in these materials ranged from 0.05 to 1.4%, while TIC ranged from 0.46 to 12.6%. Sample replicates were digested with sulfurous acid, dried at 40°C, and then combusted at 950°C using LECO or UIC instruments. For the three test materials that contained >2% TIC, incomplete acidification resulted in a systematic positive bias of TOC values reported by five of the six laboratories that used the test method. Participants did not have enough time to become proficient with the new method before they analyzed the test materials. A seventh laboratory successfully used an alternative method that analyzed separate liquid and solid fractions of the acidified sample residues.     

Simulating Large Gaussian Random Vectors Subject to Inequality Constraints by Gibbs Sampling

The Gibbs sampler is an iterative algorithm used to simulate Gaussian random vectors subject to inequality constraints. This algorithm relies on the fact that the distribution of a vector component conditioned by the other components is Gaussian, the mean and variance of which are obtained by solving a kriging system. If the number of components is large, kriging is usually applied with a moving search neighborhood, but this practice can make the simulated vector not reproduce the target correlation matrix. To avoid these problems, variations of the Gibbs sampler are presented. The conditioning to inequality constraints on the vector components can be achieved by simulated annealing or by restricting the transition matrix of the iterative algorithm. Numerical experiments indicate that both approaches provide realizations that reproduce the correlation matrix of the Gaussian random vector, but some conditioning constraints may not be satisfied when using simulated annealing. On the contrary, the restriction of the transition matrix manages to satisfy all the constraints, although at the cost of a large number of iterations.     

A multi-domain decomposition-based Fourier finite element method for the simulation of 3D marine CSEM measurements

We introduce a multi-domain decomposition Fourier finite element (MDDFFE) method for the simulation of three-dimensional (3D) marine controlled source electromagnetic measurement (CSEM). The method combines a 2D finite element (FE) method in two spatial dimensions with a hybrid discretization based on a Fourier FE method along the third dimension. The method employs a secondary field formulation rather than the total field formulation. We apply the MDDFFE method to several synthetic marine CSEM examples exhibiting bathymetry and/or multiple 3D subdomains. Numerical results show that the use of the MDDFFE method reduces the problem size by as much as 87 % in terms of the number of unknowns, without any sacrifice in accuracy.     

First outcomes in the definition of groundwater protection zones at the Viñales National Park (Cuba) and surrounding area

Due to the existence of several pollution hazards in the area and its surroundings, groundwater protection in the Viñales National Park has been defined as a priority task in its management plan for 5 years (2009–2013). Those hazards affect directly the surface waters and the karst nature of the territory allowing the contaminants to enter directly into the aquifers without any or at least small self-depuration processes. In turn, those processes can affect very large areas in a short time. In this paper several results of precedent studies related to the vulnerability, hazards, risk, and the mapping of the spatial variability of natural recharge are integrated in a GIS analysis. The methodology is based on the method of Jimenez-Madrid et al. which is slightly modified in this contribution. The analysis of hydrographs, chemographs, mathematical models and tracer tests is necessary for the definition of a response plan against any contamination events and for model validation, together with the incorporation of data related to the energetic responses from the aquifer systems to external forcing.     

基于数字图像相关法的两类岩石断裂特征研究

岩石破坏时的临界断裂特征,如过程区长度,裂缝口张开位移是运用断裂力学解决岩石断裂问题的关键所在。传统测量方法,如应变片、直线位移传感器无法获得岩石破坏时的全场变形,因此,也无法准确地获得上述断裂特征。数字图像相关法是一种光学的变形测量方法,其通过试样表面的数字图像采集及相关计算,能够获得岩石破裂过程各个阶段的高精度全场变形特征。利用该试验手段,对两类岩石,即相对较硬的大理岩和相对较软的黄砂岩开展了一系列半圆盘三点弯曲断裂试验并获得了岩石断裂时的临界变形场,对变形场进行分析,从而确定了两类岩石破坏时临界特征、过程区长度及裂缝口张开位移。结果表明,大理岩的断裂过程区长度明显小于黄砂岩的断裂过程区长度,峰值时黄砂岩COD的值均大于相同裂纹长度的大理岩,而较软岩的力学行为更容易受边界效应的影响。上述研究有助于进一步了解不同类型岩石的断裂发展过程并采用相应的方法来解决岩石破裂问题。     

南海重力异常特征及其显著的构造意义

在南海地区地震测深数据有限的情况下,利用重力异常可以研究南海大范围的深部地壳结构及地质构造展布特征。基于空间重力异常,结合最新的地形、沉积物厚度及地震测深等数据,分别从地震约束的莫霍面反演和无约束的三维相关成像两个视角研究南海的地壳结构,利用壳幔界面起伏、地壳厚度及三维等效密度分布来探讨地壳结构的纵横向变化。同时,联合采用延拓、水平梯度及线性构造增强滤波方法聚焦重力异常中的区域线性特征,突出显示了反映地壳横向变化的深断裂、洋陆转换边界、海盆扩张轴等线性构造的展布。重力解释与贯穿南海南北的广州-巴拉望地学断面对比表明,重力异常反演及异常的区域线性特征,较好地揭示了南海海域大范围的地壳结构与区域构造展布。     

Automatic extraction of sparse trees from high-resolution ortho-images

Obtaining information about tree species distribution in agricultural lands is a topic of interest for various applications, such as tree inventory, forest management, agricultural land management, crop estimation, etc. This information can be derived from images obtained from modern remote sensing technology, which is the most economical way as compare to field surveys covering large geographic areas. Therefore, in this study, a new method is proposed for extraction and counting of sparse and regular distributed individual pistachio trees from agricultural areas on large scale from high-resolution digital ortho-photo maps, which were obtained using an airborne sensor (Ultracam-X). The input images were first smoothed by applying Gaussian filter to reduce the impact of noise. Normalized difference vegetation indices (NDVI) were then derived to obtain vegetation areas followed by Otsu’s global thresholding algorithm to obtain candidate tree areas. Further, connected component (CC) analysis was applied to segregate each object. Morphological processing was performed to fill holes within tree objects and get smooth contours, which were obtained by using the Moore-neighbor tracing method (MNTM) for each CC, while geometrical constraints were applied to undermine possible non-tree elements from output image. To further improve the segmentation results for sparse trees, a new method was applied, called quadratic local analysis (QLA). QLA helped to segment the trees, which were missed by the Otsu method due to low contrast and resulted in improved accuracy (3–6%). The obtained results were compared with well-known support vector machine (SVM) classifier. Proposed method produced slightly better results (1–5%) than SVM for extraction of pistachio trees and obtained accuracy for QLA and SVM were 96 and 91% for region 1, while 91 and 90% for region 2 respectively.     

Optimised design of jet-grouted raft using response surface method

The use of jet grouting as a foundation supporting element has brought about the need to define the design strategies for jet-grouted raft (JGR) system which differs from the traditional raft design due to the effects of jet-grouted columns (JGCs). Paper tackles this important aspect by combining the previously defined design strategies on piled raft and JGC to achieve an optimised design of complete JGR system. Initially, a single JGC has been analysed by means of three dimensional finite element (3D FE) analysis accounting for the previously measured actual trial JGC’s geometrical variation with depth. The image processing technique allowing the complex 3D FE modelling of JGC system is utilised and the results are validated by the back-analysis of the well-known experimental results reported in the literature. In a parametric study, extension of the numerical analysis to the complete JGR system is accomplished by considering the mutual interaction between the foundation elements. The effects of design factors that are interlayer thickness, JGC spacing and lengths on the design responses of vertical stresses, bending moments, average and differential settlements are presented. Response Surface Method is utilised in the multi objective optimisation analysis to present the optimised design solution by accounting for the design constraints previously defined for the considered factors and the responses. The significance of design constraints and their relative influences on the optimised design of JGR are finally discussed.     

A simultaneous perturbation stochastic approximation algorithm for coupled well placement and control optimization under geologic uncertainty

Development of subsurface energy and environmental resources can be improved by tuning important decision variables such as well locations and operating rates to optimize a desired performance metric. Optimal well locations in a discretized reservoir model are typically identified by solving an integer programming problem while identification of optimal well settings (controls) is formulated as a continuous optimization problem. In general, however, the decision variables in field development optimization can include many design parameters such as the number, type, location, short-term and long-term operational settings (controls), and drilling schedule of the wells. In addition to the large number of decision variables, field optimization problems are further complicated by the existing technical and physical constraints as well as the uncertainty in describing heterogeneous properties of geologic formations. In this paper, we consider simultaneous optimization of well locations and dynamic rate allocations under geologic uncertainty using a variant of the simultaneous perturbation and stochastic approximation (SPSA). In addition, by taking advantage of the robustness of SPSA against errors in calculating the cost function, we develop an efficient field development optimization under geologic uncertainty, where an ensemble of models are used to describe important flow and transport reservoir properties (e.g., permeability and porosity). We use several numerical experiments, including a channel layer of the SPE10 model and the three-dimensional PUNQ-S3 reservoir, to illustrate the performance improvement that can be achieved by solving a combined well placement and control optimization using the SPSA algorithm under known and uncertain reservoir model assumptions.     

Forward Calculation of 2-D and 3-D Structures with a Cover by the Boundary Element Method while Using Electrical Methods

This paper develops the boundary element method,the authors employtwo-layered earth Green’s functions as the weighting functions of residual andderive boundary integral equations.The forward problems of point sources on2-D and 3-D structures with an influencing cover are solved by this method.The results show that this method markedly improves the original boundaryelement method.The features of the improved method are greater numericalaccuracy and much smaller systems of equations and thus considerable savingsfor the storage capacity of computers,allowing us to solve the above problemswith only ordinary microcomputers.The results in this paper extend the scopeof applying the boundary element method while using electrical methods forgeophysical prospecting.     

Kriging with Inequality Constraints

A Gaussian random field with an unknown linear trend for the mean is considered. Methods for obtaining the distribution of the trend coefficients given exact data and inequality constraints are established. Moreover, the conditional distribution for the random field at any location is calculated so that predictions using e.g. the expectation, the mode, or the median can be evaluated and prediction error estimates using quantiles or variance can be obtained. Conditional simulation techniques are also provided.     

模拟三维裂纹问题的多尺度扩展有限元法

扩展有限元法模拟裂纹时独立于网格,因此该方法是目前求解裂纹问题最有效的数值方法。为了在计算代价不大的情况,实现大型结构分析中考虑小裂纹或提高裂纹附近精度,在裂纹附近一般采用小尺度单元,其他区域采用大尺度单元。提出了分析三维裂纹问题的多尺度扩展有限元法,在需要的地方采用小尺度单元。基于点插值构造了六面体任意节点单元。所有尺度单元都采用8节点六面体单元,这样六面体任意节点单元可方便有效地连接不同尺度单元。采用互作用积分法计算三维应力强度因子。边裂纹和中心圆裂纹算例分析结果表明,该方法是正确和有效的。     

An interpolation method taking into account inequality constraints: II. Practical approach

Common features of models for interpolation, consistent with a finite number of inequality constraints on the range of values of a variablez, are discussed. A method based on constrained quadratic minimization yielding kriging estimates when no constraints exist, is presented. A computationally efficient formulation of quadratic minimization is obtained by using results on duality in quadratic programming. Relevant properties of the optimal interpolator are derived in a simple, self-contained way. The method is applied to mapping of horizon depth and estimation of thickness of an oil-bearing formation.