舰空导弹性能评估, 在其研发改进和定型过程中发挥着重要作用, 是验证其是否满足作战需求的重要手段^[1-3]。在真实作战场景下, 舰空导弹所处环境复杂, 作战目标多样, 为了能够对其性能作出最全面的评估, 理想的试验方法是进行全因素组合的飞行试验, 以获得足够多的性能信息。然而, 由于导弹成本昂贵, 其允许的试验样本数量极为有限。因此, 如何在验证部分必要的试验点基础上, 进行有效的补充试验设计, 在获取关键性能信息的同时尽可能节省试验样本, 已成为舰空导弹研制和试验部门关注的焦点问题。

结合实际作战环境, 舰空导弹在试验过程中, 其性能受多种因素的影响, 这些因素既包括定量因素, 如靶标的距离、高度、速度、遭遇距离、航路捷径等, 也包括定性因素, 如靶标种类、干扰样式等。同时存在的定性-定量因子给舰空导弹补充试验设计增加了不小的难度。当前, 补充试验设计方法主要可分为基于输入的试验设计方法和基于输出的试验设计方法。其中, 基于输入的试验设计方法主要有随机抽样法^[3-4]、低差异序贯法^[5]、嵌套拉丁超立方体设计法^[6-7]、泰森多边形抽样法^[8]、蒙特卡罗抽样法^[9-10]、均方误差(mean square error, MSE)法等方法^[11]; 基于输出的试验设计方法主要为贝叶斯优化法^[12]、黄金分割法^[13]、局部线性近似-泰森(local linear approximation-Voronoi, LOLA-Voronoi)方法^[14]等, 其中贝叶斯优化法主要用于黑箱函数极值搜索, 而黄金分割法和LOLA-Voronoi方法无法利用代理模型在样本信息利用上的优势, 补充试验设计效果不稳定。

上述基于输入的试验方法仅使用当前试验样本的位置来确定下一步的采样位置, 不考虑输出响应值的影响, 一味追求样本分布的均匀性, 对试验空间重点区域的开发考虑不足; 而基于输出的试验设计方法同时兼顾了输入和响应对采样的共同影响, 通过构建近似梯度来平衡采样中的探索与开发(LOLA-Voronoi方法), 不依赖代理模型。但是, 对于某些具有有效先验知识的采样空间而言, 构建代理模型可提高搜索效率, 相比于通过近似梯度进行补充试验的设计方法, 基于模型的补充试验设计方法可有效改善试验结果的适用性^[15-16]。此外, 目前对基于模型的补充试验设计方法的研究的主要关注点仍在定量因子设计, 对混合定性-定量因子的试验设计的研究较少。与仅包含定量因子的预测评估模型相比, 定性-定量因子混合的预测评估模型有其自身特点, 探索定性-定量因子混合预测评估模型基础上的补充试验设计方法具有一定价值。

事实上, 类似舰空导弹的复杂装备在试验方案设计中面临的问题普遍存在, 本文针对当前补充试验设计方法在处理包含定性-定量因子方面存在的不足, 研究了一种混合定性-定量因子模型的补充试验设计方法。该方法首先构建了定性-定量因子混合评估模型。然后, 利用探索-开发平衡策略进行补充试验设计建模, 并在此基础上确定目标试验点与最佳邻域, 再借助遗传算法搜索出最佳邻域位置。最后, 通过示例对所提出的模型和方法进行了验证。

假设需要进行的补充试验设计中, 试验输入既包含定性因子, 又包含定量因子, 可以表示为w=(x^t, z^t)^t, 其中x= (x₁, x₂, …, x_m)^t表示m个定量因子, z=(z₁, z₂, …, z_J)^t表示J个定性因子。

假设已获得n组试验数据, 在此基础上构建高斯过程模型:

(1)$y(\boldsymbol{w})=\boldsymbol{\beta}^t \boldsymbol{f}(\boldsymbol{w})+\boldsymbol{\varepsilon}(\boldsymbol{w})$

