CN102789539B - 一种飞行器升力面结构优化设计方法 - Google Patents

Abstract

一种飞行器升力面结构优化设计方法：①采用LHD在设计空间内采取试验设计样本并获得飞行器升力面高精度分析模型响应值；②构造径向基函数代理模型；③在设计空间使用LHD方法随机生成一定数量的网格点，计算其对应的径向基函数代理模型响应值；④给定全局空间缩减率M₁%和局部空间缩减率M₂%，进行2次伪缩减；⑤在每个聚类空间内利用遗传算法对径向基函数代理模型寻优；⑥得到最优点与上次空间全局近似最优点的空间距离，若满足缩减准则，则删除当前空间得到新空间；⑦获取当前全局近似最优解的真实响应值，判断其是否收敛，如收敛则停止优化，否则反复迭代直至找到最优解。本发明提出的方法提高了优化效率，节约了飞行器升力面结构优化设计成本。

Description

一种飞行器升力面结构优化设计方法

技术领域

本发明涉及一种飞行器升力面结构优化设计方法，属于飞行器优化设计中多学科耦合优化设计技术领域。

背景技术

上世纪80年代起，美国NASA Langley研究中心的Sobieski等一批航空领域的科学家和工程设计人员，提出一种新的飞行器设计方法：多学科设计优化（Multidisciplinary Design Optimization，MDO）。MDO的提出受到了各国研究机构、高校和工业界的普遍重视（如NASA、斯坦福大学、麻省理工学院、Airbus公司和Boeing公司等），掀起了一股MDO的研究热潮。一方面，为了提高分析精度和可信度，高精度分析模型被广泛用于飞行器设计中，例如结构分析中采用的有限元分析（FEA）模型、气动分析中使用的计算流体力学（CFD）分析模型等。高精度分析模型在提高分析精度和可信度的同时也带来了计算耗时的困难，虽然当今计算机软硬件技术已经有了长足的发展，然而，调用高精度分析模型完成一次分析仍然极其耗时，例如使用CFD模型完成一次气动仿真分析需要数小时甚至数十小时。另一方面，飞行器设计问题往往涉及多个相互耦合的学科。譬如，飞行器设计涉及气动、结构、动力、隐身、控制等学科，各学科相互影响，相互制约，飞行器的性能是各学科耦合的综合体现。由于学科之间的耦合关系，飞行器设计问题的系统分析表现为多学科分析。本质上，多学科分析过程是一个典型的非线性求解过程，每次多学科分析都需要进行多次迭代，计算耗时，如果各学科都采用高精度分析模型，其计算量将非常庞大。再次，飞行器优化设计过程中需要经过反复迭代方能收敛到局部或全局最优解，而每次迭代都需要进行多次飞行器设计问题的多学科分析，可见计算成本将进一步增加。因此，直接采用传统的全局优化方法与高精度飞行器分析模型来解决飞行器优化设计会进一步增加设计成本。如何解决飞行器多学科优化设计这个难题困扰许多设计专家。为了解决这样的问题，基于代理模型的飞行器优化设计方法引起了越来越多的注意。

基于代理模型的飞行器升力面多学科耦合设计优化方法本质就是构造分析结果与高精度模型接近（包括气动、结构、隐身等学科分析模型以及全系统分析模型），计算成本更低的代理模型并代替高精度分析模型用于优化。由于高精度分析模型计算一次所需时间的量级为小时或天，而代理模型计算一次所用时间的量级仅为秒甚至毫秒，因此与高精度分析模型的计算时间相比，代理模型以及基于代理模型优化的计算时间往往可以忽略不计。目前常用的代理模型方法包括多项式响应面、Kriging模型、径向基函数、移动最小二乘法以及神经网络等。

多项式响应面能够快速有效地解决低阶相对简单分析模型；Kriging模型能够对一个未知的高阶非线性分析模型进行高精度构造，但是该模型很难获得与使用。径向基函数（Radial Basis Function,RBF）是最常用的代理模型方法之一，其优点在于对于高阶非线性的飞行器高精度分析模型，径向基函数在全局近似精度较高，并且随着试验设计样本的增加，所构造的径向基代理模型的近似精度会提高；在试验设计样本附近，近似精度较高。按照代理模型在优化过程中的使用方式可以分为静态代理模型和动态代理模型。静态代理模型在优化前构造完毕，在优化过程中代理模型保持不变，而动态代理模型在优化过程中根据已知信息逐步进行更新直至优化收敛。与静态代理模型相比，动态代理模型在优化效率和结果精度方面更具有优势。

