Performance comparison of multiple and single surrogate models for pumping optimization of coastal aquifers

ABSTRACT

Single and multiple surrogate models were compared for single-objective pumping optimization problems of a hypothetical and a real-world coastal aquifer. Different instances of radial basis functions and kriging surrogates were utilized to reduce the computational cost of direct optimization with variable density and salt transport models. An adaptive surrogate update scheme was embedded in the operations of an evolutionary algorithm to efficiently control the feasibility of optimal solutions in pumping optimization problems with multiple constraints. For a set of independent optimization runs, results showed that multiple surrogates, either by selecting the best or by using ensembles, did not necessarily outperform the single surrogate approach. Nevertheless, the ensemble with optimal weights produced slightly better results than selecting only the best surrogates or applying a simple averaging approach. For all cases, the computational cost, by using single or multiple surrogate models, was reduced by up to 90% of the direct optimization.

KEYWORDS:

• coastal aquifers
• pumping optimization
• surrogate models
• multiple surrogates
• computational cost

Editor R. Woods; Associate editor L. Bouchaou

Introduction

Variable density and salt transport (VDST) numerical models provide high-fidelity simulations of seawater intrusion (SWI) at the expense of increased computational cost (Werner et al. 2013). As a result, the integration of VDST models within simulation-based frameworks, such as uncertainty analysis or optimization, is often computationally impractical. Several SWI studies have proposed the use of surrogate models to reduce the computational cost without compromising the quality of the analysis (Ketabchi and Ataie-Ashtiani 2015a, Sreekanth and Datta 2015).

Typically, surrogate models are trained on a set of input–output data produced from the computationally expensive physics-based model and then serve as efficient empirical approximators of the original model’s response (Razavi et al. 2012a, Asher et al. 2015). The selection of the surrogate model is often indicated by the type of analysis required by the expensive computer model (Broad et al. 2015). For example, Rajabi et al. (2015) employed polynomial chaos expansions as appropriate surrogate models to cope with uncertainty in SWI simulations. Artificial neural networks, genetic programming, radial basis functions, fuzzy inference systems and multivariate adaptive regression splines are examples of surrogate models that have been used as surrogates for VDST models, in deterministic single-objective or multi-objective optimization problems of coastal aquifers (e.g. Rao et al. 2004, Rao and Manju 2007, Dhar and Datta 2009, Sreekanth and Datta 2010, Kourakos and Mantoglou 2011, Grundmann et al. 2012, Ataie-Ashtiani et al. 2013, Christelis and Mantoglou 2016, Roy and Datta 2017a, 2017b, 2017c). Some studies have also compared the prediction skills and the efficiency of different surrogates for approximating the response of VDST models (e.g. Sreekanth and Datta 2011a, Lal and Datta 2018, Yadav et al. 2018).

Instead of employing a single surrogate to approximate the response of an expensive computer model, various researchers in engineering optimization have investigated the use of multiple surrogates aiming to improve the exploration of the optimal search space (e.g. Viana et al. 2009, 2013, Acar 2010, Müller and Piché 2011, Nikolos 2013, Müller and Shoemaker 2014, Jiang et al. 2015, Shi et al. 2016, Bhosekar and Ierapetritou 2017, Hou et al. 2017). The generic approach in a multiple surrogate framework is to identify a suite of reliable surrogates for the problem at hand, usually through cross-validation errors, and either use the best or construct an ensemble (e.g. weighted-average surrogate) to replace the original expensive model (Viana et al. 2010). The superiority of employing multiple against single surrogates is debatable, since the nature of each optimization problem may favour the one approach or the other (e.g. Viana et al. 2009, Babaei and Pan 2016).

In addition, the use of multiple surrogates involves several considerations that are related to the effort and computation time required to implement a surrogate-based optimization (SBO). A primary issue is if the best surrogate models are identified through a cross-validation strategy, or if it is known, based on previous experience with the optimization problem, which surrogates should work best (Viana et al. 2009). In the first case, it must be decided how many surrogates should be explored and what computational budget is available to apply an efficient and informative cross-validation framework. Furthermore, the training time of surrogate models varies with their complexity, and thus employing more complex yet possibly more accurate models may add considerable computational time.

The multiple surrogate approach is limited in coastal aquifer management and has been mostly explored using ensembles constructed from different instances of the same surrogate model (Roy and Datta 2018b). Sreekanth and Datta (2011b) first utilized an ensemble of genetic programming surrogate models in lieu of a VDST model to solve a multi-objective pumping optimization problem. Their study showed that the ensemble performed better than using a single genetic programming surrogate model. Recently, Roy and Datta (2017a, 2017c) solved multi-objective pumping optimization problems by utilizing ensembles of fuzzy inference systems and of multivariate adaptive regression splines, respectively. Their approach reduced the uncertainty in the prediction of the surrogate models while the multivariate adaptive regression splines provided a more efficient ensemble surrogate selection. Roy and Datta (2018a) also proposed the use of ensembles of Gaussian process regression models to increase the prediction accuracy of surrogate models in coastal aquifer management under parameter uncertainty. To handle the computational burden, they also used parallel computing facilities.

