## 2 Materials and Methods
### Gradient-Enhanced Kriging  
If the objective function values and sensitivity information of the corresponding points are available, incorporating gradient information into the construction of the kriging model can significantly improve the accuracy of fitting [27]. Figure 3 presents an example of fitting the two-dimensional Griewank function using both OK and GEK, with the same number of training points. Clearly, the latter achieves higher fitting accuracy. The GEK model is an  
Figure 2: The variation trend of spatial correlation between sample points under different hyperparameter values (one-dimensional correlation function). extension of the OK model. To fit the GEK model, the correlation terms are also extended to include the gradient terms. First, the correlation matrix is expanded as follows:  
\[\mathbf{R}= \begin{bmatrix}\mathbf{Q^{1,1}}&\mathbf{Q^{1,2}}\\ \left(\mathbf{Q^{1,2}}\right)^{T}&\mathbf{Q^{2,2}}\end{bmatrix}  \tag{12}\]  
where \(\mathbf{Q}^{1,1}\) is the correlation matrix that does not include gradient information:  
\[\mathbf{Q}^{1,1}= \begin{bmatrix}\psi\left(x^{1},x^{1}\right)&\ldots&\psi\left(x^{ 1},x^{n}\right)\\ \vdots&\ddots&\vdots\\ \psi\left(x^{n},x^{1}\right)&\ldots&\psi\left(x^{n},x^{n}\right)\end{bmatrix}  \tag{13}\]  
\(\mathbf{Q}^{1,2}\) is the extended correlation matrix that includes the first-order gradient information of \(\mathbf{R}\):  
\[\mathbf{Q}^{1,2}= \begin{bmatrix}\frac{\partial\psi\left(x^{1},x^{1}\right)}{ \partial x^{1}\partial x^{1}}&\ldots&\frac{\partial\psi\left(x^{1},x^{n} \right)}{\partial x^{1}\partial x^{n}}\\ \vdots&\ddots&\vdots\\ \frac{\partial\psi\left(x^{n},x^{1}\right)}{\partial x^{1}}&\ldots&\frac{ \partial\psi\left(x^{n},x^{n}\right)}{\partial x^{n}\partial x^{n}}\end{bmatrix}  \tag{14}\]  
\(\mathbf{Q}^{2,2}\) is the extended correlation matrix that includes the second-order gradient information of \(\mathbf{R}\):  
\[\mathbf{Q}^{2,2}= \begin{bmatrix}\frac{\partial^{2}\psi\left(x^{1},x^{1}\right)}{ \partial x^{1}\partial x^{1}}&\ldots&\frac{\partial^{2}\psi\left(x^{1},x^{n} \right)}{\partial x^{1}\partial x^{n}}\\ \vdots&\ddots&\vdots\\ \frac{\partial^{2}\psi\left(x^{n},x^{1}\right)}{\partial x^{n}\partial x^{1} }&\ldots&\frac{\partial^{2}\psi\left(x^{n},x^{n}\right)}{\partial x^{n} \partial x^{n}}\end{bmatrix}  \tag{15}\]  
The spatial correlation between the prediction point and the training points is also extended in the same manner as follows:  
\[\mathbf{r}= \begin{bmatrix}\psi\left(\mathbf{x}^{n+1},\mathbf{x}^{1}\right),\ldots,\psi \left(\mathbf{x}^{n+1},\mathbf{x}^{n}\right),\frac{\partial\psi\left(\mathbf{x}^{n+1},\bm {x}^{1}\right)}{\partial x^{1}},\ldots,\frac{\partial\psi\left(\mathbf{x}^{n+1}, \mathbf{x}^{n}\right)}{\partial x^{n}}\end{bmatrix}  \tag{16}\]  
Similarly, the extension of the function value vector also includes the gradient information of the response:  
\[\mathbf{y}= \begin{bmatrix}f\left(\mathbf{x}^{1}\right),\ldots f(\mathbf{x}^{n}),\frac{ \partial f(\mathbf{x}^{1})}{\partial\mathbf{x}^{1}},\ldots,\frac{\partial f(\mathbf{x}^{ n})}{\partial\mathbf{x}^{n}}\end{bmatrix}  \tag{17}\]  
Figure 3: Comparison of fitting accuracy between gradient-enhanced kriging (GEK) and ordinary kriging (OK). (**a**) Real two-dimensional Griewank function, (**b**) OK fitting, and (**c**) GEK fitting.