设计空间是飞行器升力面气动、热、结构多学科耦合设计优化中较重要的一环。一方面，设计空间过小可能使优化设计快速地收敛到局部最优解；另一方面，设计空间过大可能使优化很难收敛，甚至不能收敛。因此，从较大的设计空间放缩到合适的包括全局最优的设计空间是飞行器升力面气动、热、结构多学科耦合设计优化中重要的一步。

径向基函数（Radial Basis Function,RBF）是最常用的代理模型方法之一，其优点在于对于高阶非线性的飞行器高精度分析模型，径向基函数在全局近似精度较高；并且随着试验设计样本的增加，所构造的径向基代理模型的近似精度会提高；在试验设计样本附近，近似精度较高。按照代理模型在优化过程中的使用方式可以分为静态代理模型和自适应代理模型。静态代理模型是通过一次试验设计采取足够多的试验设计样本，然后构造代理模型，在优化过程中代理模型保持不变；而自适应代理模型是在优化设计过程序列采取试验设计样本，然后在每次优化迭代过程中根据已知信息逐步改进和更新代理模型，直至优化收敛。与静态代理模型相比，自适应代理模型在优化效率和结果精度方面更具有优势。

基于代理模型的优化技术在航空航天领域具有很广的应用前景，但是国内研究开展较晚。国外研究机构、工业界和商业软件公司纷纷给予该类算法开发了高效的优化器并用于飞行器优化设计，例如播音公司的Boeing探索器，Altair公司Hyperstudy中的ARSM优化器，以及美国Sandia国家研究实验室的DAKOTA等。

径向基函数（RBF）代理模型的基本形式为：

其中，f(x)是代理模型的响应值；B为权向系数向量，B=[β₁,…,β_m]^T且B应满足公式（2）；β_i是第i个权向系数，1≤i≤m；m是试验设计样本个数；是基函数向量，φ是基函数，||·||是二范数，x是预测点；x_i是第i个试验设计样本。

f_i=y_i             （2）

其中，f_i为第i个试验设计样本的径向基函数代理模型预测值；y_i为第i个试验设计样本所对应的分析模型响应值。于是有：

AB=Y              （3）

其中，Y为试验设计样本所对应的分析模型响应值向量。

由公式（3）可得公式（4）：

B=A^-1Y            （4）

本发明采用的另外一种重要技术是模糊聚类法。

在缩减设计空间时，模糊聚类法用于对网格点的隶属度进行分类，完成数据的模糊聚类划分，通过聚类的放缩与删减，最终得到全局最优解所在的重点区域。模糊c-均值聚类方法(Fuzzy c-means Clustering Method,FCM)是由Bezkek于1981年提出的,它是目前广泛采用的一种模糊聚类算法。考虑一个含有m个网格点的集合n_dim代表网格点的维数，此集合依据一定的准则用模糊聚类的方法分成c(2≤c≤c_max)个模糊子集，这里c是指一给定的聚类个数，c_max为聚类个数的最大值，所用的准则是优化一个用来表征聚类的性能指标的目标函数J_con(U,v)，其基本形式为

$\begin{array}{c}{J}_{\mathrm{con}}(U,v)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{con}}{\left(d_{\mathrm{ik}}\right)}^2\\ s.t.\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ik}}=1,{\mathrm{\μ}}_{\mathrm{ik}}\∈\left[\mathrm{0,1}\right]\end{array}---\left(5\right)$

其中，U为隶属度矩阵，U=[u_ik]；u_ik表示第k个网格点属于第i个聚类的隶属度；v=(v₁,v₂,…,v_c)，v_i是第i个聚类中心,1≤i≤c；con是一个大于1的常数（通常取2）；d_ik是欧几里德距离，其数学表达式如公式（6）所示。

d_ik=||x_k-v_i||           （6）

其中，x_k为第k个网格点。

使用拉格朗日乘子法将式（5）转化为无约束优化问题，如式（7）所示。

${\stackrel{\‾}{J}}_{\mathrm{con}}(U,v)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{con}}{\left(d_{\mathrm{ik}}\right)}^2+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_k(\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ik}}-1)---\left(7\right)$

其中，为无约束的目标函数；λ_k为第k个网格点对应的拉格朗日乘子，求公式（7）的极小化问题，得公式（8）：

$\underset{U,v}{\min }{\stackrel{\‾}{J}}_{\mathrm{con}}(U,v)---\left(8\right)$

由公式（8）可得：