The combination of different surrogate modelling techniques may potentially increase the accuracy in approximating the original computer model response (Viana et al. 2009, Müller and Shoemaker 2014). The present work compares the performance of single and multiple surrogates for single-objective pumping optimization problems of coastal aquifers. An adaptive SBO framework is adopted where the surrogate models are used to predict the computationally expensive constraint functions associated with multiple pumping wells. By utilizing a cross-validation error strategy, the multiple surrogate approach is applied either by identifying the best model or by forming weighted-average surrogates (ensembles) for each constraint function. Different instances of radial basis functions and kriging surrogate models are employed due to their beneficial interpolating capabilities in approximating deterministic computer simulations (Razavi et al. 2012b). The SBO frameworks are also compared with the direct optimization using the VDST model, which provided the benchmark solutions for a hypothetical and a real-world coastal aquifer management problem. To enable a more comprehensive comparison due to the probabilistic nature of the optimization algorithm, multiple optimization runs were performed for each framework.

SWI modelling and pumping optimization

Overview of the study area

An unconfined, elongated coastal aquifer along the Vathi valley, located in the central part of the Greek island of Kalymnos, is considered as the real-world application of this study (Figure 1).

Figure 1. Hydro-lithological map of Kalymnos island. Source: Hellenic Ministry of Development (2005)

Regarding the hydro-lithology of the examined aquifer, it consists mainly of highly permeable limestone, which outcrops at the valley margins. Borehole lithological data indicate that the aquifer is bounded by an impermeable schist formation, which underlies the limestones. The valley plain consists of scree, alluvium deposits and volcanic formations (tuff) of high, medium and very low permeability, respectively. Although the carbonate rocks are characterized by secondary porosity, for simplification purposes a uniform hydraulic conductivity is considered for the limestones. The aquifer is divided into four zones of uniform hydraulic conductivity, as depicted in Figure 2, and it is replenished by a uniform surface recharge of approximately 7849 m³/d.

Figure 2. Hydraulic conductivity zones of the Vathi aquifer

The aquifer was simulated assuming a horizontal bottom at a depth of −25 m. The northern, southern and western boundaries are considered impermeable. At the eastern sea boundary, first-type (Dirichlet) boundary conditions are applied. That is, a specified hydrostatic equivalent freshwater head and a constant relative concentration corresponding to the maximum fluid density of seawater are applied to represent sea boundaries. Due to the lack of hydrogeological data and hydraulic head measurements, it was not possible to perform a reliable and accurate calibration of the aquifer model. Therefore, the parameter values are rough estimates, based on a trial-and-error procedure (Mantoglou et al. 2004).

Application on a hypothetical coastal aquifer model

The optimization frameworks were all first tested for the numerical SWI simulations of an illustrative coastal aquifer model (Figure. 3). The dimensions of the coastal aquifer model are set to x = 7000 m, y = 3000 m, and the aquifer thickness is z = 25 m, while a rectangular shape is assumed. Unconfined, steady-state and saturated flow conditions are assumed. Table 1 summarizes the basic input parameters for the VDST numerical model.

Table 1. Parameters for the numerical SWI simulations. K_x, K_y, K_z: hydraulic conductivities (m/d); R_gw: total aquifer recharge (m³/d); α_L, α_T: longitudinal and transverse dispersivity (m); Δ_x, Δ_y, Δ_z: grid discretization settings (m)

Model parameters	Values
Kx, Ky, Kz	100, 100, 10
Rgw	5422
αL, αT	100, 10
Δx, Δy, Δz	100, 100, 5

Figure 3. Illustrative coastal aquifer model with boundary conditions applied

A CPU time of 30 s is required for a single run of the VDST model. All the simulation–optimization methods were conducted on the same desktop, namely, a 2.7-GHz Intel i5 processor with 8 GB of RAM in a 64-bit Windows 10 system. An initial simulation run of the VDST model was performed with no pumping present, until the head and salinity concentration fields reached steady state. These outputs provided the initial conditions for the subsequent VDST simulations related to the optimization runs.

VDST models

The phenomenon of SWI is more realistically described by VDST models which consider mixing between freshwater and salt water due to hydrodynamic dispersion. The VDST modelling is based on numerical codes which solve a coupled system of partial differential equations of flow and transport:

(1) $\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left[K_{ij}\left(\frac{\mathrm{\partial }{h}_f}{\mathrm{\partial }{x}_j}+{\rho }_r{n}_j\right)\right]+Q_{\rho }=S_s\frac{\mathrm{\partial }{h}_f}{\mathrm{\partial }t}$(1)

(2) $\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(D_{ij}\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}-q_{i}c\right)+Q_c=\frac{\mathrm{\partial }\left(\varphi c\right)}{\mathrm{\partial }t}$(2)

