梯度增强 Kriging 代理模型理论相关推导

梯度信息可用于提高 Kriging 模型精度，而如果采用Adjoint方法等快速求解梯度方法，还可提高建立Kriging模型的效率。利用梯度信息来提高Kriging模型的精度，成为一种新的代理模型方法，称为梯度增强型 Kriging（Gradient-Enhanced Kriging, GEK）模型。

4.1 问题定义

除一阶导数信息外，二阶导数信息也可用于提高代理模型的精度。而由于获得二阶导数的计算代价较大，且对于提高代理模型精度的作用有限，因而本文主要介绍在Kriging模型中引入一阶偏导数信息的情况。

对于一个具有 $m$ 个设计变量的优化问题，首先假设对目标函数（或状态变量）$y$ 在$n$ 个抽样位置，获得 $n$ 个函数值及 $n\times m$ 个偏导数值，则获得设计空间抽样位置及函数响应值分别为：
$$\begin{array}{rl}& \mathbf{S}={\left[\begin{array}{ccccccccc}{\mathbf{x}}^{(1)}& \dots & {\mathbf{x}}^{(n)}& {\mathbf{x}}^{(1)}& \dots & {\mathbf{x}}^{(1)}& \dots & {\mathbf{x}}^{(n)}\dots & {\mathbf{x}}^{(n)}\end{array}\right]}^T\in {\mathbb{R}}^{(n+nm)\times m}\\ & {\mathbf{y}}_s={\left[\begin{array}{cccccccccc}{y}^{(1)}& \dots & y^{(n)}& \frac{\partial y^{(1)}}{\partial x_1}& \dots & \frac{\partial y^{(1)}}{\partial x_n}& \dots & \frac{\partial y^{(n)}}{\partial x_1}& \dots & \frac{\partial y^{(n)}}{\partial x_m}\end{array}\right]}^T\in {\mathbb{R}}^{n+nm}\end{array}\text{\hspace{1em}(1)}$$
由上式，每一个偏导数信息都被看作一个独立的样本信息。如果在某些点处的样本信息不可用，则从上述样本数据集的相应位置上去除即可；如果没有任何偏导数信息，GEK模型将退化为传统的Kriging模型。

4.2 GEK模型的建立