$\begin{array}{c}{\mathrm{\μ}}_{\mathrm{ik}}=\frac{1}{\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{ik}}}{d_{\mathrm{jk}}}\right)}^{\frac{2}{\mathrm{con}-1}}}\\ v_i=\frac{1}{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{con}}}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{con}}{x}_k\end{array}---\left(9\right)$

对于高维设计问题（设计变量大于3），FCM方法在相似度与空间距离方面有很大的缺陷，由Bezdek和Dunn于1979年提出的一种基于最大似然估计距离的模糊聚类方法Gath-Geva clustering method（GGM）用来解决高维设计变量中优化设计的不足。该聚类与FCM聚类方式唯一不同的是距离函数不再是欧几里德距离，而是最大似然估计距离，具体形式为：

$d_{\mathrm{ik}}=\frac{\sqrt{\det \left(F_{\mathrm{\ωi}}\right)}}{{\mathrm{\α}}_i}\exp \left(\frac{1}{2}{(x_k-v_i)}^T{F_{\mathrm{\ωi}}}^{-1}(x_k-v_i)\right)---\left(10\right)$

其中，det(F_ωi)表示对F_ωi取行列式；F_ωi为第i个聚类的协方差，1≤i≤c，F_ωi的值由公式（11）可得；exp()表示底数为e的指数函数。

$F_{\mathrm{\ωi}}=\frac{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ik}}(x_k-v_i){(x_k-v_i)}^T}{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{\ω}}}---\left(11\right)$

其中，ω=2为加权幂指数。

发明内容

本发明的目的是为了解决传统的全局优化方法对飞行器高精度分析模型进行优化设计过程中存在计算分析耗时的缺点，提出了一种飞行器升力面结构优化设计方法。

本发明的一种飞行器升力面结构优化设计方法，具体实现步骤如下：

步骤1：给定飞行器升力面高精度分析模型、初始设计空间(用符号S⁽⁰⁾表示)、设计变量。

步骤2：利用拉丁超方试验设计方法（LHD）在初始设计空间S⁽⁰⁾中构造N_initial个初始试验设计样本。初始试验设计样本的个数N_initial由公式（12）确定，用符号N记录迭代次数，此时设置迭代次数N＝1。

$\left\{\begin{array}{c}{N}_{\mathrm{initial}}=6(n\≤3)\\ N_{\mathrm{initial}}=2n(n>3)\end{array}\right.---\left(12\right)$

其中，n表示设计变量的个数，即设计空间的维数。

步骤3：通过调用步骤1给定的飞行器升力面高精度分析模型，计算/仿真当前试验设计样本所对应的飞行器升力面高精度分析模型的响应值，并将试验设计样本及其相对应的响应值存储到试验设计样本数据库中；所述试验设计样本数据库中包含的内容包括：通过拉丁超方试验设计方法所得到的试验设计样本及其相对应的飞行器升力面高精度分析模型响应值、每次迭代过程中的全局近似最优解及其相对应的飞行器升力面高精度分析模型响应值。

步骤4：利用步骤3中所述试验设计样本数据库中的所有试验设计样本及其相对应的飞行器升力面高精度分析模型响应值、每次迭代过程中的全局近似最优解及其相对应的飞行器升力面高精度分析模型响应值构造径向基函数代理模型。

步骤5：在第N次迭代（当前）的设计空间S^(N-1)中，使用拉丁超方试验设计方法随机生成m_{grid_points}个网格点，网格点个数m_{grid_points}由公式（13）确定。用步骤4得到的径向基函数代理模型计算所有网格点对应的径向基函数代理模型响应值。给定整体空间缩减率M₁%，将网格点对应的径向基函数代理模型响应值中较小的M₁%网格点保留，删除其余网格点，从而完成对设计空间的第一次伪缩减。

m_{grid_points}=60n              （13）

步骤6：给定聚类空间个数c，对经过步骤5的处理后保留下的网格点使用模糊聚类法进行搜索与分类，具体为：使用公式（14）对经过步骤5的处理后保留下的网格点进行聚类，依次得到第k个网格点属于第i个聚类的隶属度u_ik，根据隶属度将该网格点分成c个聚类空间(i＝1,2,…,c)，i=1,2,…,c。

$\begin{array}{c}{\mathrm{\μ}}_{\mathrm{ik}}=\frac{1}{\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{ik}}}{d_{\mathrm{jk}}}\right)}^{\frac{2}{\mathrm{con}-1}}}\\ v_i=\frac{1}{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{con}}}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{con}}{x}_k\end{array}---\left(14\right)$