式中: 未知系数向量为β, β =(β₁, β₂, …, β_l)^t, 相关函数向量为f(w), f(w)=(f₁(w), f₂(w), …, f_l(w))^t, 两者相乘组成高斯过程模型的回归项; ε(w)~GP(0, σ²)是满足均值为0、方差为σ²的高斯过程残差项。

由上可知, 残差项ε(w)~GP(0, σ²)中相关函数的构建是关键。与定量因子相关函数的构建方法不同, 定性因子中并没有“距离”的概念, 因此, 基于距离的相关函数方法对定性因子并不适用^[17]。为此, 需要建立既包含定性因子, 又包含定量因子的相关函数模型。

为了便于理解相关函数的构建过程, 首先构建仅包含一个定性因子的相关函数模型。设定性因子z₁具有I₁等级类别水平, 进一步, 令ε_u(w)=ε((x^t, u^t)^t), u=1, 2, …, I₁, 可进一步得到:

(2)$\boldsymbol{\varepsilon}^*(\boldsymbol{x})=\left(\varepsilon_1(\boldsymbol{x}), \varepsilon_2(\boldsymbol{x}), \cdots, \varepsilon_{I_1}(\boldsymbol{x})\right)^t$

由式(2)可知, 仅需要定义ε^*(x)的相关函数和互相关函数即可。为了便于后续计算, 假设ε^*(x)=Aη(x)。其中, A=(a₁, a₂, …, a_I1)是一个具有单位行向量的I₁×I₁的非奇异矩阵(即当u=1, 2, …, I₁时, 向量乘积a_u^ta_u=1);η(x)=(η₁(x), η₂(x), …, η_I1(x))^t。其中, η₁(x), η₂(x), …, η_I1(x)为具有相同方差σ²和相关函数K_ϕ的独立随机过程。由此可知, corr (η(x_i), η(x_j))=K_ϕ(x_i, x_j)I_I1。其中, I_I1是一个I₁×I₁的单位矩阵。对于任意两个输入w_i和w_j, 其相关函数可定义为

(3)$\begin{gathered}\operatorname{corr}\left(\varepsilon\left(\boldsymbol{w}_i\right), \varepsilon\left(\boldsymbol{w}_j\right)\right)=\operatorname{corr}\left(\boldsymbol{\varepsilon}_{z_{1 i}}\left(\boldsymbol{x}_i\right), \boldsymbol{\varepsilon}_{z_{1 j}}\left(\boldsymbol{x}_j\right)\right)= \\\operatorname{corr}\left(\boldsymbol{a}_{z_{1 i}}^t \boldsymbol{\eta}\left(\boldsymbol{x}_i\right), \boldsymbol{a}_{z_{1 j}}^t \boldsymbol{\eta}\left(\boldsymbol{x}_j\right)\right)= \\\boldsymbol{a}_{z_{1 i}}^t \boldsymbol{a}_{z_{1 j}} K_\phi\left(\boldsymbol{x}_i, \boldsymbol{x}_j\right)\end{gathered}$

令τ_{r, s}=a_r^ta_s, 其中, r, s=1, 2, …, I₁, 可进一步得到定性因子相关矩阵T₁=(τ_{r, s})=AA^t, 这是一个具有单位对角元素的正定矩阵。由上可知, 对于任意单位对角元素的正定矩阵T₁=(τ_{r, s})和任意相关函数K_ϕ(x_i, x_j), 式(3)中的相关函数可进一步写为corr(ε(w_i), ε(w_j))=τ_{z1i, z1j}K_ϕ(x_i, x_j)。其中, τ_{z1j, z1i}=τ_{z1i, z1j}描述等级类别z_1j和z_1i的互相关性。定量因素响应模型构建相同, K_ϕ(x_i, x_j)的构建仍然选用高斯核函数。因此, ε(w_i)和ε(w_j)的相关函数可进一步表示为

(4)$\operatorname{corr}\left(\varepsilon\left(\boldsymbol{w}_i\right), \varepsilon\left(\boldsymbol{w}_j\right)\right)=\tau_{z_{1 i}, z_{1 j}} \exp \left\{-\sum\limits_{k=1}^m \phi_k\left(x_{i k}-x_{j k}\right)^2\right\}$