In the flow Equation (1)(1) $\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left[K_{ij}\left(\frac{\mathrm{\partial }{h}_f}{\mathrm{\partial }{x}_j}+{\rho }_r{n}_j\right)\right]+Q_{\rho }=S_s\frac{\mathrm{\partial }{h}_f}{\mathrm{\partial }t}$(1) , h_f [L] represents the equivalent freshwater head given by $h_f=\left(p/{\rho }_{f}g\right)+z$; K_ij [L T⁻¹] are the coefficients of the freshwater hydraulic conductivity tensor; S_s [L⁻¹] is the specific storage; t [T] is time; and Q_ρ [L³ L⁻³ T⁻¹] is a volumetric fluid source/sink term per unit aquifer volume. The gravity acceleration constant is denoted by g [M L⁻²] and z [L] is the elevation above the horizontal datum. The fluid pressure is denoted by p [M L⁻¹ T⁻²], while ρ_f [M L⁻³] is the reference fluid density. A linear relation is assumed between fluid density ρ [M L⁻³] and concentration c according to:

(3) $\left(\rho -{\rho }_f\right)/{\rho }_f=\left[\left({\rho }_{max}-{\rho }_f\right)/{\rho }_f\right]c$(3)

where c [-] is the dimensionless relative concentration, which varies between 0 and 1. Thus, ρ = ρ_max corresponds to c_max = 1. In the transport Equation (2)(2) $\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(D_{ij}\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}-q_{i}c\right)+Q_c=\frac{\mathrm{\partial }\left(\varphi c\right)}{\mathrm{\partial }t}$(2) , D_ij [L² T⁻¹] are the coefficients of the dispersion tensor; ϕ [-] is porosity and Q_c [M L⁻³ T⁻¹] is a solute mass source/sink term. Note that the indices i,j are the unit vectors in the x and y directions and n_j is the direction of flow, which is equal to 1 in the vertical direction and 0 in the horizontal directions. The Darcy flux term q_i is expressed as:

(4) $q_i=-K_{ij}\left(\frac{\mathrm{\partial }{h}_f}{\mathrm{\partial }{x}_j}+{\rho }_r{n}_j\right)$(4)

where ρ_r is the dimensionless relative density ((ρ – ρ_f)/ρ_f). In this study, HydroGeoSphere (HGS) code (Graf and Therrien 2005, Therrien et al. 2006) was used to simulate SWI for the case of a hypothetical and a real-world coastal aquifer model. The numerical solution is based on the control volume finite element method and adaptive time-stepping (Therrien et al. 2006). A Picard scheme is utilized by HGS to iterate between the flow and transport equations.

Pumping optimization based on SWI models

Let $f\left(\mathit{Q}\right)=\sum _{i=1}^M{Q}_i$ represent the objective function of the problem, which is linear in respect to the decision variables Q_i, with i = 1, 2, …, M; that is, the individual pumping rates of each pumping well, where Q is the decision vector of pumping rates (Q = (Q₁, Q₂, …, Q_M)). The objective is to find the maximum of f (total groundwater abstraction) subject to inequality constraints of the form g_i(Q) ≤ 0, which protect the coastal aquifer from SWI. Here, it is assumed that the salinity concentration c_i in each pumping well is not allowed to exceed a specified salinity concentration threshold c_T = 35 mg/L. This is formulated as a nonlinear constrained optimization problem as follows:

(5) $\begin{array}{c}min-\sum _{i=1}^M{Q}_i\\ \mathrm{s}.\mathrm{t}.c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T,\mathrm{\forall }i=1,2,\dots ,M\\ Q_{min}\le Q_i\le Q_{max},i=1,2,\dots ,M\end{array}$(5)

Note that the number of constraint functions is equal to the number of pumping wells. In Equation (5)(5) $\begin{array}{c}min-\sum _{i=1}^M{Q}_i\\ \mathrm{s}.\mathrm{t}.c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T,\mathrm{\forall }i=1,2,\dots ,M\\ Q_{min}\le Q_i\le Q_{max},i=1,2,\dots ,M\end{array}$(5) , Q_min and Q_max define the lower and upper limits that pumping rates can take. The negative sign in the formulation of the objective function denotes that the objective function is maximized, since the location of the maximum of f(Q) occurs at the same point, with the minimum of −f(Q).

It has been demonstrated in several coastal aquifer management studies that evolutionary algorithms can successfully handle the multimodality of pumping optimization problems in coastal aquifers (e.g. Karpouzos and Katsifarakis 2013, Ketabchi and Ataie-Ashtiani 2015b). Here, a heuristic optimization method, namely, the evolutionary annealing–simplex (EAS) algorithm (Efstratiadis and Koutsoyiannis 2002, Rozos et al. 2004), is utilized to solve problem (5). The latter was translated to a bound-constrained optimization problem using penalty terms in the objective function so that EAS can be applied, as given by Equation (6)$minf\left(\mathit{Q}\right)=\left\{\begin{array}{c}-\sum _{i=1}^M{Q}_i\mathrm{i}\mathrm{f}\mathrm{\forall }i=1,2,\dots ,M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T\\ M_v\sum _{i=1}^M{\left[max\left(\left(c_i-c_T\right),0\right)\right]}^2\mathrm{i}\mathrm{f}\mathrm{\exists }i=1,2,\dots .M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)>c_T\end{array}\right.(6)$:

