## 3 Meta-Modeling Techniques for Surrogating NARX Model Coefficients
### Kriging Meta-Model  
Kriging uses a weighted linear combination of all observed output values to estimate the system output of a random function or random process as a Gaussian process. These weights are generally described by the specified correlation function and are based on the distances between the location to be predicted and the locations already observed. The Kriging meta-model can be mathematically described as [26]:  
\[y(\mathbf{\xi})=\mathbf{\beta}^{T}\mathbf{f}(\mathbf{\xi})+\varepsilon(\mathbf{\xi}),\ \mathbf{\beta}^{T}\mathbf{f}(\mathbf{\xi})=\sum_{i=1}^{D}\beta_{i}f_{i}(\mathbf{\xi}) \tag{14}\]  
where \(\mathbf{\beta}^{T}\mathbf{f}(\mathbf{\xi})\) is the mean value of the Gaussian process; \(f_{i}(\mathbf{\xi})\) is the \(i\)th polynomial basis function with corresponding coefficients \(\beta_{i}\); \(D\) represents the total number of basis function terms which depends on the degree of the polynomials. The residual term \(\varepsilon(\mathbf{\xi})\) is assumed to have zero mean and the following covariance:  
\[\text{Cov}\left[\varepsilon(\mathbf{\xi}_{i}),\ \varepsilon(\mathbf{\xi}_{j}) \right]=\sigma_{Z}^{2}\mathbf{\mathcal{R}}\left(\theta,\ \mathbf{\xi}_{i},\ \mathbf{\xi}_{j}\right) \tag{15}\]  
where \(\mathbf{\xi}_{i}\) and \(\mathbf{\xi}_{j}\) represent the two different sample inputs; \(\sigma_{Z}^{2}\) is the process variance; and \(\mathbf{\mathcal{R}}\left(\theta,\ \mathbf{\xi}_{i},\ \mathbf{\xi}_{j}\right)\) represents the correlation function with parameter \(\theta\). Several auto-correlation functions have been used in the literature [26; 27; 28; 29]. In this study, the Gaussian correlation function is adopted, which can be expressed as:  
\[\mathbf{\mathcal{R}}\left(\xi_{i},\xi_{j}\right)=\prod_{k=1}^{d}\exp[-\theta_{k} (\xi_{i,k}-\xi_{j,k})^{2}] \tag{16}\]  
where the parameters \(\theta_{k}\) are hyper-parameters of auto-correlation functions.  
The Kriging model predicts the response for untriled input \(\mathbf{\xi}_{u}\) through the best linear unbiased predictor:  
\[\hat{y}(\mathbf{\xi}_{u})=\hat{\mathbf{\beta}}^{T}\mathbf{f}(\mathbf{\xi}_{u})+\mathbf{\mathfrak{ R}}^{T}\mathbf{\mathcal{R}}^{-1}\left(\hat{\mathbf{\beta}}^{T}\mathbf{f}(\mathbf{\xi}_{u})\right) \tag{17}\]  
\[\sigma_{\mathbf{\hat{M}}}^{2}(\mathbf{\xi}_{u})=\sigma_{Z}^{2}\left(1-\left[\mathbf{f}^{T} (\xi_{u})\ \mathfrak{R}^{T}(\mathbf{\xi}_{u})\right]\begin{bmatrix}\mathbf{0}&\mathbf{\mathcal{L}}\\ \mathbf{\mathcal{L}}&\mathbf{\mathcal{R}}\end{bmatrix}^{-1}\begin{bmatrix}f(\mathbf{\xi}_{ u})\\ \mathfrak{R}(\mathbf{\xi}_{u})\end{bmatrix}\right) \tag{18}\]where \(\mathbf{\mathfrak{R}}(\mathbf{\xi}_{u})=[\mathbf{\mathcal{R}}(\mathbf{\xi}_{u},\ \mathbf{\xi}_{1}),\ \mathbf{ \mathcal{R}}(\mathbf{\xi}_{u},\ \mathbf{\xi}_{2}),\ \ldots,\ \mathbf{\mathcal{R}}(\mathbf{\xi}_{u},\ \mathbf{\xi}_{N})]^{T}\) defines the correlation matrix between untried sample point \(\mathbf{\xi}_{u}\) and training samples; \(\mathbf{\mathcal{L}}\) is the matrix of basis functions at training points of dimension \(N\times D\); and \(\mathbf{\mathcal{R}}\) is the correlation matrix for training samples.  
The model parameters \(\mathbf{\hat{\beta}}\) and \(\hat{\sigma}_{Z}^{2}\), and the hyper-parameters \(\theta\) can be determined using the maximum-likelihood estimation (MLE) [27; 28]:  
\[\mathbf{\hat{\beta}}=\big{(}\mathbf{\mathcal{L}}^{T}\mathbf{\mathcal{R}}^{-1}\mathbf{\mathcal{ L}}\big{)}^{-1}\mathbf{\mathcal{L}}^{T}\mathbf{\mathcal{R}}^{-1}\mathbf{Y} \tag{19}\]  
\[\hat{\sigma}_{Z}^{2}=\frac{1}{N}\Big{(}\mathbf{Y}-\mathbf{\mathcal{L}}\mathbf{\hat{\beta} }\Big{)}^{T}\mathbf{\mathcal{R}}^{-1}\Big{(}\mathbf{Y}-\mathbf{\mathcal{L}}\mathbf{\hat{\beta }}\Big{)} \tag{20}\]  
More details of the Kriging meta-model can be referred to in the work of [PERSON] et al. [28]. The MATLAB-based toolbox DACE [30] is used in this study for the Kriging model.