CN105893646A - 一种基于增量Kriging的序列优化实验设计方法 - Google Patents

Abstract

本发明公开了一种基于增量Kriging的序列优化实验设计方法，其包括步骤1、初始实验设计的选择；步骤2、初始建模以及模型验证；步骤3、增量的Kriging建模(IKM)；步骤4、优化采样准则：步骤5、更新准则：6：DACE建模；Kriging方法的建模是有效地，随着采样点的增加(对于二维问题，当采样点数目大于600)，Kriging的建模时间陡然增加，因此这里提出一种改进的Kriging方法，在采样点数目非常大的情况下最大限度地降低时间消耗，提出的基于增量Kriging的序列优化实验设计方法在减少建模时间上将会在工程仿真中得到广泛的应用。

Description

一种基于增量Kriging的序列优化实验设计方法

技术领域

本发明涉及一种序列优化设计方法来提高工程应用中的建模效率。

背景技术

(1).现代机电产品日益复杂，对其进行建模分析(特别是仿真模型估值)需要耗费大量的计算时间。尽管计算机的计算速度不断提高，但对于越来越复杂的仿真模型，如有限元分析、计算流体动力学等计算密集型的仿真模型需大量计算时间，因此不能满足工业界对仿真分析的需求。为应对巨大的挑战，在过去的几十年中，Kriging方法应运而生并在工业界得到了普遍应用。该方法能够在不影响仿真模型精度的情况下最大程度地减少优化迭代过程中源模型的仿真次数，从而降低对计算资源的消耗。

(2).Kriging方法是一种通过已知点来预测位置观察点的一种插值方法。Kriging方法利用方差的变化来表达空间的变化，而且可以保证有空间分布得到的预测值的误差最小。Kriging方法源于南非的一位矿业工程师Daniel Gerhard Krige。他率先将统计学应用于地质，矿业的分析与评估。20世纪70年代，法国的数学家Georges Matheron对Krige D.G.的研究成果进行系统化、理论化分析，进而提出了一种插值和外推的理论。随着计算机技术的飞速发展，该方法又被运用在计算科学，产生Kriging模型，该模型作为一种代理模型可以大大提高计算速度。后来人们将试验设计方法与使用Kriging模型的全过程称为计算机试验设计与分析(DACE)，被广泛应用于多个行业，包括采矿业、水文地质学、自然资源、环境科学、遥感、工程分析、机电产品设计和计算机实验中的黑箱模型。

(3).当采样点比较少的时候，Kriging方法的建模是有效地，随着采样点的增加(对于二维问题，当采样点数目大于600)，Kriging的建模时间陡然增加，因此这里提出一种改进的Kriging方法，在采样点数目非常大的情况下最大限度地降低时间消耗。此外，尽管一次实验设计是非常直观的，但是利用它构造Kriging模型是不合适的。首先，对于大量的仿真估值，一次实验设计是非常耗时的；其次，当采样大量的数据点时，它或许能够导致模型的不正常使用。因此，序列实验设计或许是一种更好的选择。与一次实验设计相比，序列实验设计能控制整个采样过程并充分利用先前的模型信息来决定下面需要采样的有效点，所以，它是一种稳定、有效而又精确的实验设计方法，并在许多工程设计中得到了广泛的应用。

(4).以上述为背景，所提出的基于增量Kriging的序列优化实验设计方法在减少建模时间上将会在工程仿真中得到广泛的应用。

所以，需要改经。

发明内容

本发明所要解决的技术问题是提供一种新的基于增量Kriging的序列优化实验设计方法。

本发明的技术方案如下：一种基于增量Kriging的序列优化实验设计方法，该方法包括步骤如下：

步骤1、初始实验设计的的选择：

初始实验设计采用拉丁方空间填充采样(LHS)方法，在整个设计空间均匀统一地获取10n+1(n为模型的维度)个初始采样点，并获取相应的仿估值；

步骤2、初始建模以及模型验证：

根据初始采样点和相应的估值，并利用计算机试验设计与分析(DACE)来建立初始的Kriging模型，当获得Kriging近似模型之后，需要进行模型精度的验证，交叉验证法是的验证方法，这里采用了一种高效而又简单的改进留一交叉验证方法实现Kriging模型精度的验证，主要利用当前所有采样点，经过N-1次(N为当前的总样本数)建模和验证之后来决定Kriging模型是否需要继续进行更新，这种验 证方法是当前最行之有效且节约时间的模型精度验证方法；

步骤3、增量的Kriging建模(IKM)