$minf\left(\mathit{Q}\right)=\left\{\begin{array}{c}-\sum _{i=1}^M{Q}_i\mathrm{i}\mathrm{f}\mathrm{\forall }i=1,2,\dots ,M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T\\ M_v\sum _{i=1}^M{\left[max\left(\left(c_i-c_T\right),0\right)\right]}^2\mathrm{i}\mathrm{f}\mathrm{\exists }i=1,2,\dots .M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)>c_T\end{array}\right.(6)$

where M_v represents the number of pumping wells for which the constraint is violated. Equation (6)$minf\left(\mathit{Q}\right)=\left\{\begin{array}{c}-\sum _{i=1}^M{Q}_i\mathrm{i}\mathrm{f}\mathrm{\forall }i=1,2,\dots ,M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T\\ M_v\sum _{i=1}^M{\left[max\left(\left(c_i-c_T\right),0\right)\right]}^2\mathrm{i}\mathrm{f}\mathrm{\exists }i=1,2,\dots .M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)>c_T\end{array}\right.(6)$ attributes a separate score for each violated constraint. Furthermore, the penalized objective function is multiplied by M_v to include the number of constraint violations for the case of a non-feasible vector Q. The pumping optimization problem defined in Equation (6)$minf\left(\mathit{Q}\right)=\left\{\begin{array}{c}-\sum _{i=1}^M{Q}_i\mathrm{i}\mathrm{f}\mathrm{\forall }i=1,2,\dots ,M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T\\ M_v\sum _{i=1}^M{\left[max\left(\left(c_i-c_T\right),0\right)\right]}^2\mathrm{i}\mathrm{f}\mathrm{\exists }i=1,2,\dots .M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)>c_T\end{array}\right.(6)$ can be directly solved by combining the VDST model with the EAS algorithm. Guidelines for the EAS parameters can be found in Efstratiadis and Koutsoyiannis (2002) and Tsoukalas et al. (2016). Accordingly, the initial population size was set to n_pop = 10M, the two annealing schedule control parameters to λp = 0.95 and ψ = 2, the mutation probability to mp = 0.1, and the convergence criterion to ε = 10⁻⁴. For all optimization frameworks, the termination criteria are met if the convergence criterion ε equals its pre-set value, or the number of objective function evaluations reaches the maximum of n_max = 100n_pop.

Surrogate-based pumping optimization

The surrogate models

The nonlinear constraints described in Equation (5)(5) $\begin{array}{c}min-\sum _{i=1}^M{Q}_i\\ \mathrm{s}.\mathrm{t}.c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T,\mathrm{\forall }i=1,2,\dots ,M\\ Q_{min}\le Q_i\le Q_{max},i=1,2,\dots ,M\end{array}$(5) are computationally expensive to evaluate, since the VDST runs first to calculate the salinity concentration field at the end of the simulated management period. Given the computational cost of VDST simulations and the thousands of runs required by an evolutionary algorithm for a moderate dimensionality optimization problem, the direct approach may result in excessive computational burden. To reduce this computational burden, surrogate models are employed to replace the VDST simulation and to enable an efficient simulation–optimization routine. Furthermore, it is assumed that there is no a priori information for selecting a specific type of surrogate model, but instead it is computationally feasible to explore different surrogate model formulations for implementing a SBO framework. After an appropriate training procedure, the surrogate models can approximate the original salinity output at the corresponding observation points of the VDST model.

A suite of radial basis functions and kriging models is considered here to investigate the efficacy of multiple surrogates for the coastal aquifer management problem defined in the previous section. Both radial basis functions and kriging surrogate models have been successfully applied in water resources optimization problems (e.g. Baú and Mayer 2006, Shoemaker et al. 2007, Razavi et al. 2012b, Pan et al.2014, Tsoukalas and Makropoulos 2015, Christelis and Mantoglou 2016, Tsoukalas et al. 2016, Christelis et al. 2016, 2018). As interpolating surrogate models, they pass through all the previously evaluated points with the original model and thus can be more accurate as new input–output data become available from deterministic computer simulations (Forrester et al. 2008, Wang et al. 2014). A Latin hypercube sampling (LHS) design (McKay et al. 1979) was utilized to uniformly sample the decision vector space, evaluate the VDST model and create an initial set of m training patterns for each surrogate model.

Radial basis functions (RBF)

For the construction of the RBF surrogate models we used the codes developed within the MATSuMoTo MATLAB toolbox developed by Müller (2014). A unique RBF model is associated with each of the M pumping wells. The set of decision vectors ${\mathit{Q}}^{\left(1\right)},{\mathit{Q}}^{\left(2\right)},\dots ,{\mathit{Q}}^{\left(m\right)}\in R^M$ obtained from the LHS design and the values $c_i\left({\mathit{Q}}^{\left(1\right)}\right),c_i\left({\mathit{Q}}^{\left(2\right)}\right),\dots ,c_i\left({\mathit{Q}}^{\left(m\right)}\right)$ (i = 1, …, M) obtained from the VDST model define a RBF model, with a linear polynomial tail, of the following form (Powell 1992):