其中，d_id为距离函数，con是一个大于1的常数，v_i为第i个聚类的距离中心，m为网格点的个数， 表示向上取整，x_k为第k个网格点。

当设计变量小于等于3时，距离函数d_ik通过公式（15）得到：

d_ik=||x_k-v_i||            （15）

当设计变量大于3时，距离函数d_ik通过公式（16）得到：

$d_{\mathrm{ik}}=\frac{\sqrt{\det \left(F_{\mathrm{\ωi}}\right)}}{{\mathrm{\α}}_i}\exp \left(\frac{1}{2}{(x_k-v_i)}^T{F_{\mathrm{\ωi}}}^{-1}(x_k-v_i)\right)---\left(16\right)$

其中，det(F_ωi)表示对F_ωi取行列式；F_ωi为第i个聚类的协方差，1≤i≤c，F_ωi的值由公式（17）可得；exp()表示底数为e的指数函数。

$F_{\mathrm{\ωi}}=\frac{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ik}}(x_k-v_i){(x_k-v_i)}^T}{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left({\mathrm{\μ}}_{\mathrm{ik}}\right)}^{\mathrm{\ω}}}---\left(17\right)$

其中，ω=2为加权幂指数。

步骤7：给定局部空间缩减率M₂%，依次将c个聚类空间中径向基函数代理模型响应值较小的M₂%个网格点保留，删除其余网格点，完成对设计空间的第二次伪缩减，用(i＝1,2,…,c)表示每个聚类空间中删除(1-M₂%)个网格点后得到的c个新聚类空间。

步骤8：采用遗传算法在每个新聚类空间对步骤4中所述径向基函数代理模型进行优化，得到第N次迭代中各个新聚类空间的近似最优点（用(i＝1,2,…,c)表示）。将对应的径向基函数代理模型响应值中的最小值对应的近似最优点作为寻优过程中的当前全局近似最优解（用表示）。

步骤9：根据步骤8中所述与通过公式（18）定义两种欧氏距离，分别称为第一欧氏距离（用D₀表示）和第二欧氏距离（用D_i表示），由公式（19）得到第N次迭代的缩减因子（用λ表示）；由公式（20）得到第N次迭代的控制因子（用σ表示）。用λ₀表示初始缩减因子，其为人为设定值，λ₀在区间[0,1]中；用σ₀表示初始控制因子；如果λ<λ₀并且σ>σ₀，表明当前设计空间S^(N-1)需要缩减，进行步骤10的操作；否则，进行步骤12的操作。

