## 2 Methodology
### System Modeling Based on Machine Learning (ML) Method-Gaussian Process Regression (GPR)  
A Gaussian process (GP) is a stochastic process, which consists of random variables with normal distributions. Thus, the GP distribution is a joint distribution of any number of random variables (\(\mathbf{X}=[x_{1},\ x_{2},\ \ldots,\ x_{n}]^{T}\)) [34]. From the function-space view, a GP is denoted by a mean function \(\mu(x)\) as well as covariance function \(k(x,x^{\prime})\), and as such, the GP could be expressed in Equation (3) as follows:  
\[f(x)\sim GP\big{(}\mu(x),k\big{(}x,x^{\prime}\big{)}\big{)} \tag{3}\]  
where \(\mu(x)=E[f(x)]\), \(k(x,x^{\prime})=E[(f(x)-\mu(x))(f(x^{\prime})-\mu(x^{\prime}))]\), in which the mean function is usually zero for notational simplicity. The covariance function is considered as a crucial ingredient that defines nearness or similarity. In GPR, the squared exponential (SE) covariance function is probably the most commonly-used kernel, that corresponds to a Bayesian linear regression model. The SE covariance function is expressed as follows:  
\[k\big{(}x,x^{\prime}\big{)}=\sigma_{f}^{2}\exp(-\frac{(x-x^{\prime})^{2}}{2l^{ 2}}) \tag{4}\]  
where the amplitude \(\sigma_{f}\) and characteristic length-scale \(l\) are two hyperparameters to model a smooth rising trend. For the standard linear regression model with Gaussian noise, the Bayesian analysis is written as follows:  
\[\mathbf{y}=f(\mathbf{x})+\varepsilon \tag{5}\]  
where \(\mathbf{x}\) is the input vector, \(\mathbf{y}\) is the observed value considering noise pollution, \(f\) is the objective function, \(f(\mathbf{x})\) is the function value. It can be shown that \(\mathbf{y}\) differs from \(f(\mathbf{x})\) by \(\varepsilon\). Here, suppose that \(\varepsilon\) follows an independent, uniformly Gaussian distribution (\(\mu=0\), \(\sigma_{n}^{2}\)), it has the following distribution:  
\[\varepsilon\sim N(0,\sigma_{n}^{2}) \tag{6}\]  
Subsequently, for multivariate input variables, the covariance function for the standard linear regression model is expressed in Equation (7):  
\[k\Big{(}\mathbf{x},\mathbf{x}^{\prime}\Big{)}=\sigma_{f}^{2}\exp(-\sum_{s=1}^ {n}\frac{(x_{s}-x_{s}^{\prime})^{2}}{2l_{s}^{2}}) \tag{7}\]  
where \(n\) is the dimension of vector \(\mathbf{x}\) (\(x_{s}\), \(x_{s}\)', in which \(s=1,2,\,\ldots,\,n\)), \(l_{s}\) is the characteristic length-scale of the \(s\)-th input variable. Assume a training sample set \(\mathbf{S}=\{(x_{i},y_{i}),i=1,\,2,\,\ldots,t\}\), \(t\) is the number of the training sample. According to the GP described in Equation (3), the function values in the form of a vector is \(\mathbf{F}=[f(x_{1}),f(x_{2}),\,\ldots,f(x_{t})]\), which abides by a joint Gaussian distribution. The prior distribution of \(\mathbf{y}\) can be described in Equation (8):  
\[\mathbf{y}\sim N\Big{(}0,\mathbf{K}(\mathbf{x},\mathbf{x})+\sigma_{t}^{2} \mathbf{I}_{t}\Big{)} \tag{8}\]  
where \(\mathbf{I}_{t}\) is the identity vector of size \(t\), \(\mathbf{K}(\mathbf{x},\mathbf{x})\) donates the \(t\times t\) order symmetric positive definite covariance vector. Subsequently, given a new test input \(\mathbf{x}_{*}\), the test function value \(\mathbf{f}_{*}\) is generated from the joint posterior distribution by the corresponding mean and covariance vector. The joint distribution between the observed values and predictive values can be written as follows:  
\[\left[ \begin{array}{c}\mathbf{y}\\ \mathbf{f}_{*}\end{array} \right]\sim N\bigg{(} \begin{array}{cc}0,&\left[ \begin{array}{c}\mathbf{K}(\mathbf{x},\mathbf{x})+\sigma_{t}^{2}\mathbf{I}_ {t}&\mathbf{K}(\mathbf{x},\mathbf{x}_{*})\\ \mathbf{K}(\mathbf{x}_{*},\mathbf{x})&\mathbf{K}(\mathbf{x}_{*},\mathbf{x}_{* })\end{array}\right]\end{array}