(7) $S_m\left(\mathit{Q}\right)=\sum _{k=1}^m{\lambda }_k\varphi \left(\left|\left(\mathit{Q}-{\mathit{Q}}^{\left(k\right)}\right)\right|\right)+p\left(\mathit{Q}\right)$(7)

where ${\lambda }_1,\dots ,{\lambda }_m\in R$ are coefficients to be determined, and p(Q) is a linear polynomial whose coefficients also need to be determined such that S_m passes through all the design points. Here, three types of basis functions are considered, that is: the linear form where $\varphi \left(r\right)=r$, the cubic form where $\varphi \left(r\right)=r^3$ and a thin plate spline (TPS) where $\varphi \left(r\right)=r^{2}lnr$. The construction time for the RBF models is negligible in comparison to the VDST simulation time even for a large number of training patterns.

Kriging (KRG)

The MATLAB toolbox DACE (Design and Analysis of Computer Experiments) (Lophaven et al. 2002) was used to construct KRG surrogate models. KRG was first introduced as a surrogate model of deterministic computer simulations in Sacks et al. (1989). Details about the mathematical background of KRG can be found elsewhere in the literature (e.g. Jones et al. 1998, Lophaven et al. 2002, Forrester et al. 2008). Here, we only briefly present the general concept for a KRG surrogate model as applied in our problem.

Kriging treats the deterministic outputs $c_i\left({\mathit{Q}}^{\left(1\right)}\right),$ $c_i\left({\mathit{Q}}^{\left(2\right)}\right),\dots ,c_i\left({\mathit{Q}}^{\left(m\right)}\right)$ (i = 1, …, M) obtained from the VDST model as if they were generated from a stochastic process defined by a known regression function f(Q) and a zero mean Gaussian process Z(Q), with variance σ² and a correlation model R (Jones et al. 1998):

(8) $Y\left(\mathit{Q}\right)=f\left(\mathit{Q}\right)+Z\left(\mathit{Q}\right)$(8)

The concept utilized by the KRG models is that, if the difference between two vectors Q and Q′ is small, then the scalar resulting predictions $\stackrel{\wedge }{c_i}\left(\mathit{Q}\right)$ and $\stackrel{\wedge }{c_i}\left({\mathit{Q}}^{\prime }\right)$ of the “true”, black-box function (the VDST numerical model in this case) should be closely correlated (Forrester et al. 2008). The DACE toolbox allows for both regression models of different order polynomials and correlation models of different structures. The correlation model involves a set of parameters θ, which are identified using the maximum likelihood estimation method and numerical optimization techniques (Forrester et al. 2008). Here, the initial guess for θ values was set according to Viana (2011). In total six KRG surrogate models were utilized by selecting exponential and Gaussian correlation models combined with zero-order, first-order and second-order polynomial models.

Adaptive SBO using single surrogates

As previously mentioned, a total of nine surrogate models (three RBF and six KRG) were considered in this study. Each was used separately to develop an adaptive SBO framework similar to that presented in Christelis and Mantoglou (2016). The steps of the SBO approach can be summarized as follows:

1. Use LHS design to provide the initial m training points ${\mathit{Q}}^{\left(1\right)},{\mathit{Q}}^{\left(2\right)},\dots ,{\mathit{Q}}^{\left(m\right)}$ and obtain $c_i\left({\mathit{Q}}^{\left(1\right)}\right),$ $c_i\left({\mathit{Q}}^{\left(2\right)}\right),\dots ,c_i\left({\mathit{Q}}^{\left(m\right)}\right)$(i = 1, …, M) through m VDST simulations.

2. Construct M surrogate models (e.g. M cubic RBFs) corresponding to the M pumping wells, and create an external archive of training patterns.

3. Run the EAS algorithm based on the surrogate models and, if during optimization a better optimum is found, do the following:

4. Re-evaluate the current decision vector with the VDST model;

5. Replace the objective function value with that obtained from the VDST model;

6. Store the new input–output data to the archive and re-train the surrogate models.

7. Are the stopping criteria for EAS algorithm met? If yes, return the final solution; otherwise go back to Step 3.

Whereas the above adaptive strategy favours the fast improvement of the surrogate model at the region of the current optimum (local exploitation), it is not considered as a global SBO approach. However, it quickly locates good optimal solutions and has been successfully applied in other single-objective pumping optimization studies of coastal aquifers (Kourakos and Mantoglou 2009, Papadopoulou et al. 2010, Christelis and Mantoglou 2016). Another very active field of research found in the adaptive SBO literature is the development of sampling strategies that effectively utilize the expensive original model simulations to update the surrogate and increase its global accuracy (e.g. Regis 2011, 2014, Zhou et al. 2017).

Adaptive SBO using multiple surrogates

