Adaptive in situ model refinement for surrogate-augmented population-based optimization

Abstract

In surrogate-based optimization (SBO), the deception issues associated with the low fidelity of the surrogate model can be dealt with in situ model refinement that uses infill points during optimization. However, there is a lack of model refinement methods that are both independent of the choice of surrogate model (neural networks, radial basis functions, Kriging, etc.) and provides a methodical approach to preserve the fidelity of the search dynamics, especially in the case of population-based heuristic optimization processes. This paper presents an adaptive model refinement (AMR) approach to fill this important gap. Therein, the question of when to refine the surrogate model is answered by a novel hypothesis testing concept that compares the distribution of model error and distribution of function improvement over iterations. These distributions are respectively computed via a probabilistic cross-validation approach and by leveraging the probabilistic improvement information uniquely afforded by population-based algorithms such as particle swarm optimization. Moreover, the AMR method identifies the size of the batch of infill points needed for refinement. Numerical experiments performed on multiple benchmark functions and an optimal (building energy) planning problem demonstrate AMR’s ability to preserve computational efficiency of the SBO process while providing solutions of more attractive fidelity than those provisioned by a standard SBO approach.

Appendices

Appendix A: Kernel density estimation (KDE)

KDE is a non-parametric model to estimate the probability density function of random variables. Here, it is assumed that $\varDelta f = (\varDelta f_{1}, \varDelta f_{2}, ..., \varDelta f_{N_{\text {pop}}})$ is an independent and identically distributed sample drawn from a distribution with an unknown density Θ_Δf. The kernel density estimator can then be used to determine Θ_Δf, as given by:

$$ \hat{\varTheta}_{\varDelta f}(x; H) = \frac{1}{N_{\text{pop}}} \sum\limits_{i=1}^{N_{\text{pop}}} K_{H}(x-x_{i}) $$

(12)

Here, the kernel K(x) is a symmetric probability density function, H is the bandwidth matrix which is symmetric and positive-definite, and K_H(x) = |H|^− 1/2K(H^− 1/2x). The choice of K is not as crucial as the choice of the H estimator for the accuracy of the KDE (Epanechnikov 1969). In this article, we consider $K(x) = (2\pi )^{-d/2}\exp (-0.5x^{T}x)$, the standard normal throughout. The mean integrated squared error (MISE) method is used as a criterion for selecting the bandwidth matrix H (Duong and Hazelton 2003) as follows:

$$ \text{MISE}(H) = \mathbb{E}\Big[\int\Big(\varTheta_{\varDelta f}(x; H)-\hat{\varTheta}_{\varDelta f}(x; H)\Big)^{2} \Big] $$

(13)

Appendix B: Particle Swarm Optimization

Particle swarm optimization (PSO) is a population-based optimization method introduced by Kennedy and Eberhart (1995). In this method, each particle’s movement is described in terms of its velocity (v_i(t)) and its location (x_i(t)), where i denotes the i th particle and t denotes the t th iteration. Here, we specifically exploit the MDPSO algorithm developed by Chowdhury et al. (2013), which includes explicit diversity preservation in addition to the standard PSO dynamics, in order to provide greater robustness. In MDPSO, the velocity and location of particles are updated as follows:

$$ \mathbf{x}_{i}(t+1) = \mathbf{x}_{i}(t) + \mathbf{v}_{i}(t+1) $$

(14)

$$ \begin{array}{@{}rcl@{}} \mathbf{v}_{i}(t+1) &=& \omega \mathbf{v}_{i}(t) + r_{1}C_{1}(\mathbf{P}_{i}^{l}(t)-\mathbf{x}_{i}(t))\\ &&+ r_{2}C_{2}(\mathbf{P}^{g}(t)-\mathbf{x}_{i}(t))+r_{3}\gamma_{c}\hat{\nu}_{i}(t) \end{array} $$

(15)

Here, x_i(t) and v_i(t) respectively denote the position and the velocity of particle i at the t^th iteration; ω, C₁, and C₂ represent the inertial weight, the individual search, and the global search coefficients, respectively; these are used to balance the local search (exploitation) and the global search (exploration); $\mathbf {P}_{i}^{l}(t)$ is the local leader of particle i at the t^th iteration, which represents the best local solution found in the motion history of particle i; P^g(t) is the global leader of the entire swarm at the t^th iteration, which is determined by comparing the local leaders of all particles; γ_c is the coefficient used to weigh the explicit diversity preservation component; $\hat {\nu }_{i}(t)$ is the explicit diversity preservation vector; and r₁, r₂, and r₃ are random real numbers between 0 and 1.

Appendix C: Predictive estimation of model fidelity

Predictive estimation of model fidelity (PEMF) method (Mehmani et al. 2015) can be perceived as a novel sequential implementation of k-fold cross-validation, with carefully constructed error measures that are significantly less sensitive to outliers and the DoE (compared with mean or root mean square error measures). The PEMF method predicts the error by capturing the variation of the surrogate model error with an increasing density of training points (without investing any additional test points). Algorithm 1 summarizes the PEMF method.

Appendix D: The settings of COSMOS and MDPSO

Table 5 The COSMOS and MDPSO settings for the analytical problem and application problem

Appendix E: Description of six-hump camel back function

Figure 12 shows the six-hump camel back function, a well-known 2D test function used for benchmarking global optimization and surrogate modeling methods (Molga and Smutnicki 2005). It has six local minima, with two of them being the global minima (Molga and Smutnicki 2005). The global minima are located at (− 0.0898, 0.7126) and (0.0898,− 0.7126), and gives a minimum function value of f(x^∗) = − 1.0316.

Fig. 12

The response surface of the six-hump camel back function

$$ f(x_{1}, x_{2}) = \left( \frac{{x_{1}^{4}}}{3} - 2.1 {x_{1}^{2}} + 4\right){x_{1}^{2}} + x_{1} x_{2} + 4\left( {x_{2}^{2}} - 1\right){x_{2}^{2}} $$

(16)

where − 3 ≤ x₁ ≤ 3 and − 2 ≤ x₂ ≤ 2.

Appendix F: Further performance illustration on the six-hump problem

The following illustrations in Fig. 13 respectively represent the best individual run of PSO-AMR, PSO-AMR-Local, PSO-SBO-k1, and PSO-SBO-k2.

Fig. 13

The selected best run of each optimization approach, for the six-hump camel back function. Left figures show filled contours and dashed line contours, respectively, representing the true function and the surrogate model (SM) approximation; in the PSO-AMR and PSO-AMR-Local cases, there are two SM contours, corresponding to the models before and after refinement. Right figures show the convergence and SM error histories over the optimization processes