$\left\{\begin{array}{c}{D}_0=\left|\right|{\hat{x}}^{^\left(N\right)}-x^{(N-1)}\left|\right|\\ D_i=\left|\right|{\hat{x}}_i^{^(N-1)}-x^{\left(N\right)}\left|\right|\end{array}\right.---\left(18\right)$

其中，||·||是二范数；i＝1,2,…,c。

$\mathrm{\&lambda;}=\frac{D_0}{L_{\max }}---\left(19\right)$

其中，L_max为S^(N-1)的最大边长，设计空间S^(N-1)被认为是超矩形。

$\mathrm{\&sigma;}=\frac{V}{V_0}---\left(20\right)$

其中，V为空间S^(N-1)的超体积，V₀为初始空间S⁽⁰⁾的超体积。

步骤10：对公式（18）中距离D_i(i＝1,2,…,c)从小到大排序，得到新的距离D′_i，如式（21）所示；然后给定需删除聚类数（用n_del表示，n_del<c），则距离D′_i对应的聚类空间(i=c,c-1,…,c-n_del+1)则被删除。

D′_i=sort(D_i)              （21）

其中，i＝1,2,…,c，sort()表示从小到大排序函数。

步骤11：删除当前设计空间S^(N-1)，重组新的设计空间S^(N)。将经过步骤10的操作后剩下的所有网格点的各维的最小值和最大值分别作为新的设计空间S^(N)各维的最小值和最大值。

步骤12：获取径向基函数代理模型的当前全局近似最优解对应的飞行器升力面高精度分析模型响应值（用y^*(N)表示），并将此全局近似最优解及其飞行器升力面高精度分析模型响应值y^*(N)添加到试验设计样本数据库中，并设置迭代次数N自增1。

步骤13：通过公式（22）判断飞行器升力面高精度分析模型响应值y^*(N)是否收敛，如果收敛则将步骤12中所述全局近似最优解作为飞行器升力面高精度分析模型的全局最优解输出，结束整个优化过程；否则，重复步骤4至步骤13，直到找到飞行器升力面高精度分析模型的全局最优解。

$\left\{\begin{array}{c}\left|\frac{y^{*\left(N\right)}-y^{*(N-1)}}{y^{*(N-1)}}\right|<{\mathrm{\&epsiv;}}_1\\ \left|\frac{y^{*\left(N\right)}-{\hat{y}}^{*\left(N\right)}}{y^{*\left(N\right)}}\right|<{\mathrm{\&epsiv;}}_2\end{array}\right.---\left(22\right)$

其中，为第N次全局近似最优解对应的径向基函数代理模型响应值；ε₁、ε₂为人为确定的小量，ε₁取[0.0001,0.01]区间中的值，ε₂=10ε₁。

有益效果

本发明提出的飞行器升力面结构优化设计方法克服了传统的全局优化方法在求解飞行器升力面高精度分析模型设计问题时存在的计算耗时的缺点，能够有效的降低计算成本，提高优化效率，有助于缩短飞行器升力面高精度分析模型设计的周期。用更少的调用飞行器升力面高精度分析模型的次数得到精度更高的全局最优解。

附图说明

图1为本发明具体实施例中的一种飞行器升力面结构优化设计方法的操作流程图；

图2为本发明具体实施例中的飞行器升力面外形布局图；

图3为本发明具体实施例中的飞行器升力面内部结构图及设计变量分布图。

具体实施方式

为了更好的说明本发明的技术方案，下面结合附图和具体实施例，对本发明做进一步说明，并通过与传统飞行器升力面结构优化设计方法结果比较，对本发明的综合性能进行验证分析。

本实施例优化设计的目标是在满足颤振速度约束条件下飞行器升力面质量最小。对于高速飞行器来讲，颤振速度是一个很重要的非线性指标，当飞行器飞行速度达到颤振临界速度时将会引起非线性自激振动，从而使飞行器遭到破坏。

一种飞行器升力面结构优化设计方法，其操作流程如图1所示，具体为：

步骤1：给定飞行器升力面高精度分析模型、初始设计空间(用符号S⁽⁰⁾表示)、设计变量。

本实施例中，优化模型是低展弦比飞行器升力面，其气动外形如图2所示。其中，翼根长4米,翼梢长1.8米,前缘后掠角为34度，后缘后掠角为18度。内部结构分为蒙皮、翼肋与梁，将它们的厚度分别设为三个设计变量，范围为[0.5,3]毫米，建立有限元模型，共有956个壳单元和540个节点，如图3所示。其中，设计变量1为蒙皮厚度，设计变量2为翼肋厚度，设计变量3为梁厚度。飞行器结构材料使用铝，材料参数如表1示。

表1铝材料参数

结构质量由公式（23）得到：

$\mathrm{mass}=\underset{i=1}{\overset{956}{\mathrm{\&Sigma;}}}{S}_i{\mathrm{\&tau;}}_i\mathrm{\&rho;}---\left(23\right)$

其中，mass结构质量；为S_i指第i个单元的面积；τ_i为第i个单元的厚度；ρ为材料密度。

颤振速度可由公式（24）得到：

$\left[M^{*}\right]\left\{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&xi;}}\right\}+\left[K^{*}\right]\left\{\mathrm{\&xi;}\right\}=\frac{1}{2}{\mathrm{\&rho;V}}^2\left[A^{*}\right]\left\{\mathrm{\&xi;}\right\}---\left(24\right)$

其中，[M^*]为模态质量矩阵；{ξ}为模态坐标；[K^*]为模态刚度矩阵；ρ为空气密度；V为颤振速度；[A^*]为广义的模态坐标气动力系数。假定临界颤振速度为1800m/s。

将飞行器升力面高精度分析模型转化为如下的数学模型：

min    mass

st.V_f≥1800m/s                          （25）

0.5mm≤x₁、x₂ x₃≤3mm

其中，V_f为颤振速度，x₁，x₂，x₃为飞行器升力面蒙皮、翼肋和梁的厚度，为优化模型的设计变量。通过拉格朗日乘子法，将公式（25）转化为无约束优化问题如公式（26）所示。

YI=FI+0.5M((max(0,GI+μ/M))²-(μ/M)²)    （26）

其中，FI,GI分别是结构质量mass与颤振速度约束方程(1800-V_f)的归一化。

步骤2：因为设计变量个数n＝3，故在初始设计空间S⁽⁰⁾中，使用拉丁超方试验设计方法生成初始样本点个数N_initial＝6，并设置迭代次数N＝1。