A cross-validation strategy is typically employed to assess the performance of the surrogate models and then either the best one is selected, or a set of the best available surrogates is used to construct an ensemble. Cross-validation is a practical approach to obtain information about the accuracy of the surrogates based on an initial training set. However, whether the combination provides better results over the selection of a single best surrogate model is inconclusive in the existing literature (Viana et al. 2010).

Three approaches are followed to make use of the multiple surrogate concept. First, given a relatively small number of m input–output patterns from the VDST model, a leave-one-out strategy is employed to find the best surrogate model. Here, the available surrogates emulate the scalar response of the computer model and thus a total of M surrogate models need to be constructed for each of the constraint functions corresponding to the pumping wells (Equation (5)(5) $\begin{array}{c}min-\sum _{i=1}^M{Q}_i\\ \mathrm{s}.\mathrm{t}.c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T,\mathrm{\forall }i=1,2,\dots ,M\\ Q_{min}\le Q_i\le Q_{max},i=1,2,\dots ,M\end{array}$(5) ). That is, a cross-validation error is calculated by fitting the available surrogate models to m − 1 training points; then the response of the surrogate models to the point left out from the fitting process is predicted. This process is repeated m times until all points have been evaluated separately. The root mean squared error metric is used as a cross-validation score, which is defined as follows for the ith pumping well:

(9) $\mathrm{R}\mathrm{M}\mathrm{S}{\mathrm{E}}^{\left(i\right)}=\sqrt{\frac{1}{m}\sum _{j=1}^m{\left(c_i^j\left(Q\right)-{\hat{c}}_i^j\left(Q\right)\right)}^2},i=1,\dots ,M$(9)

The surrogate models with the lowest RMSE values for each pumping well are identified and then used in the adaptive SBO scheme discussed above. Within the adaptive SBO framework, the surrogate models are reconstructed as new training points become available through the optimization operations.

Apart from selecting the best surrogate model, two ensemble formulations are also implemented. The two best surrogate models identified for each constrained function from the cross-validation strategy are used to form an ensemble of surrogates (ESM). Here, the multiple surrogates are formed as ensembles by constructing a weighted-average surrogate. Thus, the prediction of the ESM is formulated as:

(10) ${\hat{y}}_{\mathrm{E}\mathrm{S}\mathrm{M}}=\sum _{k=1}^{N_{\mathrm{S}\mathrm{M}}}{w}_k{\hat{y}}_k$(10)

where ${\hat{y}}_k$ is the prediction produced by the kth surrogate model, N_SM is the number of surrogate models that form the ensemble and w_k is the weight attributed to the kth surrogate model, and should be $\sum _{k=1}^{N_{\mathrm{S}\mathrm{M}}}{w}_k=1$. The first ESM approach implemented in this work is based on assigning uniform weights to the surrogates, which means that they contribute equally to the ensemble response. This is known as a simple averaging technique and has been adopted for coastal aquifer management by Sreekanth and Datta (2011b) and Roy and Datta (2017c). The second ensemble approach involves the construction of a weighted-average surrogate model where optimal weights are obtained through optimization. It is noted that the optimization task for the optimal weights must be performed for each one of the constraint functions to form M ensembles of surrogates. The relevant optimization problem is defined as (Zhou et al. 2013, Jiang et al. 2015):

(11) $\begin{array}{l}\min f\left(w\right)=\sqrt{\frac{1}{m}{\displaystyle \sum _{j=1}^m{\left(w_1{y}_1+w_2{y}_2-y_{\text{VDST}}\right)}^2}}\\ \text{s}.\text{t}.{\displaystyle \sum _{k=1}^{N_{SM}}{w}_k}=1\end{array}$(11)

where m is the number of the initial training sample; y₁ and y₂ are the individual predictions of the best two surrogate models for each prediction; and y_VDST is the corresponding response of the VDST model. Thus, the weights for the surrogate models are the decision variables of the above optimization problem, which was solved using the Sequential Quadratic Programming (SQP) method in MATLAB. Since the number of constraint functions may be large, and the construction of the weighted-average surrogate with optimal weights may add considerable computation time, another approach could be to use surrogates, which are known to be accurate for the problem at hand. The reasoning is that we consider experience with the available surrogate models, and the additional computational time is utilized to optimize the weights instead of searching for the best surrogates using a cross-validation strategy. The above multiple surrogate frameworks are sample dependent, which means that using a different training sample may result in a different selection/combination of the best surrogate model. If an ensemble is to be constructed by different instances of the same surrogate model, resampling techniques may be beneficial to increase the diversity of the surrogate model response (e.g. Sreekanth and Datta 2011b, Roy and Datta 2017c). The adaptive multiple SBO framework implemented here is summarized in Figure 4.

Figure 4. Workflow of the multiple SBO frameworks

The adaptive SBO frameworks are compared with the optimal results from the VDST-based optimization, which provided the benchmark solution. Due to the probabilistic nature of the EAS algorithm, 30 independent runs were performed for each of the SBO methods for a more thorough comparison. The VDST-based optimization was run only once to provide the benchmark solution. Table 2 summarizes the different surrogate modelling techniques compared in this study.