IKM方法是在参数θ保持不变的情况下，根据矩阵分割和QR瘦型分解理论来获取的一种基于Kriging模型的高效建模方法。主要通过忽略相关矩阵、目标函数矩阵、样本点和估值点的增量矩阵来计算Kriging模型的关键参数的，这一点改进大幅度地减少建模时间。通过一个5杆桁架实例进行了IKM方法与计算机试验设计与分析(DACE)建模时间的比较，很好说明了IKM方法的建模高效性；

步骤4、优化采样准则：

优化采样准则要求利用有效的方法来实现下一个最佳采样点的选择。在优化采样过程中，受到新采样点在整个设计空间均匀分布性和定位采样点需要最大的估计误差因素的影响，利用优化最大化均方误差(MSE)的方法来获取最优采样点，以便所定位的采样点能够在最短的时间内为模型更新带来最大的潜在信息，从而在最少采样点的情况下最大限度地提高模型精度和建模效率；

步骤5、更新准则：

为了有效判断参数θ对Kriging建模的效果和参数θ在后期建模过程中的稳定性，从而决定是否需要对参数θ进行修改或更新，对于具有最大均方误差新采样点的增加对模型中的具有一个很大的影响，因此，在序列优化采样中，利用六西格玛 更新准则的连续概率方法和统计技术来判断是否需要再次对参数θ进行优化更新，实践证明，这一准则取得了良好的效果；

步骤6：DACE建模，如果步骤5中的更新准则满足的话，利用IKM重新建立Kriging模型，否则的话，需要利用DACE重新建立Kriging模型，随着模型精度的逐渐提高，利用IKM方法建模的概率将大大提高，

步骤7：重新回到步骤3，进行下次的循环。

优选方案，其中，增量Kriging方法(IKM)的实现

给定m设计点X＝[x₁，...，x_m]^T，其响应Y＝[y₁，...，y_m]^T，则Kriging模型的回归方程可以表示为：

Y(x)＝Fβ+Z(x) (1)

即一个线性回归部分和一个随机过程部分，其中F为已知的回归模型的基函数(共p个)，β为对应基函数的系数；Z(x)为一随机过程，具有以下性质：

E[Z(x)]＝0

E[Z(x)Z(w)]＝σ²R(θ，ω，x)

其中σ²为该随机过程的方差，R(θ，ω，x)为点x和点ω之间的关联函数，θ为关联参数；

回归函数的基函数有多种选法，常用的有：

(4)常数，即p＝1，f₁(x)＝1

(5)一次函数，p＝n+1，f₁(x)＝1，f₂(x)＝x，...，f_n+1(x)＝x_n

(6)二次函数，p＝(n+1)(n+2)/2

f₁(x)＝1，

f₂(x)＝x₁，...，f_n+1(x)＝x_n

$f_{n+2}\left(x\right)=x_1^2,...,f_{2n+1}\left(x\right)=x_1{x}_n$

$f_{2n+2}\left(x\right)=x_2^2,...,f_{3n}\left(x\right)=x_2{x}_n$

$...,f_p\left(x\right)=x_n^2$

关联函数模型的基本形式为

$R(\mathrm{\θ},\mathrm{\ω},x)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{R}_i({\mathrm{\θ}}_i,{\mathrm{\ω}}_i-x_i)---\left(2\right)$

最常用的高斯关联函数模型可以表示为：高斯模型，R_i(θ，d_i)＝exp(-θ_i，|d_i|²)

关联函数模型以及θ＝[θ₁，θ₂，...，θ_n]；的选取直接影响到Kriging模型的精度。模型的选取需要认为凭经验选择，而θ的选择则可以通过后面的算法实现；

根据以上公式，相关矩阵R和回归函数F可表示为；

根据无偏估计理论，Fβ≈Y有一个最小二乘解和一个过程方差，分别可表示为：

$\widehat{\mathrm{\β}}={\left(F^T{R}^{-1}F\right)}^{-1}{F}^T{R}^{-1}Y---\left(4\right)$

${\mathrm{\σ}}^2=\frac{1}{m}{(Y-F\widehat{\mathrm{\β}})}^T{R}^{-1}(Y-F\widehat{\mathrm{\β}}).---\left(5\right)$

根据公式(2)，矩阵R及相应的和σ2都依赖于θ，基于最大似然估计理论，我们可以对下式子的最大化来获得最优的θ值，

-(mlnσ²+ln|R|)/2 (6)

此外，R是一个对称正定的矩阵，因此R的乔里斯因子分解可由R＝CC^T获取，这里的C即为乔里斯因子，令则公式(4)可以变形为

$\widehat{\mathrm{\β}}={\left({\tilde{F}}^{T}F\right)}^{-1}{\tilde{F}}^T\tilde{Y}.---\left(7\right)$

为了阻止R出现病态矩阵的情况，瘦型QR分解可由

$\tilde{F}=Q{G}^T---\left(8\right)$

