Multi-fidelity surrogate-based optimization for decomposed buffer allocation problems

Abstract

The buffer allocation problem (BAP) for flow lines has been extensively addressed in the literature. In the framework of iterative approaches, algorithms alternate an evaluative method and a generative method. Since an accurate estimation of system performance typically requires high computational effort, an efficient generative method reducing the number of iterations is desirable, for searching for the optimal buffer configuration in a reasonable time. In this work, an iterative optimization algorithm is proposed in which a highly accurate simulation is used as the evaluative method and a surrogate-based optimization is used as the generative method. The surrogate model of the system performance is built to select promising solutions so that an expensive simulation budget is avoided. The performance of the surrogate model is improved with the help of fast but rough estimators obtained with approximated analytical methods. The algorithm is embedded in a problem decomposition framework: several problem portions are solved hierarchically to reduce the solution space and to ease the search of the optimum solution. Further, the paper investigates a jumping strategy for practical application of the approach so that the algorithm response time is reduced. Numerical results are based on balanced and unbalanced flow lines composed of single-machine stations.

Appendix

This section describes briefly how to build the surrogate model using the EKR method. A MATLAB EKR toolbox (both the EKR code and a manual) can be found at Lin et al. (2020) (https://doi.org/10.13140/RG.2.2.16632.19206). More details about the method can be found in Lin et al. (2019).

Given a system configuration $\varvec{x}$ under which the high-fidelity system performance $y_h(\varvec{x})$ is unknown, its low-fidelity outputs $y_{l_j}(\varvec{x}), \forall j$ are firstly corrected by scaling functions. Two scaling functions are considered:

• 1. Additive scaling function: ${\tilde{y}}^{l_j}_i(\varvec{x})=y_{l_j}(\varvec{x})+(y_h(\varvec{x}^0_i)-y_{l_j}(\varvec{x}^0_i)),\forall i\in {\mathcal {N}}, \forall j\in {\mathcal {J}};$

• 2. Multiplicative scaling function: ${\tilde{y}}^{l_j}_i(\varvec{x})=\frac{y_h(\varvec{x}^0_i)}{y_{l_j}(\varvec{x}^0_i)} \cdot y_{l_j}(\varvec{x}),\forall i\in {\mathcal {N}}, \forall j\in {\mathcal {J}}.$

where $y_h(\varvec{x}_i^0)$ and $y_{l_j}(\varvec{x}^0_i)$ are the outputs of the high-fidelity model and the jth low-fidelity model at the ith design point $\varvec{x}^0_i$, respectively. Different scaling functions can be used for different low-fidelity models.

These corrected outputs are expected to have more reliable prediction performance if their scaling functions are estimated by the initial design points close to the unobserved point. Therefore, for the jth low-fidelity model, Kernel regression (Wand and Jones 1995) is used to locally fit a polynomial on the corrected data ${\tilde{y}}^{l_j}_i(\varvec{x})$ with distance-based weights. The system performance estimate at the unobserved point $\varvec{x}$, using the jth low-fidelity model’s corrected data, has a closed form:

$$\begin{aligned} \hat{y}_{l_j}(\varvec{x})=\varvec{e}_1^T(\varvec{X_x}^T\varvec{W_x}\varvec{X_x})^{-1}\varvec{X_x}^T\varvec{W_x}\tilde{\varvec{Y}}_{l_j}, \forall j \in {\mathcal {J}}, \end{aligned}$$

(15)

where $\varvec{e}_1$ is a $(dp+1)$-dimensional vector whose first element is 1 and the rest are 0,

$$\begin{aligned} {\varvec{X}}_{\varvec{x}}=\left[ \begin{matrix} 1 &{}\quad (\varvec{x}^0_1-\varvec{x})^T &{}\quad \cdots &{}\quad \left[ (\varvec{x}^0_1-\varvec{x})^p\right] ^T \\ 1 &{}\quad (\varvec{x}^0_2-\varvec{x})^T &{}\quad \cdots &{}\quad \left[ (\varvec{x}^0_2-\varvec{x})^p\right] ^T \\ \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ 1 &{}\quad (\varvec{x}^0_n-\varvec{x})^T &{}\quad \cdots &{}\quad \left[ (\varvec{x}^0_n-\varvec{x})^p\right] ^T \\ \end{matrix} \right] \end{aligned}$$

is an $n \times (dp+1)$ matrix, p is the degree of the fitted polynomial and $\tilde{\varvec{Y}}_{l_j}=[{\tilde{y}}^{l_j}_1(\varvec{x}),\ldots ,{\tilde{y}}^{l_j}_n(\varvec{x})]^T$. $\varvec{W_x}=\mathrm{diag}\{K_{1,\varvec{\varTheta }_1}(\varvec{x}^0_1-\varvec{x}),\ldots ,K_{1,\varvec{\varTheta }_1}(\varvec{x}^0_n-\varvec{x})\}$ is an $n \times n$ diagonal matrix where

$$\begin{aligned} K_{1,\varvec{\varTheta }_1}(\varvec{x}^0_i-\varvec{x})=\prod _{k=1}^{d}\exp \left\{ -\frac{1}{2\theta _{1,k}}(x_{ik}^0-x_k)^2\right\} , \forall i \in {\mathcal {N}}, \end{aligned}$$

$\varvec{\varTheta }_1=\mathrm{diag}\{\theta _{1,1},\ldots ,\theta _{1,d}\}$, and $\theta _{1,k}>0, k=1,\cdots ,d$ are parameters to be selected.

Finally, the estimates from different low-fidelity models are combined with the weights related to the estimated weighted square error:

$$\begin{aligned} \hat{y}_{\text {EKR}}(\varvec{x})=\sum _{j \in {\mathcal {J}}}w_{l_j}(\varvec{x})\hat{y}_{l_j}(\varvec{x}), \end{aligned}$$

where

$$\begin{aligned} w_{l_j}(\varvec{x})= & {} \frac{K_{2,\theta _2}(\hat{WSE}_{l_j}(\varvec{x}))}{\sum _{i \in {\mathcal {J}}} K_{2,\theta _2}(\hat{WSE}_{l_i}(\varvec{x}))},\\ \hat{WSE}_{l_j}(\varvec{x})= & {} (\mathrm{tr}(\varvec{W_x}))^{-1}\tilde{\varvec{Y}}_{l_j}^T(\varvec{W_x}-\varvec{W_x}^T\varvec{X_x}(\varvec{X_x}^T\varvec{W_x}\varvec{X_x})^{-1}\varvec{X_x}^T\varvec{W_x})\tilde{\varvec{Y}}_{l_j}, \forall j \in {\mathcal {J}}, \end{aligned}$$

and $\mathrm{tr}(\varvec{W_x})$ is the trace of $\varvec{W_x}$. $K_{2,\theta _2}(\cdot )$ has the following form:

$$\begin{aligned} K_{2,\theta _2}(\hat{WSE}_{l_j}(\varvec{x}))=\exp \left\{ {-\frac{\hat{WSE}_{l_j}(\varvec{x})}{2\theta _2WSE_\mathrm{min}(\varvec{x})}}\right\} , \forall j \in {\mathcal {J}}, \end{aligned}$$

where

$$\begin{aligned} WSE_\mathrm{min}(\varvec{x})=\min _{j\in {\mathcal {J}}}\{\hat{WSE}_{l_j}(\varvec{x})\}, \end{aligned}$$

and $\theta _2$ is an unknown parameter to be selected.

The model parameters $\theta _{1,k}, \theta _2$ are selected according to the cross-validation. Where $p=1$, the estimated root square error of $\hat{y}_{\text {EKR}}(\varvec{x})$ as an estimate of the high-fidelity response is provided as:

$$\begin{aligned} \hat{s}_{\text {EKR}}(\varvec{x})=\sqrt{\hat{WSE}(\varvec{x})\left( 1+\frac{1}{2^{d/2}\mathrm{tr}(\varvec{W_x})}\right) } \end{aligned}$$

where

$$\begin{aligned} \hat{WSE}(\varvec{x})=(\mathrm{tr}(\varvec{W_x}))^{-1}\tilde{\varvec{Y}}^T(\varvec{W_x}-\varvec{W_x}^T\varvec{X_x}(\varvec{X_x}^T\varvec{W_x}\varvec{X_x})^{-1}\varvec{X_x}^T\varvec{W_x})\tilde{\varvec{Y}} \end{aligned}$$

and:

$$\begin{aligned} \tilde{\varvec{Y}}=\left[ \sum _{j\in {\mathcal {J}}}w_{l_j}(\varvec{x}){\tilde{y}}^{l_j}_1(\varvec{x}),\ldots ,\sum _{j\in {\mathcal {J}}}w_{l_j}(\varvec{x}){\tilde{y}}^{l_j}_n(\varvec{x})\right] ^T. \end{aligned}$$