Table 2. Summary of the surrogate models (SM) and the abbreviated names given in this work

Abbreviation	Surrogate model
SM1	Linear RBF
SM2	Thin plate spline RBF
SM3	Cubic RBF
SM4	KRG with Gaussian correlation and zero-order polynomial
SM5	KRG with Gaussian correlation and first-order polynomial
SM6	KRG with Gaussian correlation and second-order polynomial
SM7	KRG with exponential correlation and zero-order polynomial
SM8	KRG with exponential correlation and first-order polynomial
SM9	KRG with exponential correlation and second-order polynomial
SM10	Multiple surrogates
SM11	Ensemble (uniform weights)
SM12	Ensemble (optimal weights)

Optimization results

Hypothetical application

The leave-one-out strategy was used to identify the best surrogate model corresponding to each constraint function. Based on the cross-validations scores, the following SMs were identified as best and second best for the total of M = 10 pumping wells and are presented in Table 3.

Table 3. The surrogate models with the lowest RMSE values for each pumping well (PW) – hypothetical application

	PW1	PW2	PW3	PW4	PW5	PW6	PW7	PW8	PW9	PW10
Best	SM5	SM8	SM3	SM2	SM9	SM2	SM3	SM3	SM3	SM3
Second best	SM8	SM5	SM9	SM3	SM6	SM3	SM9	SM2	SM2	SM9

The results of the optimal solutions obtained from all the optimization frameworks are presented in Figure 5 by box-plot visualization, whereas some output statistics are presented in Table 4.

Table 4. Statistics of the best feasible objective function values for the 30 independent runs (best results are in bold)

	Worst	Best	Median	Mean
EAS-VD	4.8666E+03	4.8666E+03	4.8666E+03	4.8666E+03
EAS-SM1	4.3259E+03	4.5216E+03	4.4741E+03	4.4415E+03
EAS-SM2	4.3896E+03	4.8293E+03	4.6207E+03	4.6161E+03
EAS-SM3	4.4811E+03	4.8471E+03	4.7680E+03	4.7498E+03
EAS-SM4	4.6530E+03	4.8332E+03	4.7572E+03	4.7553E+03
EAS-SM5	4.6684E+03	4.8571E+03	4.7802E+03	4.7824E+03
EAS-SM6	4.1570E+03	4.8486E+03	4.7608E+03	4.6920E+03
EAS-SM7	4.6350E+03	4.7178E+03	4.6707E+03	4.6699E+03
EAS-SM8	4.4108E+03	4.7975E+03	4.6948E+03	4.6784E+03
EAS-SM9	4.2503E+03	4.7927E+03	4.5198E+03	4.5284E+03
EAS-SM10	4.6531E+03	4.8471E+03	4.7837E+03	4.7736E+03
EAS-SM11	4.5509E+03	4.8079E+03	4.7774E+03	4.7512E+03
EAS-SM12	4.6569E+03	4.8401E+03	4.7958E+03	4.7863E+03

Figure 5. Performance comparison between direct optimization and SBO frameworks

It is obvious from the results in the previous section that the performance of the surrogate models may vary in several aspects. For example, surrogate models that were identified as best for some constraint functions failed to provide good optimal solutions and a robust performance over the 30 independent runs when they were used alone in a single SBO framework. The box plots in Figure 5 demonstrate that the SBO using multiple surrogates produced good optimal solutions, which compare satisfactorily with the benchmark solution from the EAS-VD model. In addition, the solutions from the ensemble with optimal weights (EAS-SM12) exhibit a low dispersion, and they also have the best median and mean values among all the SBO frameworks (Table 4).

The notches in the box plots of EAS-SM12 and EAS-SM10 appear to overlap more clearly than with the EAS-SM11, indicating that the simple averaging technique may have a different and lower true median from the other two multiple SBO frameworks (Figure. 5). Nevertheless, EAS-SM12 exhibits the most symmetrical distribution and appears as a more robust approach than the other two multiple SBO approaches. Of the nine single surrogates, three (EAS-SM3, EAS-SM4 and EAS-SM5) also produced good optimal solutions and demonstrate a similar performance to the multiple SBO approaches. Interestingly, among all frameworks, the highest maximum and minimum optimal solutions were obtained from SBO with a KRG model of Gaussian correlation and first-order polynomial. Based on the initial training sample, the latter (SM5) was identified as among the best surrogate models only two times (Table 3) compared to other surrogate modelling techniques. This could mean that the size of the initial training sample used to perform cross-validation was not adequate to identify this KRG model as best using other constraint functions. It is not always feasible though, to design large training datasets to find the best surrogates due to the prohibitive computational cost of the VDST simulations. In such cases, the only choice may be to select a surrogate model solely on previous experience with the optimization problem.