根据公式(5)，(7)和(8)，最终能够得到

$G^T\widehat{\mathrm{\β}}=Q^T\tilde{Y}---\left(9\right)$

${\mathrm{\σ}}^2=\frac{1}{m}{\left|\right|\tilde{Y}-\tilde{F}\widehat{\mathrm{\β}}\left|\right|}^2$

优选方案，其中，根据需要增加了k个观测点，(假设对于θ值的影响很小)记之前的各项下标为0，有：

X＝[X₀ ΔX]^T，Y＝[Y₀ ΔY]^T，F＝[F₀ ΔF]^T

$R=\left[\begin{array}{cc}{R}_0& L\\ L^T& \mathrm{\ΔR}\end{array}\right]={\mathrm{CC}}^T=\left[\begin{array}{cc}{C}_1& 0\\ C_2& C_3\end{array}\right]\left[\begin{array}{cc}{C}_1^T& C_2^T\\ 0& C_3^T\end{array}\right]=\left[\begin{array}{cc}{C}_1{C}_1^T& C_1{C}_2^T\\ C_2{C}_1^T& C_2{C}_2^T+C_3{C}_3^T\end{array}\right]---\left(11\right)$

解得：

$\left\{\begin{array}{c}{C}_1=C_0\\ C_2=L^T{C}_0^{-T}\\ C_3=\mathrm{\ΔC}\end{array}\right.$ 中的乔里斯因子

下三角矩阵C的逆矩阵易得，

$C^{-1}=\left[\begin{array}{cc}{C}_0^{-1}& 0\\ N& \mathrm{\Δ}{C}^{-1}\end{array}\right],N=-\mathrm{\Δ}{C}^{-1}{C}_2{C}_0^{-1}$

$\tilde{Y}=C^{-1}Y=\left[\begin{array}{cc}{C}_0^{-1}& 0\\ N& \mathrm{\Δ}{C}^{-1}\end{array}\right]\left[\begin{array}{c}{Y}_0\\ \mathrm{\ΔY}\end{array}\right]=\left[\begin{array}{c}{\tilde{Y}}_0\\ \mathrm{\Δ}\tilde{Y}\end{array}\right],\mathrm{\Δ}\tilde{Y}=N{Y}_0+\mathrm{\Δ}{C}^{-1}\mathrm{\ΔY}$

$\tilde{F}=C^{-1}F=\left[\begin{array}{cc}{C}_0^{-1}& 0\\ N& \mathrm{\Δ}{C}^{-1}\end{array}\right]\left[\begin{array}{c}{F}_0\\ \mathrm{\ΔF}\end{array}\right]=\left[\begin{array}{c}{\tilde{F}}_0\\ \mathrm{\Δ}\tilde{F}\end{array}\right],\mathrm{\Δ}\tilde{F}=N{F}_0+\mathrm{\Δ}{C}^{-1}\mathrm{\ΔF}$

对进行瘦型QR分解：

$\tilde{F}=\left[\begin{array}{c}{\tilde{F}}_0\\ \mathrm{\Δ}\tilde{F}\end{array}\right]=Q{G}^T=\left[\begin{array}{c}{Q}_1\\ Q_2\end{array}\right]G^T---\left(12\right)$

由于瘦型分解具有唯一性，故

${\tilde{F}}_0=Q_1{G}^T,\mathrm{\Δ}\tilde{F}=Q_2{G}^T$

所以

$Q=\left[\begin{array}{c}{Q}_0\\ \mathrm{\ΔQ}\end{array}\right],$ 这里的ΔQ为的QR分解，

求解得到：

$G_0^T\widehat{\mathrm{\β}}=Q^T\tilde{Y}=\left[\begin{array}{cc}{Q}_0^T& \mathrm{\Δ}{Q}^T\end{array}\right]\left[\begin{array}{c}{\tilde{Y}}_0\\ \mathrm{\Δ}\tilde{Y}\end{array}\right]=Q_0^T{\tilde{Y}}_0+\mathrm{\Δ}{Q}^T\mathrm{\Δ}\tilde{Y}=G_0^T{\widehat{\mathrm{\β}}}_0+G_0^T\mathrm{\Δ}\widehat{\mathrm{\β}}$

故得到

$\widehat{\mathrm{\β}}={\widehat{\mathrm{\β}}}_0+\mathrm{\Δ}\widehat{\mathrm{\β}}.$

有益效果，Kriging方法的建模是有效地，随着采样点的增加(对于二维问题，当采样点数目大于600)，Kriging的建模时间陡然增加，因此这里提出一种改进的Kriging方法，在采样点数目非常大的情况下最大限度地降低时间消耗。此外，尽管一次实验设计是非常直观的，但是利用它构造Kriging模型是不合适的。首先，对于大量的仿真估值，一次实验设计是非常耗时的；其次，当采样大量的数据点时，它或许能够导致模型的不正常使用。因此，序列实验设计或许是一种更好的选择。与一次实验设计相比，序列实验设计能控制整个采样过程并充分利用先前的模型信息来决定下面需要采样的有效点，所以，它是一种稳定、有效而又精确的实验设计方法，并在许多工程设计中得到了广泛的应用。

(4).以上述为背景，所提出的基于增量Kriging的序列优化实验设计方法在减少建模时间上将会在工程仿真中得到广泛的应用。SIEDA方法在失去可接受精度范围内已经大幅度提高了Kriging的全局建模效率。

附图说明

图1SIEDA方法的设计框图；

图2SOEDK方法的设计框图；

图3-1为Kriging模型中向量θ中各分量的变化趋势之一；

图3-2为Kriging模型中向量θ中各分量的变化趋势之二；

图4SIEDA和SOEDK的消耗时间比较图；

图5SIEDA、SOEDK和LHS的精度比较图；

图6IKM方法与计算机试验设计与分析(DACE)建模时间的比较。

具体实施方式

为了便于理解本发明，下面结合附图和具体实施例，对本发明进行更详细的说明。本说明书及其附图中给出了本发明的较佳的实施例，但是，本发明可以以许多不同的形式来实现，并不限于本说明书所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。

需要说明的是，当某一元件固定于另一个元件，包括将该元件直接固定于该另一个元件，或者将该元件通过至少一个居中的其它元件固定于该另一个元件。当一个元件连接另一个元件，包括将该元件直接连接到该另一个元件，或者将该元件通过至少一个居中的其它元件连接到该另一个元件。

1.增量Kriging的序列优化实验设计(SIEDA)的详细过程

增量Kriging的序列优化实验设计框图如图1所示，该设计是在传统的SOEDK方法(如 图2)的基础上进行了很大的改进。一般来说，SIEDA可以按照如下七步进行。

步骤1：初始实验设计。对于全局近似问题，利用拉丁方采样(LHS)的实验设计方法(一种空间填充实验设计)来有效获取初始采样点。此外，对采样点所对应的目标函数或仿真模型进行估值。对于自变量个数为n的问题，为了确保空间填充分布的的稳定性和一致性，一般需要选取10n+1个初始采样点来创建初始Kriging模型。

步骤2：Kriging初始建模。利用步骤1中的采样点和相应的响应值并通过DACE方法进行初始Kriging的构造。

步骤3：模型验证。一般来说，方差或者均方根误差是衡量模型精度的主要参数。然而，由于函数的估值需要消耗大量的时间，因此，通过综合分析模型维度、模型精度、运行时间和计算机运算能力等因素，这里使用留一交叉验证的方法进行模型验证。这样，可以大大减少因模型验证带来的大量时间消耗。

步骤4：优化采样。本发明的目的之一就是利用有效的优化采样准则来指导并确定下一个新的采样点。综合考虑消耗时间、精度和采样的准确定位，排除了最大积分均方误差(IMSE)和最大熵值的方法，而采用更有效的最大均方误差方法。

步骤5：更新准则。本发明的更新准则是用来决定DACE和IKM方法中的哪一个方法将被用来进行接下来的Kriging模型更新。实质就是判断参数θ是否发生变化。

步骤6：DACE建模。如果步骤5中的更新准则满足的话，利用IKM重新建立Kriging 模型，否则的话，需要利用DACE重新建立Kriging模型，随着模型精度的逐渐提高，利用IKM方法建模的概率将大大提高。

步骤7：重新回到步骤3，进行下次的循环。

步骤4、5、6是该序列实验设计的核心，对于每一步的具体实现将在下面具体阐述。

2.增量Kriging方法(IKM)的实现

2.1Kriging方法

给定m设计点X＝[x₁，...，x_m]^T，其响应Y＝[y₁，...，y_m]^T，则Kriging模型的回归方程可以表示为：

Y(x)＝Fβ+Z(x) (1)

即一个线性回归部分和一个随机过程部分，其中F为已知的回归模型的基函数(共p个)，β为对应基函数的系数；Z(x)为一随机过程，具有以下性质：

E[Z(x)]＝0

E[Z(x)Z(w)]＝σ²R(θ，ω，x)

其中σ²为该随机过程的方差，R(θ，ω，x)为点x和点ω之间的关联函数，θ为关联参数。

回归函数的基函数有多种选法，常用的有：

(7)常数，即p＝1，f₁(x)＝1

(8)一次函数，p＝n+1，f₁(x)＝1，f₂(x)＝x，...，f_n+1(x)＝x_n

(9)二次函数，p＝(n+1)(n+2)/2

f₁(x)＝1，

f₂(x)＝x₁，...，f_n+1(x)＝x_n

$f_{n+2}\left(x\right)=x_1^2,...,f_{2n+1}\left(x\right)=x_1{x}_n$

$f_{2n+2}\left(x\right)=x_2^2,...,f_{3n}\left(x\right)=x_2{x}_n$

$...,f_p\left(x\right)=x_n^2$

关联函数模型的基本形式为

$R(\mathrm{\θ},\mathrm{\ω},x)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{R}_i({\mathrm{\θ}}_i,{\mathrm{\ω}}_i-x_i)---\left(2\right)$

最常用的高斯关联函数模型可以表示为：高斯模型，R_i(θ，d_i)＝exp(-θ_i，|d_i|²)

关联函数模型以及θ＝[θ₁，θ₂，...，θ_n]；的选取直接影响到Kriging模型的精度。模型的选取需要认为凭经验选择，而θ的选择则可以通过后面的算法实现。

根据以上公式，相关矩阵R和回归函数F可表示为；

根据无偏估计理论，Fβ≈Y有一个最小二乘解和一个过程方差，分别可表示为：

$\widehat{\mathrm{\β}}={\left(F^T{R}^{-1}F\right)}^{-1}{F}^T{R}^{-1}Y---\left(4\right)$

${\mathrm{\σ}}^2=\frac{1}{m}{(Y-F\widehat{\mathrm{\β}})}^T{R}^{-1}(Y-F\widehat{\mathrm{\β}}).---\left(5\right)$

根据公式(2)，矩阵R及相应的和σ²都依赖于θ，基于最大似然估计理论，我们可以对下式子的最大化来获得最优的θ值。

-(mlnσ²+ln|R|)/2. (6)

此外，R是一个对称正定的矩阵，因此R的乔里斯因子分解可由R＝CC^T获取，这里的C即为乔里斯因子，令则公式(4)可以变形为

$\widehat{\mathrm{\β}}={\left({\tilde{F}}^{T}F\right)}^{-1}{\tilde{F}}^T\tilde{Y}.---\left(7\right)$

为了阻止R出现病态矩阵的情况，瘦型QR分解可由

$\tilde{F}=Q{G}^T---\left(8\right)$

根据公式(5)，(7)和(8)，最终能够得到

$G^T\widehat{\mathrm{\β}}=Q^T\tilde{Y}---\left(9\right)$

${\mathrm{\σ}}^2=\frac{1}{m}{\left|\right|\tilde{Y}-\tilde{F}\widehat{\mathrm{\β}}\left|\right|}^2---\left(10\right)$

2.2增量Kriging方法的实现

假设根据需要增加了k个观测点，(假设对于θ值的影响很小)记之前的各项下标为0，有：

X＝[X₀ ΔX]^T，Y＝[Y₀ ΔY]^T，F＝[F₀ ΔF]^T

$R=\left[\begin{array}{cc}{R}_0& L\\ L^T& \mathrm{\ΔR}\end{array}\right]={\mathrm{CC}}^T=\left[\begin{array}{cc}{C}_1& 0\\ C_2& C_3\end{array}\right]\left[\begin{array}{cc}{C}_1^T& C_2^T\\ 0& C_3^T\end{array}\right]=\left[\begin{array}{cc}{C}_1{C}_1^T& C_1{C}_2^T\\ C_2{C}_1^T& C_2{C}_2^T+C_3{C}_3^T\end{array}\right]---\left(11\right)$

解得：

$\left\{\begin{array}{c}{C}_1=C_0\\ C_2=L^T{C}_0^{-T}\\ C_3=\mathrm{\ΔC}\end{array}\right.$ 中的乔里斯因子

下三角矩阵C的逆矩阵易得，

$C^{-1}=\left[\begin{array}{cc}{C}_0^{-1}& 0\\ N& \mathrm{\Δ}{C}^{-1}\end{array}\right],N=-\mathrm{\Δ}{C}^{-1}{C}_2{C}_0^{-1}$

$\tilde{Y}=C^{-1}Y=\left[\begin{array}{cc}{C}_0^{-1}& 0\\ N& \mathrm{\Δ}{C}^{-1}\end{array}\right]\left[\begin{array}{c}{Y}_0\\ \mathrm{\ΔY}\end{array}\right]=\left[\begin{array}{c}{\tilde{Y}}_0\\ \mathrm{\Δ}\tilde{Y}\end{array}\right],\mathrm{\Δ}\tilde{Y}=N{Y}_0+\mathrm{\Δ}{C}^{-1}\mathrm{\ΔY}$

$\tilde{F}=C^{-1}F=\left[\begin{array}{cc}{C}_0^{-1}& 0\\ N& \mathrm{\Δ}{C}^{-1}\end{array}\right]\left[\begin{array}{c}{F}_0\\ \mathrm{\ΔF}\end{array}\right]=\left[\begin{array}{c}{\tilde{F}}_0\\ \mathrm{\Δ}\tilde{F}\end{array}\right],\mathrm{\Δ}\tilde{F}=N{F}_0+\mathrm{\Δ}{C}^{-1}\mathrm{\ΔF}$

对进行瘦型QR分解：

$\tilde{F}=\left[\begin{array}{c}{\tilde{F}}_0\\ \mathrm{\Δ}\tilde{F}\end{array}\right]=Q{G}^T=\left[\begin{array}{c}{Q}_1\\ Q_2\end{array}\right]G^T---\left(12\right)$

由于瘦型分解具有唯一性，故

${\tilde{F}}_0=Q_1{G}^T,\mathrm{\Δ}\tilde{F}=Q_2{G}^T$

所以

$Q=\left[\begin{array}{c}{Q}_0\\ \mathrm{\ΔQ}\end{array}\right],$ 这里的ΔQ为的QR分解，

求解得到：

$G_0^T\widehat{\mathrm{\β}}=Q^T\tilde{Y}=\left[\begin{array}{cc}{Q}_0^T& \mathrm{\Δ}{Q}^T\end{array}\right]\left[\begin{array}{c}{\tilde{Y}}_0\\ \mathrm{\Δ}\tilde{Y}\end{array}\right]=Q_0^T{\tilde{Y}}_0+\mathrm{\Δ}{Q}^T\mathrm{\Δ}\tilde{Y}=G_0^T{\widehat{\mathrm{\β}}}_0+G_0^T\mathrm{\Δ}\widehat{\mathrm{\β}}$

故得到

$\widehat{\mathrm{\β}}={\widehat{\mathrm{\β}}}_0+\mathrm{\Δ}\widehat{\mathrm{\β}}---\left(13\right)$

3.基于DIRECT算法的优化采样准则的实现

在优化采样过程中，新数据点的确定受到以下两个因素的影响：

(1)优化采样需要尽可能使得新数据点均匀地填充到整个设计空间中；

(2)所优化搜索的数据点要具有尽量小的估计值和尽量大的估计误差。

因此，对于获取新采样点的序列优化采样中，有三种方法是比较适合的。一种是用于开发计算机实验设计的最大熵值原理法，一种是应用到确定计算机模型中的最大均方误差方法，而最后一种方法类似于第二种，它是最大化均方误差的积分。

对于最大熵值采样来说，设计者通常需要在一些尽可能远离当前设计点的位置增加采样点，因此，设计决策过程中的模型响应值的相关信息将会被忽视。此外，由于最大熵值采样不隶属于具体的仿真或函数，因此它不具有更大的灵活性。我们希望有一种优化采样方法，它被认为能够从先前估值或Kriging模型中获取或确定具有最大潜在有效信息的新数据点。或许，最大化均方误差和最大化均方误差的积分都是不错的优化采样方法。比较，除了最大均方误差的积分增加了一个权重函数和一个积分过程之外，这两种优化方法是非常相似的。最后，综合考虑应用范围 和时间消耗，我们选择最大化均方误差的优化方法作为最终的优化采样准则。

对于Kriging模型来说，估计的军方误差的定义可用下式来表示：

那么，相应的采样准则能够用如下的表达式表示：

这里，φ(x)表示点x处的均方误差。通过搜索均方误差，能够从当前的Kriging模型中探索到更多潜在和有用的优化建模信息。此外，相关函数r(x)和公式(14)是与新采样点和当前采样点之间的距离相关的，因此，新的采样点将会均匀统一地分布到整个设计空间中。

这个优化过程是一个序列探索的过程，直到找到满足要求的优化采样点。为了使得Kriging模型所对应均方误差的全局优化过程以一个快速的收敛，由Jones发明的一种有效的全局优化方法(DIRECT算法)被应用到公式(15)中来获取最优的采样点。当在五次连续迭代过程中，一种满足约束条件|φ⁽ⁱ⁾-φ^(i-1)|＜0.001情况下，我们将终止DIRECT算法的执行。

上述的优化采样策略有如下两个优点：

(1)每一个新采样点通常都是Kriging继续建模所需要的最优采样点；

(2)随着新采样点的不断增加，DIRECT算法将使优化模型的精度快速下降。

4.更新准则的实现

在建模的初始阶段，具有最大均方误差的新采样点的增加或许对参数有一个较大的影响。但是，随着新采样点的连续不断地增加，参数θ的稳定性将会变得越来越好。当样本总量达到一定的程度后，参数θ的状态往往是趋于稳定的。在这种情况下，即使参数θ的值有一个轻微的改变，这种变化对Kriging模型的精度来几乎是可以被忽略的。图3以二维GP函数为例通过DACE建模方法给出Kriging模型中向量θ中各分量的变化趋势，从中我们可以清楚地看出，向量θ中各分量的变化趋势基本符合我们上述的分析，这从另外一个角度也说明了SIEDA方法的可行性。

依据上述，如何在SIEDA方法中引入一个合适的更新准则来确定是否参数θ在Kriging的连续建模中需要更新将是我们即将解决的问题。六西格玛准则是一种应用连续随机方法和统计技术的一种估值策略，通常用于判断和改进产品的质量。Jone曾经将该准则应用到Kriging模型的留一交叉验证中。在这个验证过程中，如果目标值位于均值加减3倍的标准差区间范围内，我们认为Kriging模型具有99.7％置信度。因此，这个估值原理也是适合于SIEDA算法的更新准则的。

使用更新准则的目的就是使用当前的Kriging模型来有效地指导或判断下一个实际相应值是否位于上述的置信区间内(均值加减3倍的标准差区间)。如果响应值 位于该置信区间内，则保持参数θ的不变，否则的话，利用公式(6)来更新参数θ的值。

已知当前Kriging模型具有k个样本点，一个具有最大均方误差S²(x_k+1)的新采样点(x_k+1，y_k+1)和当前的Kriging模型的估计值(由公式(1)获取)，那么，以及六西格玛准则的更新准则可以表述如下：

${\mathrm{UR}}_{k+1}=\left|\frac{y_{k+1}-{\hat{y}}_k}{S\left(x_{k+1}\right)}\right|>3.---\left(16\right)$

如果不等式(16)得到满足，这意味在UR_k+1不在区间[-3，+3]内，即(x_k+1，y_k+1)对于当前模型是不合适的或无效的，那么，Kriging模型将需要重新构造通过DACE建模。否则的话，Kriging模型仅仅需要通过IKM方法进行模型的重新构造。

值得注意的是，这个更新准则也具有一定的随机性。经过当了的试验测试发现，在Kriging建模的初始阶段，更新准则位于区间[-3，+3]的概率是相当低，但是，仍然存在某些点满足上述的更新准则。同样，当样本总量达到一定数量后，URk+1位于区间[-3，+3]的概率几何是百分之百，然而，也存在少部分采样点不在该区间内。总的来说，随着采样点的不断增加，参数θ对Kriging模型的全局近似影响会很小，那么，由于IKM建模方法的实用，将会节约大量的建模时间。

5.SIEDA方法的验证

尽管我们已经讨论的SIEDA方法，但是它的可行性、建模效率和建模精度需要通过进一步的测试来获得验证。这里选用二维的Golden and Price(GP)函数来验证了序列优化采样中的IKM方法的可行性，如图3所示。接下来，对SIEDA方法和SOEDK方法的时间消耗进行了比较，来说明SIEDA方法的高效性，这里以四个二维函数为例进行的说明，如图4所示，图中明显可以看出SIEDA方法在后期的建模中，利用了很少的时间，仅占SOEDK方法所消耗时间的20％左右。最后，进行了SIEDA、SOEDK和基于LHS序列采样这三种方法的精度进行了比较，这里仍然以四个二维函数为例进行的说明，如图5、图6所示，比较结果发现，基于LHS序列采样的精度最低，而SIEDA和SOEDK方法在初始阶段有一定的差别，然而，这个差异随着采样点的增加而逐渐减小甚至消失。总的来说，SIEDA方法在失去可接受精度范围内已经大幅度提高了Kriging的全局建模效率。

需要说明的是，上述各技术特征继续相互组合，形成未在上面列举的各种实施例，均视为本发明说明书记载的范围；并且，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (3)

1.一种基于增量Kriging的序列优化实验设计方法，其特征在于，其特征在于，该方法包括步骤如下：

步骤1、初始实验设计的选择：

初始实验设计采用拉丁方空间填充采样(LHS)方法，在整个设计空间均匀统一地获取10n+1(n为模型的维度)个初始采样点，并获取相应的仿估值；

步骤2、初始建模以及模型验证：

根据初始采样点和相应的估值，并利用计算机试验设计与分析(DACE)来建立初始的Kriging模型；当获得Kriging近似模型之后，需要进行模型精度的验证，交叉验证法是的验证方法，这里采用了一种高效而又简单的改进留一交叉验证方法实现Kriging模型精度的验证；主要利用当前所有采样点，经过N-1次(N为当前的总样本数)建模和验证之后来决定Kriging模型是否需要继续进行更新；这种验证方法是当前最行之有效且节约时间的模型精度验证方法；

步骤3、增量的Kriging建模(IKM)

IKM方法是在参数θ保持不变的情况下，根据矩阵分割和QR瘦型分解理论来获取的一种基于Kriging模型的高效建模方法；主要通过忽略相关矩阵、目标函数矩 阵、样本点和估值点的增量矩阵来计算Kriging模型的关键参数β的，这一点改进大幅度地减少建模时间；通过一个5杆桁架实例进行了IKM方法与计算机试验设计与分析(DACE)建模时间的比较，很好说明了IKM方法的建模高效性；

步骤4、优化采样准则：

优化采样准则要求利用有效的方法来实现下一个最佳采样点的选择。在优化采样过程中，受到新采样点在整个设计空间均匀分布性和定位采样点需要最大的估计误差因素的影响，利用优化最大化均方误差(MSE)的方法来获取最优采样点，以便所定位的采样点能够在最短的时间内为模型更新带来最大的潜在信息，从而在最少采样点的情况下最大限度地提高模型精度和建模效率；

步骤5、更新准则：

为了有效判断参数θ对Kriging建模的效果和参数θ在后期建模过程中的稳定性，从而决定是否需要对参数θ进行修改或更新；对于具有最大均方误差新采样点的增加对模型中的具有一个很大的影响；因此，在序列优化采样中，利用六西格玛更新准则的连续概率方法和统计技术来判断是否需要再次对参数θ进行优化更新，实践证明，这一准则取得了良好的效果；

步骤6：DACE建模；如果步骤5中的更新准则满足的话，利用IKM重新建立Kriging模型， 否则的话，需要利用DACE重新建立Kriging模型，随着模型精度的逐渐提高，利用IKM方法建模的概率将大大提高；

步骤7：重新回到步骤3，进行下次的循环。

2.根据权利要求1所述基于增量Kriging的序列优化实验设计方法，其特征在于，增量Kriging方法(IKM)的实现，

给定m设计点X＝[x₁，...，x_m]^T，其响应Y＝[y₁，...，y_m]^T，则Kriging模型的回归方程可以表示为：

Y(x)＝Fβ+Z(x) (1)

即一个线性回归部分和一个随机过程部分，其中F为已知的回归模型的基函数(共p个)，β为对应基函数的系数；Z(x)为一随机过程，具有以下性质：

E[Z(x)]＝0

E[Z(x)Z(w)]＝σ²R(θ，ω，x)

其中σ²为该随机过程的方差，R(θ，ω，x)为点x和点ω之间的关联函数，θ为关联参数，

回归函数的基函数有多种选法，常用的有：

(1)常数，即p＝1，f₁(x)＝1

(2)一次函数，p＝n+1，f₁(x)＝1，f₂(x)＝x，...，f_n+1(x)＝x_n

(3)二次函数，p＝(n+1)(n+2)/2

f₁(x)＝1，

f₂(x)＝x₁，...，f_n+1(x)＝x_n

关联函数模型的基本形式为

最常用的高斯关联函数模型可以表示为：高斯模型，R_i(θ，d_i)＝exp(-θ_i，|d_i|²)

关联函数模型以及θ＝[θ₁，θ₂，...，θ_n]；的选取直接影响到Kriging模型的精度，模型的选取需要认为凭经验选择，而θ的选择则可以通过后面的算法实现，

根据以上公式，相关矩阵R和回归函数F可表示为；

根据无偏估计理论，Fβ≈Y有一个最小二乘解和一个过程方差，分别可表示为：

根据公式(2)，矩阵R及相应的和σ2都依赖于θ，基于最大似然估计理论，我们可以对下式子的最大化来获得最优的θ值；

-(m lnσ²+ln|R|)/2 (6)

此外，R是一个对称正定的矩阵，因此R的乔里斯因子分解可由R＝CC^T获取，这里的C即为乔里斯因子，令则公式(4)可以变形为

为了阻止R出现病态矩阵的情况，瘦型QR分解可由

根据公式(5)，(7)和(8)，最终能够得到

3.根据权利要求1所述基于增量Kriging的序列优化实验设计方法，其特征在于，

根据需要增加了k个观测点，(假设对于θ值的影响很小)记之前的各项下标为0，有：

X＝[X₀ ΔX]^T，Y＝[Y₀ ΔY]^T，F＝[F₀ ΔF]^T

解得：

中的乔里斯因子

下三角矩阵C的逆矩阵易得，

对进行瘦型QR分解：

由于瘦型分解具有唯一性，故

所以

这里的ΔQ为的QR分解，

求解得到：

故得到