步骤3：通过调用步骤1给定的飞行器升力面高精度分析模型，计算/仿真当前试验设计样本所对应的飞行器升力面高精度分析模型的响应值，并将试验设计样本及其相对应的响应值存储到试验设计样本数据库中。计算一次该飞行器升力面高精度分析模型大约需要花50秒。

步骤4：利用步骤3中所述试验设计样本数据库中的所有试验设计样本及其相对应的飞行器升力面高精度分析模型响应值、每次迭代过程中的全局近似最优解及其相对应的飞行器升力面高精度分析模型响应值构造径向基函数代理模型。

步骤5：在第N次设计空间S^(N-1)中，使用拉丁超方试验设计方法随机生成60×3=180个网格点。用步骤4得到的径向基函数代理模型计算所有网格点对应的径向基函数代理模型响应值。给定整体空间缩减率M₁%=70%，将网格点对应的径向基函数代理模型响应值中较小的70%网格点保留，删除其余网格点，从而完成对设计空间的第一次伪缩减。

步骤6：给定聚类空间个数c＝3，对经过步骤5的处理后保留下的网格点使用模糊聚类法进行搜索与分类，具体为：使用公式（14）对经过步骤5的处理后保留下的网格点进行聚类，依次得到第k个网格点属于第i个聚类的隶属度u_ik，根据隶属度将该网格点分成3个聚类空间(i＝1,2,3)

步骤7：给定局部空间缩减率M₂%=70%，依次将3个聚类空间中径向基函数代理模型响应值较小的70%个网格点保留，删除其余网格点，完成对设计空间的第二次伪缩减，得到3个新的聚类空间(i＝1,2,3)。

步骤8：采用遗传算法在每个新聚类空间对步骤4中所述径向基函数代理模型进行优化，得到第N次迭代中各个新聚类空间的近似最优点(i＝1,2,3)。将对应的径向基函数代理模型响应值中的最小值对应的近似最优点作为寻优过程中的当前全局近似最优解

步骤9：根据步骤8中所述与通过公式（18）定义两种欧氏距离，分别称为第一欧氏距离D₀和第二欧氏距离D_i，由公式（19）得到第N次迭代的缩减因子λ表示；由公式（20）得到第N次迭代的控制因子σ表示。给定初始缩减因子λ₀=0.1与初始控制因子σ₀=0.001，如果λ<λ₀并且σ>σ₀满足，表明当前设计空间S^(N-1)需要缩减，进行步骤10的操作；否则，进行步骤12的操作。

步骤10：对公式（18）中距离D_i(i＝1,2,…,c)从小到大排序，得到新的距离D′_i，如式（21）所示；给定需要删除的聚类个数n_del=1，因此，D_s=max(D₁,D₂,D₃)对应的即被删除。

步骤11：删除当前设计空间S^(N-1)，重组新的设计空间S^(N)。将经过步骤10的操作后剩下的所有网格点的各维的最小值和最大值分别作为新的设计空间S^(N)各维的最小值和最大值。

步骤12：获取径向基函数代理模型的当前全局近似最优解对应的飞行器升力面高精度分析模型响应值y^*(N)，并将此全局近似最优解及其飞行器升力面高精度分析模型响应值y^*(N)添加到试验设计样本数据库中，并设置迭代次数N自增1。

步骤13：通过公式（22）判断飞行器升力面高精度分析模型响应值y^*(N)是否收敛，如果收敛则将步骤12中所述全局近似最优解作为飞行器升力面高精度分析模型的全局最优解输出，结束整个优化过程；否则，重复步骤4至步骤13，直到找到飞行器升力面高精度分析模型的全局最优解。公式（22）中的ε₁=0.01，ε₂=0.1。

经过上述步骤的操作，即完成对低展弦比飞行器升力面的优化，为了保证设计问题的鲁棒性，对该设计问题连续进行了11次优化，取11次中调用分析模型次数中的中间优化结果。同时，采用传统的飞行器优化方法序列二次规划（SQP）和遗传算法（GA）对该飞行器升力面进行优化，其优化结果如表2所示。

表2本发明方法、SQP、GA升力面优化结果对比表

由表2可知，由于该模型的特殊性，序列二次规划（SQP）不能对该有限元模型进行优化，不能进行寻优，颤振速度不能满足设计要求；遗传算法（GA）在调用了758次有限元分析模型后得到了目标结构质量为81.70kg，每次调用分析模型大约需要50s，因此遗传算法优化所消耗的总时间为50×758=37900sec=10.53h。而本发明只调用了69次分析模型，目标结构质量为72.38kg，所消耗的总时间为50×69=3450sec=57.5min。显然，本发明用了大约是遗传算法十分之一的时间，得到了目标更小的结构质量，大大提高了优化效率，节约设计成本。