In terms of computational savings, the single SBO frameworks required less than 2 hours to converge compared to the approximate time of 2 days for the VDST-based optimization. However, in the case of multiple SBO, additional time was required to apply the cross-validation approach, which took approx. 20 min to calculate RMSE scores for all constraint functions and for all surrogate models. Furthermore, the calculation of optimal weights for the ensembles constructed for each constraint function required an additional total time of 1.3 h. Therefore, whereas EAS-SM12 has probably the most promising performance among the multiple SBO frameworks, it involves a higher computational cost than the single surrogates. Nevertheless, the additional computational time is rather inexpensive in comparison to the direct optimization approach. Overall, of the 12 approaches tested, the results above favour the six SBO approaches: three multiple SBO approaches (EAS-SM10, EAS-SM11 and EAS-SM12) and the single ones based on a cubic RBF (EAS-SM3), a Gaussian KRG with zero-order polynomial (EAS-SM4) and a Gaussian KRG with first-order polynomial (EAS-SM5).

Real-world application

The VDST-based optimization provided the benchmark optimal solution and computational time. Figure 6 presents the salinity distribution simulated by applying the optimal pumping rates obtained from the VDST-based optimization. A single run of the VDST simulation requires approximately 1.8 min. There are 11 main pumping boreholes installed in the aquifer and the simulation output indicates a relatively narrow seaward front.

Figure 6. VDST simulation output, with the locations of the 11 pumping wells. Salinity is expressed in relative concentrations (-)

The SBO framework was run also only once to represent the more realistic case where the computational cost of the VDST simulations does not allow multiple optimization runs. The same optimization problem defined in Equation (6)$minf\left(\mathit{Q}\right)=\left\{\begin{array}{c}-\sum _{i=1}^M{Q}_i\mathrm{i}\mathrm{f}\mathrm{\forall }i=1,2,\dots ,M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)\le c_T\\ M_v\sum _{i=1}^M{\left[max\left(\left(c_i-c_T\right),0\right)\right]}^2\mathrm{i}\mathrm{f}\mathrm{\exists }i=1,2,\dots .M;c_i\left(Q_1,Q_2,\dots ,Q_M\right)>c_T\end{array}\right.(6)$ was solved. Table 5 shows the best surrogate models identified from the cross-validation approach.

Table 5. The surrogate models with the lowest RMSE values for each of the 11 pumping wells (PW) – cross-validation approach

	PW1	PW2	PW3	PW4	PW5	PW6	PW7	PW8	PW9	PW10	PW11
Best	SM3	SM3	SM3	SM3	SM3	SM1	SM6	SM9	SM3	SM1	SM1
Second best	SM3	SM2	SM2	SM2	SM2	SM8	SM9	SM6	SM2	SM5	SM5

Figure 7 presents the distribution of the optimal pumping rates obtained from all optimization frameworks that were tested for the real coastal aquifer problem. The best optimal solution compared to the benchmark solution from the VDST model was produced by the SBO that utilizes a single KRG model with Gaussian correlation and first-order polynomial (EAS-SM5). The second-best was the solution obtained by the cubic RBF model, and the ensemble with optimal weights (EAS-SM12) provided the third-best solution.

Figure 7. Distribution of the optimal pumping rates for all the optimization frameworks. The numbers in parentheses represent the total optimal pumping

In terms of computation times, the direct optimization using the VDST model required approximately 11 days. In contrast, all the SBO approaches required less than 10 h. It is noted that, with the adaptive SBO used in this study, the number of calls of the VDST model may vary for each case and for each optimization run. Therefore, there is no point in generalizing conclusions about the individual computation times for each surrogate model simply based on a single run. It is obvious, however, that the multiple surrogate approaches will add computation time, while the construction of KRG models is more expensive than the RBF models. However, considering that the adaptive SBO framework utilized in this study reduces the required training patterns significantly, the construction of the surrogate models did not contribute much to the overall computation time.

Conclusions

The performance of single and multiple surrogate models was assessed for a hypothetical and a real-world coastal aquifer management problem. The multiple surrogate approach was applied by using either the best or ensembles of surrogates. The application was focused on a single-objective pumping optimization problem where multiple constraint functions must be treated separately.

In both optimization problems, the multiple surrogate approach did not result in a consistent outperformance over the use of single surrogates. However, the weighted-average surrogate (ensemble), based on optimal weights, appeared as the most promising multiple surrogate approach. Notwithstanding the good results from the surrogate ensemble, the single cubic RBF and the KRG models of Gaussian and a first-order polynomial trend performed equally well. The computational gains over the direct optimization were significant, leading to more than 90% reduction for all SBO approaches. However, the adaptive SBO frameworks using multiple surrogates were associated with additional computational cost.

In this work, only the best and the second-best surrogate models were utilized to construct a weighted-average heterogeneous ensemble based on optimal weights. Although using a larger number of surrogate models to form the ensemble may increase the computational cost, it is worthwhile to investigate this approach. Since the efficiency and robustness of surrogate models may vary depending on the optimization problem at hand, further studies should explore other heterogeneous ensembles by using promising surrogate models already used in coastal aquifer management.

Acknowledgments

The authors would like to thank the anonymous reviewers for their suggestions to improve the quality of the paper.

Disclosure statement

No potential conflict of interest was reported by the authors.