GEK模型对未知函数的预估值定义为所有抽样函数值和偏导数值的线性加权，即：
$$\hat{y}(\mathbf{x})=\sum _{i=1}^n{\omega }^{(i)}{y}^{(i)}+\sum _{j=1}^m\sum _{i=1}^n{\lambda }_j^{(i)}\frac{\partial y^{(i)}}{\partial x_j}\text{\hspace{1em}(2)}$$
上式中： ${\omega }^{(i)}$ 为第 $i$ 个（抽样位置）函数值的加权系数；${\lambda }_j^{(i)}$ 为 $\partial y^{(i)}/\partial x_j$ （第 $i$ 个抽样位置处函数对 $j$ 维设计变量的偏导数）的加权系数。如果未知函数为物理量，则 ${\lambda }_j^{(i)}$ 的量纲与 $\omega \Delta x_k$ 相同。与传统的 Kriging 模型类似，引入静态随机过程假设：
$$Y(\mathbf{x})={\beta }_0+Z(\mathbf{x})\text{\hspace{1em}(3)}$$
上式中，${\beta }_0$ 为未知常数，也被称为全局趋势模型，代表 $Y(\mathbf{x})$的数学期望值；$Z(\cdot )$ 是均值为零、方差为 ${\sigma }^2$ 的静态随机过程。在设计空间不同位置处，对应的随机变量之间的协方差满足：
$$\begin{array}{rl}& Cov[Z({\mathbf{x}}^{(i)}),Z({\mathbf{x}}^{(j)})]={\sigma }^{2}R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)})\\ & Cov[\frac{\partial Z({\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_k^{(i)}},Z({\mathbf{x}}^{(j)})]={\sigma }^2\frac{\partial R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)})}{\partial {\mathbf{x}}_k^{(i)}}\\ & Cov[Z({\mathbf{x}}^{(i)}),\frac{\partial Z({\mathbf{x}}^{(j)})}{\partial {\mathbf{x}}_l^{(j)}}]={\sigma }^2\frac{\partial R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)})}{\partial {\mathbf{x}}_k^{(j)}}\\ & Cov[\frac{\partial Z({\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_l^{(i)}},\frac{\partial Z({\mathbf{x}}^{(j)})}{\partial {\mathbf{x}}_k^{(j)}}]={\sigma }^2\frac{{\partial }^{2}R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(j)})}{\partial {\mathbf{x}}_l^{(i)}\partial {\mathbf{x}}_k^{(j)}}\end{array}\text{\hspace{1em}(4)}$$
GEK模型寻找 $\mathbf{\omega }$ 和 $\lambda$ 使得下述均方差最小：
$$MSE[\hat{y}(\mathbf{x})]=E\left[{\left(\sum _{i=1}^n{\omega }^{(i)}{y}^{(i)}+\sum _{j=1}^n\sum _{i=1}^n{\lambda }_j^{(i)}\frac{\partial y^{(i)}}{\partial x_j}-Y(\mathbf{x})\right)}^2\right]\text{\hspace{1em}(5)}$$
并满足无偏估计条件：
$$E\left[Y(\mathbf{x})\right]=E\left[\sum _{i=1}^n{\omega }^{(i)}{y}^{(i)}+\sum _{j=1}^n\sum _{i=1}^n{\lambda }_j^{(i)}\frac{\partial y^{(i)}}{\partial x_j}-Y(\mathbf{x})\right]\text{\hspace{1em}(6)}$$
无偏估计条件可简化如下：
$$\sum _{i=1}^n{\omega }^{(i)}=1\text{\hspace{1em}(7)}$$
构建拉格朗日函数：
$$L=E\left[{\left(\sum _{i=1}^n{\omega }^{(i)}+\sum _{j=1}^n\sum _{i=1}^n{\lambda }_j^{(i)}\frac{\partial y^{(i)}}{\partial x_j}-Y(\mathbf{x})\right)}^2\right]+\mu \left[\sum _{i=1}^n{\omega }^{(i)}-1\right]\text{\hspace{1em}(8)}$$
令：
$$\{\begin{array}{l}\frac{\partial L}{\partial {\omega }^{(k)}}=E\left[2\left({\sum }_{i=1}^n{\omega }^{(i)}{y}^{(i)}+{\sum }_{j=1}^n{\sum }_{i=1}^n{\lambda }_j^{(i)}\frac{\partial y^{(i)}}{\partial x_j}-Y(\mathbf{x})\right)y({\mathbf{x}}^{(k)})\right]+\mu =0\\ \frac{\partial L}{\partial {\lambda }_l^{(k)}}=E\left[2\left({\sum }_{i=1}^n{\omega }^{(i)}{y}^{(i)}+{\sum }_{j=1}^n{\sum }_{i=1}^n{\lambda }_j^{(i)}-Y(\mathbf{x})\right)\frac{\partial y^{(k)}}{\partial x_l^{(k)}}\right]=0\\ \frac{\partial L}{\partial \mu }={\sum }_{i=1}^n{\omega }^{(i)}-1=0\end{array}\text{\hspace{1em}(9)}$$
又：
$$\begin{array}{rl}& E[Y({\mathbf{x}}^{(i)})\cdot Y({\mathbf{x}}^{(k)})]={\beta }_0^2+{\sigma }^{2}R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})\\ & E\left[\frac{\partial Y({\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_j^{(i)}}\cdot Y({\mathbf{x}}^{(k)})\right]={\sigma }^2\frac{\partial R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_j^{(i)}}\\ & E\left[Y(\mathbf{x})\cdot Y({\mathbf{x}}^{(k)})\right]={\beta }_0^2+{\sigma }^{2}R({\mathbf{x}}^{(k)},\mathbf{x})\\ & E\left[Y({\mathbf{x}}^{(i)})\cdot \frac{\partial Y({\mathbf{x}}^{(k)})}{\partial {\mathbf{x}}_l^{(k)}}\right]={\beta }_0^2+{\sigma }^2\frac{\partial R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(k)})}{\partial {\mathbf{x}}_l^{(k)}}\\ & E\left[\frac{\partial Y({\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_j^{(i)}}\cdot \frac{\partial Y({\mathbf{x}}^{(k)})}{\partial {\mathbf{x}}_l^{(k)}}\right]={\beta }_0^2+{\sigma }^2\frac{{\partial }^{2}R({\mathbf{x}}^{(i)},{\mathbf{x}}^{(k)})}{\partial {\mathbf{x}}_j^{(i)}\partial {\mathbf{x}}_l^{(k)}}\\ & E\left[Y(\mathbf{x})\cdot \frac{\partial Y({\mathbf{x}}^{(k)})}{\partial {\mathbf{x}}_l^{(k)}}\right]={\beta }_0^2+{\sigma }^2\frac{R({\mathbf{x}}^{(k)},\mathbf{x})}{\partial {\mathbf{x}}_l^{(k)}}\end{array}\text{\hspace{1em}(10)}$$
代入式 (8) 可得：
$$\{\begin{array}{l}{\sum }_{i=1}^n{\omega }^{(i)}R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})+{\sum }_{j=1}^n{\sum }_{i=1}^n{\lambda }_j^{(i)}\frac{\partial R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_j^{(i)}}+\frac{\mu }{2{\sigma }^2}=R({\mathbf{x}}^{(k)},\mathbf{x})\\ {\sum }_{i=1}^n{\omega }^{(i)}{y}^{(i)}\frac{\partial R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_l^{(k)}}+{\sum }_{j=1}^n{\sum }_{i=1}^n{\lambda }_j^{(i)}\frac{{\partial }^{2}R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_l^{(k)}\partial {\mathbf{x}}_j^{(i)}}=\frac{\partial R({\mathbf{x}}^{(k)},{\mathbf{x}}^{(i)})}{\partial {\mathbf{x}}_l^{(k)}}\\ {\sum }_{i=1}^n{\omega }^{(i)}=1\end{array}\text{\hspace{1em}(11)}$$
将上式改写为矩阵形式：
$$\left[\begin{array}{cc}\overline{\mathbf{R}}& \overline{\mathbf{F}}\\ {\overline{\mathbf{F}}}^T& 0\end{array}\right]\left[\begin{array}{c}\mathbf{\lambda }\\ \tilde{\mu }\end{array}\right]=\left[\begin{array}{c}\overline{\mathbf{r}}\\ 1\end{array}\right]\text{\hspace{1em}(12)}$$
上式中：
$$\overline{\mathbf{R}}=\left[\begin{array}{cc}\mathbf{R}& \partial \mathbf{R}\\ \partial {\mathbf{R}}^T& {\partial }^2\mathbf{R}\end{array}\right]\text{\hspace{1em}(13)}$$
其中，$\mathbf{R}$ 与 Kriging 模型中一致：
$$\begin{array}{rl}& \mathbf{R}=\left[\begin{array}{ccc}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})& \dots & R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(n)})\\ ⋮& \ddots & ⋮\\ R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})& \dots & R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})\end{array}\right]\\ & \partial \mathbf{R}=\left[\begin{array}{ccccccc}\frac{\partial R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_1^{(1)}}& \dots & \frac{\partial R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_m^{(1)}}& \dots & \frac{\partial R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{\partial R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_1^{(n)}}\\ ⋮& \ddots & ⋮& \ddots & ⋮& \ddots & ⋮\\ \frac{\partial R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{\partial R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_m^{(n)}}& \dots & \frac{\partial R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{\partial R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_m^{(n)}}\end{array}\right]\\ & {\partial }^2\mathbf{R}=\left[\begin{array}{ccccccc}\frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{\partial ({\mathbf{x}}_1^{(1)}{)}^2}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_1^{(1)}\partial {\mathbf{x}}_m^{(1)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_1^{(1)}\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_1^{(1)}{\mathbf{x}}_m^{(n)}}\\ ⋮& \ddots & ⋮& \ddots & ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_m^{(1)}\partial {\mathbf{x}}_1^{(1)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(1)})}{(\partial {\mathbf{x}}_m^{(1)}{)}^2}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_m^{(1)}\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(1)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_m^{(1)}\partial {\mathbf{x}}_m^{(n)}}\\ ⋮& \ddots & ⋮& \ddots & ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_1^{(n)}\partial {\mathbf{x}}_1^{(1)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_1^{(n)}\partial {\mathbf{x}}_m^{(1)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})}{(\partial {\mathbf{x}}_1^{(n)}{)}^2}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_1^{(n)}\partial {\mathbf{x}}_m^{(n)}}\\ ⋮& \ddots & ⋮& \ddots & ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_m^{(n)}\partial {\mathbf{x}}_1^{(1)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(1)})}{\partial {\mathbf{x}}_m^{(n)}\partial {\mathbf{x}}_m^{(1)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})}{\partial {\mathbf{x}}_m^{(n)}\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{{\partial }^{2}R({\mathbf{x}}^{(n)},{\mathbf{x}}^{(n)})}{(\partial {\mathbf{x}}_m^{(n)}{)}^2}\end{array}\right]\end{array}\text{\hspace{1em}(14)}$$
式 (11) 中：
$$\overline{\mathbf{r}}=\left[\begin{array}{c}\mathbf{r}\\ \partial \mathbf{r}\end{array}\right]\text{\hspace{1em}(15)}$$
其中，$\mathbf{r}$ 与 Kriging 模型中一致：
$$\begin{array}{rl}& \mathbf{r}={\left[\begin{array}{ccc}R({\mathbf{x}}^{(1)},\mathbf{x})& \dots & R({\mathbf{x}}^{(n)},\mathbf{x})\end{array}\right]}^T\\ & \partial \mathbf{r}={\left[\begin{array}{ccccccc}\frac{\partial R({\mathbf{x}}^{(1)},\mathbf{x})}{\partial {\mathbf{x}}_1^{(1)}}& \dots & \frac{\partial R({\mathbf{x}}^{(1)},\mathbf{x})}{\partial {\mathbf{x}}_m^{(1)}}& \dots & \frac{\partial R({\mathbf{x}}^{(n)},\mathbf{x})}{\partial {\mathbf{x}}_1^{(n)}}& \dots & \frac{\partial R({\mathbf{x}}^{(n)},\mathbf{x})}{\partial {\mathbf{x}}_m^{(n)}}\end{array}\right]}^T\end{array}\text{\hspace{1em}(16)}$$
式 (11) 中：
$$\begin{array}{rl}& \mathbf{\lambda }={\left[\begin{array}{ccccccccc}{\omega }^{(1)}& \dots & {\omega }^{(n)}& {\lambda }_1^{(1)}& \dots & {\lambda }_m^{(1)}& {\lambda }_1^{(n)}& \dots & {\lambda }_m^{(n)}\end{array}\right]}^T\\ & \overline{\mathbf{F}}=[\underset{n}{\underbrace{1\dots 1}}\underset{nm}{\underbrace{0\dots 0}}{]}^T\in {\mathbb{R}}^{n+nm}\end{array}\text{\hspace{1em}(17)}$$
综上得到 GEK 模型的预测值为：
$$\hat{y}(\mathbf{x})={\left[\begin{array}{c}\overline{\mathbf{r}}\\ 1\end{array}\right]}^T\left[\begin{array}{cc}\overline{\mathbf{R}}& \overline{\mathbf{F}}\\ {\overline{\mathbf{F}}}^T& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{y}}_s\\ 0\end{array}\right]={\beta }_0+{\overline{\mathbf{r}}}^T{\overline{\mathbf{R}}}^{-1}({\mathbf{y}}_s-{\beta }_0\overline{\mathbf{F}})\text{\hspace{1em}(18)}$$
上式中：
$${\beta }_0=({\overline{\mathbf{F}}}^T{\overline{\mathbf{R}}}^{-1}\overline{\mathbf{F}}{)}^{-1}{\overline{\mathbf{F}}}^T{\overline{\mathbf{R}}}^{-1}{\mathbf{y}}_s\text{\hspace{1em}(19)}$$
预测的均方误差为：
$$MSE[\hat{y}(\mathbf{x})]={\sigma }^2\left(1-{\left[\begin{array}{c}\overline{\mathbf{r}}\\ 1\end{array}\right]}^T\left[\begin{array}{cc}\overline{\mathbf{R}}& \overline{\mathbf{F}}\\ {\overline{\mathbf{F}}}^T& 0\end{array}\right]\left[\begin{array}{c}\overline{\mathbf{r}}\\ 1\end{array}\right]\right)={\sigma }^2\left(1-{\overline{\mathbf{r}}}^T{\overline{\mathbf{R}}}^{-1}\overline{\mathbf{r}}+\frac{(1-{\overline{\mathbf{F}}}^T{\overline{\mathbf{R}}}^{-1}\overline{\mathbf{r}}{)}^2}{{\overline{\mathbf{F}}}^T{\overline{\mathbf{R}}}^{-1}\overline{\mathbf{F}}}\right)\text{\hspace{1em}(20)}$$

梯度增强 Kriging 模型的相关函数，模型训练过程与Kriging模型十分相似，可以参考之前的博客文章(链接: link)，在此不再赘述。

参考文献

[1] HAN Z H, Görtz S, Zimmermann R. Improving Variable-fidelity Surrogate Modeling via Gradient-enhanced Kriging and a Generalized Hybrid Bridge Function[J]. Aerospace Science and Technology, 2013, 25(1):177-189.