根据飞行器升力面优化实例分析可见，相比于传统的优化设计方法，本发明有助于提高飞行器升力面优化设计结果与设计质量；另一方面，涉及飞行器升力面高精度分析模型的优化问题，本发明还能大大提高的优化效率，降低飞行器优化设计成本，缩短优化设计周期。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种飞行器升力面结构优化设计方法，其特征在于：其实现过程包括步骤1至步骤13，具体为：

步骤1：给定飞行器升力面高精度分析模型、初始设计空间S⁽⁰⁾、设计变量；

步骤2：利用拉丁超方试验设计方法在初始设计空间S⁽⁰⁾中构造初始试验设计样本；用符号N记录迭代次数，此时设置迭代次数N＝1；

步骤3：通过调用步骤1给定的飞行器升力面高精度分析模型，计算/仿真当前试验设计样本所对应的飞行器升力面高精度分析模型的响应值，并将试验设计样本及其相对应的响应值存储到试验设计样本数据库中；所述试验设计样本数据库中包含的内容包括：通过拉丁超方试验设计方法所得到的试验设计样本及其相对应的飞行器升力面高精度分析模型响应值、每次迭代过程中的全局近似最优解及其相对应的飞行器升力面高精度分析模型响应值；

步骤4：利用步骤3中所述试验设计样本数据库中的所有试验设计样本及其相对应的飞行器升力面高精度分析模型响应值、每次迭代过程中的全局近似最优解及其相对应的飞行器升力面高精度分析模型响应值构造径向基函数代理模型；

步骤5：在第N次迭代的设计空间S^(N-1)中，使用拉丁超方试验设计方法随机生成网格点；用步骤4得到的径向基函数代理模型计算所有网格点对应的径向基函数代理模型响应值；给定整体空间缩减率M₁％，将网格点对应的径向基函数代理模型响应值中较小的M₁％网格点保留，删除其余网格点；

步骤6：给定聚类空间个数c，对经过步骤5的处理后保留下的网格点使用模糊聚类法进行搜索与分类；

步骤7：给定局部空间缩减率M₂％，依次将c个聚类空间中径向基函数代理模型响应值较小的M₂％个网格点保留，删除其余网格点，用表示每个聚类空间中删除(1-M₂％)个网格点后得到的c个新聚类空间；

步骤8：采用遗传算法在每个新聚类空间对步骤4中所述径向基函数代理模型进行优化，得到第N次迭代中各个新聚类空间的近似最优点将对应的径向基函数代理模型响应值中的最小值对应的近似最优点作为寻优过程中的当前全局近似最优解

步骤9：根据步骤8中所述与通过公式(18)定义两种欧氏距离，分别称为第一欧氏距离D₀和第二欧氏距离D_i，由公式(19)得到第N次迭代的缩减因子λ；由公式(20)得到第N次迭代的控制因子σ；用λ₀表示初始缩减因子，其为人为设定值，λ₀在区间[0,1]中；用σ₀表示初始控制因子；如果λ＜λ₀并且σ＞σ₀，表明当前设计空间S^(N-1)需要缩减，进行步骤10的操作；否则，进行步骤12的操作；