由式(4)可知, τ_{z1j, z1i}互相关模型的构建和求解是定量-定性因子混合响应模型构建的关键。目前, 常用的解决方法包括可交换相关函数法、乘法相关函数法、组相关函数法、序数定性因子相关函数法和超球面分解法等。相比而言, 超球面分解法常用来描述复杂估计过程模型中的相关性, 可以用投影的方式直观表现定性因子的位置关系^[18], 可有效克服其他方法未知参数多、推导困难、计算量巨大的缺点, 便于理解和计算。因此, 本文选用超球面分解法进行定性因子相关性建模, 主要分以下两步进行。

步骤1  已知T₁=(τ_{r, s})是一个具有单位对角元素的正定矩阵, 对其进行Cholesky分解, 可得

(5)$\boldsymbol{T}_1=\boldsymbol{L} \boldsymbol{L}^t$

式中: L={l_{r, s}}为具有严格正对角项的下三角矩阵。

步骤2  对下三角矩阵L={l_{r, s}}的每个行向量[l_{r, 1}, l_{r, 2}, …, l_{r, r}]进行描述, 将其定义为一个多维单位超球面上的坐标, 如下所示:

(6)$\left\{\begin{array}{l}l_{r, 1}=\cos \theta_{r, 1}, s=1 \\l_{r, s}=\sin \theta_{r, 1} \sin \theta_{r, 2} \cdots \sin \theta_{r, s-1} \cos \theta_{r, s}, 2 \leqslant s \leqslant r-1 \\l_{r, r}=\sin \theta_{r, 1} \sin \theta_{r, 2} \cdots \sin \theta_{r, r-2} \sin \theta_{r, r-1}, s=r\end{array}\right.$

式中: θ_{r, s}∈(0, π)。由此可知, 对角矩阵L={l_{r, s}}中对角线中元素l_{r, r}均为正值, 可进一步推导, L={l_{r, s}}为正定矩阵。此外, 式(6)中$ \tau_{r, r}=\sum_{s=1}^r l_{r, s}^2=1$, 进一步可推导出, T₁=(τ_{r, s})具有单位对角线元素, 满足单位对角元素正定矩阵的条件。下面详细给出等级水平J=1, 2, 3时对应的定性因子相关矩阵T₁。

当J=1时, T₁=1, 即与仅包含定量因素的相关函数构建方法相同。

当J=2时, $ \boldsymbol{T}_1=\left[\begin{array}{cc} 1 & \tau_{1, 2} \\ \tau_{2, 1} & 1 \end{array}\right]$, 其中τ_{1, 2}=τ_{2, 1}; 对其进行Cholesky分解, 可得

$\boldsymbol{T}_1=\boldsymbol{L} \boldsymbol{L}^t=\left[\begin{array}{cc}1 & 0 \\l_{2, 1} & l_{2, 2}\end{array}\right]\left[\begin{array}{ll}1 & l_{1, 2} \\0 & l_{2, 2}\end{array}\right], l_{1, 2}=l_{2, 1}$

进一步构建多维单位球面坐标, 可得

$\left\{\begin{array}{l}l_{2, 1}=\cos \theta_{2, 1} \\l_{2, 2}=\sin \theta_{2, 1}\end{array}\right.$

由l_{2, 1}≥0可知, θ_{2, 1}∈(0, π), 计算方法如下所示:

$\tau_{2, 1}=\tau_{1, 2}=\cos \left(\theta_{2, 1}\right)$

当J=3时, $ \boldsymbol{T}_1=\left[\begin{array}{ccc} 1 & \tau_{1, 2} & \tau_{1, 3} \\ \tau_{2, 1} & 1 & \tau_{2, 3} \\ \tau_{3, 1} & \tau_{3, 2} & 1 \end{array}\right]$, 其中, τ_{1, 2}=τ_{2, 1}, τ_{1, 3}=τ_{3, 1}, τ_{3, 2}=τ_{2, 3}; 对其进行Cholesky分解, 可得

$\boldsymbol{T}_1=\boldsymbol{L}^t=\left[\begin{array}{ccc}1 & 0 & 0 \\l_{2, 1} & l_{2, 2} & 0 \\l_{3, 1} & l_{3, 2} & l_{3, 3}\end{array}\right]\left[\begin{array}{ccc}1 & l_{1, 2} & l_{1, 3} \\0 & l_{2, 2} & l_{2, 3} \\0 & 0 & l_{3, 3}\end{array}\right]$

其中, l_{1, 2}=l_{2, 1}, l_{1, 3}=l_{3, 1}, l_{3, 2}=l_{2, 3}; 进一步构建多维单位球面坐标可得

$\begin{gathered}\left\{\begin{array}{l}l_{2, 1}=\cos \theta_{2, 1} \\l_{2, 2}=\sin \theta_{2, 1}\end{array}\right. \\\left\{\begin{array}{l}l_{3, 1}=\cos \theta_{3, 1} \\l_{3, 2}=\sin \theta_{3, 1} \cos \theta_{3, 2} \\l_{3, 3}=\sin \theta_{3, 1} \sin \theta_{3, 2}\end{array}\right.\end{gathered}$

式中: θ_{i, j}∈(0, π)。计算方法如下所示:

$\left\{\begin{array}{l}\tau_{1, 2}=\cos \theta_{2, 1} \\\tau_{1, 3}=\cos \theta_{3, 1} \\\tau_{2, 3}=\cos \theta_{2, 1} \cos \theta_{3, 1}+\sin \theta_{2, 1} \sin \theta_{3, 1} \cos \theta_{3, 2}\end{array}\right.$

当定性因素数目u＞1时, 构造任意两个输入的相关函数模型:

(7)$\begin{gathered}\operatorname{corr}\left(\varepsilon\left(\boldsymbol{w}_i\right), \varepsilon\left(\boldsymbol{w}_j\right)\right)=\operatorname{corr}\left(\varepsilon\left(\boldsymbol{x}_i, \boldsymbol{z}_i\right), \varepsilon\left(\boldsymbol{x}_j, \boldsymbol{z}_j\right)\right)= \\{\left[\prod\limits_{p=1}^J \tau_{p, z_{p, i}, z_{p, j}}\right] \exp \left\{-\sum\limits_{k=1}^m \phi_k\left(x_{i, k}-x_{j, k}\right)^2\right\}}\end{gathered}$

与定量因素预测响应模型构建相似^[18], 定量-定性因素混合预测响应模型中包含的未知相关系数为ϕ和θ。

任意未知输入w₀的预测响应模型为

(8)$\hat{y}\left(\boldsymbol{w}_0\right)=\boldsymbol{f}^{\mathrm{T}}\left(\boldsymbol{w}_0\right) \hat{\boldsymbol{\beta}}+\boldsymbol{r}^t\left(\boldsymbol{w}_0\right) \boldsymbol{R}^{-1}(\boldsymbol{Y}-\boldsymbol{F} \hat{\boldsymbol{\beta}})$

式中: r^t(w₀)=(corr(ε(w₀), ε(w₁)), corr(ε(w₀), ε(w₂)), …, corr(ε(w₀), ε(w_n))); F=(f(w₁), f(w₂), …, f(w_n)); R是包含所有corr(ε(w_r), ε(w_s))(r, s=n)的相关函数矩阵。

与定量因素预测响应模型中未知参数估计值的求取方法相似^[19], 利用极大似然估计可进一步求得

(9)$\left\{\begin{array}{l}\hat{\boldsymbol{\beta}}=\left(\boldsymbol{F}^t \boldsymbol{R}^{-1} \boldsymbol{F}\right)^{-1} \boldsymbol{F}^{\mathrm{T}} \boldsymbol{R}^{-1} \boldsymbol{Y} \\\hat{\sigma}^2=\frac{1}{n}(\boldsymbol{Y}-\boldsymbol{F} \hat{\boldsymbol{\beta}})^{-\mathrm{T}} \boldsymbol{R}^{-1}(\boldsymbol{Y}-\boldsymbol{F} \hat{\boldsymbol{\beta}})\end{array}\right.$

(10)$(\hat{\boldsymbol{\phi}}, \hat{\boldsymbol{\theta}})=\underset{\phi, \theta}{\operatorname{argmin}}\left\{\frac{-\left(n \ln \hat{\sigma}^2+\ln |\boldsymbol{R}|\right)}{2}\right\}$

假设当前已获得的试验响应值构成的样本集合为N={W=(w₁, w₂, …, w_n)^t, Y=(y₁, y₂, …, y_n)^t}, 在该样本集合中依次选择一组试验样本(w_i, y_i)(其中i=1, 2, …, n)作为验证样本, 利用交叉验证方法对模型预测精度进行评估验证, 预测值均方根误差(root mean square error, RMSE)为

(11)$\operatorname{RMSE}(n)=\sqrt{\frac{\sum\limits_{i=1}^n\left\{\left[\hat{y}\left(\boldsymbol{w}_i\right)-y\left(\boldsymbol{w}_i\right)\right] / y\left(\boldsymbol{w}_i\right)\right\}^2}{n}}$

式中: 当RMSE(n)≤δ时, 初始样本数量满足对预测模型精度的要求, 不需要再进行补充试验设计; 当RMSE(n)＞δ时, 初始样本数量无法满足对预测模型精度的要求, 需要进行后续补充试验设计。其中, δ为预测精度误差允许值。

当预测响应模型的预测精度无法满足要求时, 需要增加试验样本的数量来提高预测模型的准确性, 在满足预测精度的前提下尽可能少地补充试验样本点, 是补充试验模型构建的关键。由预测模型构建的基本特征可知, 在进行补充采样时, 应对两个因素进行平衡, 这两个因素分别为探索和开发。探索是对全局未开发区域不确定性的度量; 开发是对感兴趣区域的度量。为此, 本文在已构建预测响应模型的基础上设计补充试验。

由上文可知, 预测点的不确定性由预测均方差表示, 预测均方差的大小表征着预测结果的不确定性。在进行补充试验时，应探索选择预测均方差较大的点, 以降低预测值的不确定性, 提高预测模型的稳定性。

对于同时包含定性-定量因子的试验, 对任意未测点w, 该未测点的预测均方差(mean square error of prediction, MSEP)可表示为

(12)$\mathrm{MSEP}=s^2(\boldsymbol{w})=\mathrm{E}\left[(y(\boldsymbol{w})-\hat{y}(\boldsymbol{w}))^2\right]$

利用多元正态条件分布可知:

$\left( {\begin{array}{*{20}{l}}{{{\mathit{\boldsymbol{W}}}_1}}\\{{{\mathit{\boldsymbol{W}}}_2}}\end{array}} \right)\sim {{\rm{N}}_{m + n}}\left( {\left( {\begin{array}{*{20}{l}}{{{\mathit{\boldsymbol{\mu }}}_1}}\\{{{\mathit{\boldsymbol{\mu }}}_2}}\end{array}} \right), \left[ {\begin{array}{*{20}{l}}{{{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{1, 1}}}&{{{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{1, 2}}}\\{{{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{2, 1}}}&{{{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{2, 2}}}\end{array}} \right]} \right)$

式中: μ₁是一个m×1向量; μ₂是一个n×1向量; Σ_{1, 1}是一个m×m矩阵; Σ_{1, 2}=Σ_{2, 1}^t是一个m×n矩阵; Σ_{2, 2}是一个n×n矩阵，则可以得到条件分布[W₁|W₂]为

(13)$\begin{array}{l}\left[ {{{\mathit{\boldsymbol{W}}}_1}\mid {{\mathit{\boldsymbol{W}}}_2}} \right]\sim {{\rm{N}}_m}\left( {{{\mathit{\boldsymbol{\mu }}}_1} + {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{1, 2}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{2, 2}^{ - 1}\left( {{{\mathit{\boldsymbol{W}}}_2} - {{\mathit{\boldsymbol{\mu }}}_2}} \right), } \right.\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{1, 1}} - {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{2, 2}^{ - 1}{{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}}_{2, 1}}} \right)\end{array}$

由于[y(w), Y]均满足多元正态分布的条件, 利用多元正态条件分布公式可进一步推导出:

$[y(\boldsymbol{w}) \mid \boldsymbol{Y}] \sim \mathrm{N}\left(\mu_2, \sigma_2^2\right)$

式中: MSEP为

(14)$\begin{gathered}s^2(\boldsymbol{w})=\sigma_2^2=\sigma^2\left(1-\boldsymbol{r}_\boldsymbol{w}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}+\right. \\\left.\left(\boldsymbol{f}_\boldsymbol{w}-\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}\right)^t\left(\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{F}\right)^{-1}\left(\boldsymbol{f}_\boldsymbol{w}-\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}\right)\right)\end{gathered}$

式中: F=(f(w₁), f(w₂), …, f(w_n)); R_w是包含所有corr(ε(w_r), ε(w_s))(r, s=n)的相关矩阵; r_w=(corr(ε(w), ε(w₁)), corr(ε(w), ε(w₂)), …, corr(ε(w), ε(w_n)))^t; f_w=f(w)。

通过构建MSEP模型可以对任意未测点进行MSEP估计, 并且可以通过优化算法找出MSEP最大的点, 并在该位置进行补充试验。通过观察MSEP可知, 未测点MSEP的取值只与输入有关, 与响应值无关。因此, 在进行补充试验时, 每组试验均可以补充多个样本点。

现实中存在着某些响应曲面局部难以近似的问题, 在补充试验中需要对该区域进行更加深入的开发, 以进一步提高预测模型精度。与文献[12]中的局部线性近似(local linear approximation, LOLA)方法不同, 本节考虑预测模型评估值对局部区域的描述, 在此基础上研究探索-开发策略。基于构建的评估预测模型, 探讨补充试验设计新方式。实现过程分为以下3步: 在当前样本空间中确定出目标试验点; 构建目标试验点的最佳邻域; 通过极大化MSEP在邻域中搜索补充样本点。

步骤1  确定包含重点区域的目标试验点。已知当前已获得试验样本空间为N, N={W=(w₁, w₂, …, w_n)^t, Y=(y₁, y₂, …, y_n)^t}, 依次将试验点{w_i, y_i}作为候选目标点, 利用样本空间中的剩余试验样本构建预测响应模型, 并对该候选目标试验点w_i的预测值和MSEP进行求解。

$\begin{gathered}\hat{y}\left(\boldsymbol{w}_i\right)=\boldsymbol{f}^{\mathrm{T}}\left(\boldsymbol{w}_i\right) \hat{\boldsymbol{\beta}}+\boldsymbol{r}^t\left(\boldsymbol{w}_i\right) \boldsymbol{R}^{-1}(\boldsymbol{Y}-\boldsymbol{F} \hat{\boldsymbol{\beta}}) \\s_{n_1-1}^2\left(\boldsymbol{w}_i\right)=\sigma^2\left(1-\boldsymbol{r}_\boldsymbol{w}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}+\right.\\\left.\left(\boldsymbol{f}_\boldsymbol{w}-\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}\right)^t\left(\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{F}\right)^{-1}\left(\boldsymbol{f}_\boldsymbol{w}-\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}\right)\right)\end{gathered}$

式中: i=1, 2，…, n。

在此基础上, 给出目标试验点w_i的预测误差(prediction error, PE), 分别为

$\operatorname{PE}\left(\boldsymbol{w}_i\right)=\left(y_i-\hat{y}\left(\boldsymbol{w}_i\right)\right)^2$

由上可知, 目标试验点是由探索项和开发项共同确定的, 因此, 计算出每个候选目标试验点的MSEP和PE, 并将MSEP与PE相加的和作为候选目标试验点的得分, 将得分最高的点作为搜索的目标试验点, 如下所示:

(15)$\boldsymbol{w}_{\mathrm{obj}}=\arg \max \limits_{\boldsymbol{w}}\left(s_{n_1-1}^2(\boldsymbol{w})+\operatorname{PE}(\boldsymbol{w})\right)$

步骤2  构建目标试验点的最佳邻域。由式(15)确定目标样本位置后, 接下来就需要构建最优目标邻域。

为不失一般性, 将最佳邻域中试验点的数量确定为d=2m, 其中m为定量因素的数目。依照定性因子对当前试验样本空间进行划分, 可得目标样本点邻域中试验点的数量为k, 则候选邻域的数目为C_k^d。与LOLA方法中构建邻域的方式相似^[14], 设定N(w_{obj, i})=(w_{i, 1}, w_{i, 2}, …, w_{i, m})表示目标试验点w_obj的第i个邻域, 在此基础上获得多定量和单定量条件下的交叉多面体比率。

在多定量条件下，有

$R\left(N\left(\boldsymbol{p}_r\right)\right)=\frac{A\left(N\left(\boldsymbol{p}_r\right)\right)}{\sqrt{2} C\left(N\left(\boldsymbol{p}_r\right)\right)}$

式中: p_r表示候选目标试验点w_obj中的定量N(因子x_obj; p_r)表示候选邻域; C(N(p_r))表示候选邻域的内聚值; A(N(p_r))为候选邻域的黏合值。

在单定量条件下，有

$R\left(N\left(\boldsymbol{p}_r\right)\right)=1-\frac{\left|p_{r 1}+p_{r 2}\right|}{\left|p_{r 1}\right|+\left|p_{r 2}\right|+\left|p_{r 1}-p_{r 2}\right|}$

接下来给出候选邻域的分数:

$S\left(N\left(\boldsymbol{p}_r\right)\right)=\frac{R\left(N\left(\boldsymbol{p}_r\right)\right)}{C\left(N\left(\boldsymbol{p}_r\right)\right)}$

对S(N(p_r))按照分数值大小进行排序, 选出最优邻域。接下来，在邻域范围内搜索出新样本的位置, 完成补充试验设计。

步骤3  利用构建的MSEP模型在邻域中搜索下一个补充样本点。

已知邻域N(p_r), 就可以确定邻域取值范围V_obj, 并在此基础上利用上文构建的MSEP模型, 搜索最优值作为补充试验输入。

由式(14)可知:

$\begin{aligned}s_n^2(\boldsymbol{w})= & \sigma^2\left(1-\boldsymbol{r}_\boldsymbol{w}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}+\left(\boldsymbol{f}_\boldsymbol{w}-\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}\right)^t \cdot\right. \\& \left.\left(\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{F}\right)^{-1}\left(\boldsymbol{f}_\boldsymbol{w}-\boldsymbol{F}^t \boldsymbol{R}_\boldsymbol{w}^{-1} \boldsymbol{r}_\boldsymbol{w}\right)\right)\end{aligned}$

在邻域取值范围V_obj内, 利用遗传算法搜索出MSEP的最大值:

(16)$\boldsymbol{w}_{n+1}=\arg \max \limits_{w_{\mathrm{obj}}} s_n^2(\boldsymbol{w})$

将式(16)作为第一组补充试验输入值, 然后将该值分别代入式(8)计算补充试验输入条件下的预测响应值, 并将该试验点补充到当前试验样本空间, 继续迭代上述过程, 直至达到该批试验次数要求。在此基础上, 将上述补充试验设计点代入响应模型, 并获得对应的响应值。

将所有的数值模拟试验数据补充到试验样本空间, 重新构建评估模型, 然后利用式(11)中的交叉验证方法对模型预测精度进行评估验证, 判定评估精度是否满足要求, 以确定是否需要继续进行补充试验。

为了验证本文所提方法的有效性, 以及更加形象地展示补充试验设计方法的采样特征, 将该方法应用于如下数值算例(无量纲)研究, 并与均匀试验设计^[20]和MSE方法^[11]的结果进行比较。

算例1

$f=-\sin x-\exp (x / 100)+10, x \in[0, 10]$

数值算例1中, 测试函数仅包含定量因子。首先, 利用均匀试验设计给出初始样本, 如图 1所示。然后, 分别利用MSE方法和本文所提方法对该算例进行试验设计, 并在此基础上进行评估建模, 试验结果如图 2~图 5所示。此外, 图 6表示补充4个样本的条件下, 两种方法的预测误差值。

为了能够对本文所提方法与均匀试验设计方法进行比较, 首先, 利用均匀试验设计方法进行试验样本采集(采集样本数量为3个), 并代入测试函数求取响应值, 在此基础上构建的预测响应模型如图 1所示。然后, 将图 1中通过均匀试验设计的样本作为初始样本点, 利用MSE方法进行补充试验设计(分别采集2个样本和4个样本), 并在此基础上构建预测响应模型, 如图 2和图 3所示。利用本文所提方法进行补充试验设计(分别采集2个样本和4个样本), 并在此基础上构建预测响应模型, 如图 4和图 5所示。图 6的分析结果表明, 本文所提方法的预测精度更高, 模型更稳定。

算例2

$f=\left\{\begin{array}{l}\exp (1.4 x) \cos (7 \pi x / 2), z=1 \\\exp (3 x) \cos (7 \pi x / 2), z=2\end{array}, x \in[0, 1]\right.$

在数值算例2中, 测试函数包含一个定量因素和一个定性因素, 其中定性因素包含两个水平。首先, 利用均匀试验设计获得10个初始样本, 如图 7所示。

然后, 利用本文提出的补充试验设计方法再补充10个样本, 如图 8所示。最后, 与均匀试验设计方法进行对比, 如图 9所示。

图 9为预测误差箱型图。分析结果表明, 本文所提方法对于解决定性-定量因子混合补充实验设计问题, 对精度的改善效果更好。

通过上述两个算例可以直观地看出, 相比于MSE方法和均匀试验设计方法, 利用本文所提方法构建的预测响应模型与算例函数更为接近。由图 5和图 8可知, 本文所提方法在预测模型急剧变化的位置, 可采集更多的样本点。通过图 6和图 9的对比分析可知, 所提方法对改善预测模型精度的贡献更大。

舰空导弹是舰艇主要的防空武器, 其作战性能受到目标、背景、干扰等复杂战场环境的影响, 涉及多个定性和定量因子。为了更好地验证本文所提方法的有效性, 下面对某型舰空导弹的试验影响因素进行抽象和简化, 在假设其他影响因素固定的前提下, 仅考虑目标因素的影响。选择脱靶量为响应值, 靶标类型、靶标速度和靶标高度为影响因子, 其取值范围如表 1所示。

由上可知, 靶标种类为定性因子, 靶标高度和靶标速度为定量因子。假设初始试验样本共36组, 分别对每个靶标进行了12组初始试验, 然后在此基础上再进行10组补充试验, 以提高评估模型精度。根据本文方法给出的补充试验方案如表 2所示, 在此基础上构建的评估响应模型如图 10~图 12所示。

由图 10~图 12可知, 补充试验的样本数目由大到小排序依次为亚音速导弹、常规飞机和直升机, 这主要与脱靶量的变化范围有关, 即变化范围越大, 补充试验样本的数目就应该越多, 补充试验方案设计结果符合这一基本假设。此外, 补充试验样本的位置多处于边界或变化较大的点周围, 这也与补充试验探索-开发的基本思想一致。

进一步, 利用交叉验证方法对补充试验后的模型预测精度进行验证, 模型预测误差如图 13所示。

由图 13可知, 在补充试验过程中, 预测误差随样本数目的增加表现为非单调性递减的特征, 这与补充试验的探索特性有关。当补充试验点的位置处于某些关键未知区域时, 交叉验证误差会出现增加的情况, 试验结果与补充试验设计方法期望的探索方式一致。

本文结合舰空导弹等复杂装备在补充试验方案设计中存在的问题, 基于定性-定量因子混合模型, 提出了改善预测模型精度的补充试验设计方法。该方法主要平衡了对补充试验空间的探索与开发, 增强了捕捉重要特征样本位置的能力, 改善了预测模型的构建精度。通过算例分析, 验证了本文提出的补充试验设计方法的有效性。最后, 结合某型舰空导弹的补充试验方案设计对方法进行了深入的验证和分析。事实上, 影响舰空导弹等复杂装备的试验因素包含定性因素和定量因素是常态, 并且影响因素数量很多, 甚至因素间存在复杂的约束关系, 这就给模型的复杂度以及算法的求解效率提出了更大的挑战。本文仅是针对此类问题进行了初步的探索, 以最大限度地提高对试验对象性能描述的准确性, 后续将结合实际工程问题进行更加深入的研究。