$\left\{\begin{array}{c}{D}_0=\left|\right|{\hat{x}}^{\left(N\right)}-{\hat{x}}^{(N-1)}\left|\right|\\ D_i=\left|\right|{\hat{x}}_i^{(N-1)}-{\hat{x}}^{\left(N\right)}\left|\right|\end{array}\right.---\left(18\right)$

其中，||·||是二范数；i＝1,2,…,c；

$\mathrm{\&lambda;}=\frac{D_0}{L_{\max }}---\left(19\right)$

其中，L_max为S^(N-1)的最大边长，设计空间S^(N-1)被认为是超矩形；

$\mathrm{\&sigma;}=\frac{V}{V_0}---\left(20\right)$

其中，V为空间S^(N-1)的超体积，V₀为初始空间S⁽⁰⁾的超体积；

步骤10：对公式(18)中距离D_i(i＝1,2,…,c)从小到大排序，得到新的距离D′_i，如式(21)所示；然后给定需删除聚类数n_del，n_del＜c，距离D′_i对应的聚类空间被删除；

D′_i＝sort(D_i)                         (21)

其中，i＝1,2,…,c，sort()表示从小到大排序函数；

步骤11：删除当前设计空间S^(N-1)，重组新的设计空间S^(N)；将经过步骤10的操作后剩下的所有网格点的各维的最小值和最大值分别作为新的设计空间S^(N)各维的最小值和最大值；

步骤12：获取径向基函数代理模型的当前全局近似最优解对应的飞行器升力面高精度分析模型响应值y^*(N)，并将此全局近似最优解及其飞行器升力面高精度分析模型响应值y^*(N)添加到试验设计样本数据库中，并设置迭代次数N自增1；

步骤13：通过公式(22)判断飞行器升力面高精度分析模型响应值y^*(N)是否收敛，如果收敛则将步骤12中所述全局近似最优解作为飞行器升力面高精度分析模型的全局最优解输出，结束整个优化过程；否则，重复步骤4至步骤13，直到找到飞行器升力面高精度分析模型的全局最优解；

$\left\{\begin{array}{c}\left|\frac{y^{*\left(N\right)}-y^{*(N-1)}}{y^{*(N-1)}}\right|<{\mathrm{\&epsiv;}}_1\\ \left|\frac{y^{*\left(N\right)}-{\hat{y}}^{*\left(N\right)}}{y^{*\left(N\right)}}\right|<{\mathrm{\&epsiv;}}_2\end{array}\right.---\left(22\right)$

其中，为第N次全局近似最优解对应的径向基函数代理模型响应值；ε₁、ε₂为人为确定的小量。

2.如权利要求1所述的一种飞行器升力面结构优化设计方法，其特征在于：其步骤2中所述初始试验设计样本的个数N_initial由公式(12)确定；

$\left\{\begin{array}{c}{N}_{\mathrm{initial}}=6(n\&le;3)\\ N_{\mathrm{initial}}=2n(n>3)\end{array}\right.---\left(12\right)$

其中，n表示设计变量的个数，即设计空间的维数。

3.如权利要求1或2所述的一种飞行器升力面结构优化设计方法，其特征在于：其步骤5中所述网格点个数m_{grid_points}由公式(13)确定；

m_{grid_points}＝60n               (13)

其中，n表示设计变量的个数。

4.如权利要求1或2所述的一种飞行器升力面结构优化设计方法，其特征在于：其步骤6中所述对经过步骤5的处理后保留下的网格点使用模糊聚类法进行搜索与分类的具体方法为：使用公式(14)对经过步骤5的处理后保留下的网格点进行聚类，依次得到第k个网格点属于第i个聚类的隶属度μ_ik，根据隶属度将该网格点分成c个聚类空间 i＝1,2,…,c；

$\left\{\begin{array}{c}{\mathrm{\&mu;}}_{\mathrm{ik}}=\frac{1}{\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(\frac{d_{\mathrm{ik}}}{d_{\mathrm{jk}}}\right)}^{\frac{2}{\mathrm{con}-1}}}\\ v_i=\frac{1}{\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left({\mathrm{\&mu;}}_{\mathrm{ik}}\right)}^{\mathrm{con}}}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left({\mathrm{\&mu;}}_{\mathrm{ik}}\right)}^{\mathrm{con}}{x}_k\end{array}\right.---\left(14\right)$

其中，d_ik为距离函数，con是一个大于1的常数，v_i为第i个聚类的距离中心，m为网格点的个数， 表示向上取整，x_k为第k个网格点；

当设计变量小于等于3时，距离函数d_ik通过公式(15)得到：

d_ik＝||x_k-v_i||                    (15)

当设计变量大于3时，距离函数d_ik通过公式(16)得到：

$d_{\mathrm{ik}}=\frac{\sqrt{\det \left(F_{\mathrm{\&omega;i}}\right)}}{{\mathrm{\&alpha;}}_i}\exp \left(\frac{1}{2}{(x_k-v_i)}^T{F_{\mathrm{\&omega;i}}}^{-1}(x_k-v_i)\right)---\left(16\right)$

其中，det(F_ωi)表示对F_ωi取行列式；F_ωi为第i个聚类的协方差，1≤i≤c，F_ωi的值由公式(17)可得；exp()表示底数为e的指数函数；

$F_{\mathrm{\&omega;i}}=\frac{\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ik}}(x_k-v_i){(x_k-v_i)}^T}{\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left({\mathrm{\&mu;}}_{\mathrm{ik}}\right)}^{\mathrm{\&omega;}}}---\left(17\right)$

其中，ω＝2为加权幂指数。

5.如权利要求1或2所述的一种飞行器升力面结构优化设计方法，其特征在于：其步骤13中所述ε₁取[0.0001,0.01]区间中的值，ε₂＝10ε₁。