CN105046021A - 非定常气动力最小状态有理近似的非线性优化算法 - Google Patents

Abstract

一种非定常气动力最小状态有理近似的非线性优化算法。本发明对D矩阵的关键模态行初始化再顺序求解E矩阵的各列和D矩阵的其余行，并充分利用当前迭代点处的函数值信息，有效提高了逆Hessian阵的近似精度，提高了外层非线性优化算法的效率。本发明有效提高了广义气动力矩阵中关键模态行元素的拟合精度，有利于提高时域颤振分析结果精度和突风响应分析结果精度。将频域下离散的非定常气动力系数矩阵以MS法的形式延拓至拉氏域，且能有效提高计算效率和关键模态项的精度。

Description

非定常气动力最小状态有理近似的非线性优化算法

技术领域

本发明涉及气动弹性力学领域，具体是一种用于频域非定常气动力有理近似的非线性优化方法。

背景技术

为便于气动伺服弹性系统的多学科优化设计，需要将弹性飞行器的运动方程转换为时域的一阶时不变状态空间方程，而结构在任意运动下的非定常气动力模型则是其中一个关键组成部分。在亚音速或超音速情况下，以偶极子格网法为基础计算得到的广义非定常气动力系数矩阵(GAF)是一系列给定减缩频率的函数，代表结构在简谐振荡时所受的广义气动力。为了将广义气动力变换到时域空间，需要将其近似延拓为拉氏域的有理函数形式，再经过整理，结合弹性飞行器的运动方程即可得到气动弹性系统的状态空间模型。

现在工程中最常用的有理函数近似(RFA)方法均以最小二乘法为基础，包括Roger法、修正矩阵(MMP)法和最小状态(MS)法等。从这些方法出发得到的状态空间方程中均包含由气动滞后根产生的气动状态扩充项。Roger法产生的气动扩充项数量为结构模态数与气动滞后根数的乘积，MMP法对应的气动力扩充项数量为广义气动力系数矩阵各列对应的气动滞后根数量之和，而MS法对应的气动力扩充项数等于气动滞后根的数量。大量的工程实践表明，在气动力扩充项的数量相同时，MS法的拟合精度最高，但其中的非线性最小二乘迭代过程计算量很大。文献1“陈青.一种建立非定常气动力频域模型的简单方法[J].空气动力学学报，1988年，第6卷第4期”吸取了Roger近似式中同一矩阵各元素独立确定时精度较高的特点，对选定模态受到的非定常气动力进行了精确拟合，并保留了MS法的形式，改善了拟合精度，但对拟合项无合理加权且对气动滞后根的一维优化方法效果不佳。针对MS法计算量大、初值的选取依赖于工程经验等缺点，文献2“何程.一种改进的非定常气动力拟合方法[J].航空学报，1993年，第14卷第7期”通过调整拟合式最后一项的形式并将其中气动滞后根取为大的互异负实数的方式将问题转化为线性拟合问题，但要求滞后根的数量必须与结构模态数一致且拟合精度不高。文献3“E.Nissim.OntheFormulationofMinimum-StateApproximationasaNonlinearOptimizationProblem1.JournalofAircraft,Vol.43,No.4,2006,pp.1007–1013.”在保留MS法拟合形式的基础上释放了等式约束，并通过解析法得到了自变量的梯度算法，再结合对广义气动力系数矩阵的比例缩放使拟合效率和精度都得到了提高，但非线性优化项的数量远大于气动滞后根的个数，应用难度较大。目前非线性优化算法已经被应用到各类拟合方法中来优化气动滞后根以减小拟合误差，这就进一步对拟合方法和非线性优化算法的效率提出了要求。

发明内容

为克服已有方法精度有限或效率不高的缺点，本发明提出了一种非定常气动力最小状态有理近似的非线性优化算法。

本发明的具体过程是：

步骤1，建立机翼有限元模型和气动力模型，计算机翼在给定减缩频率下的广义气动力系数矩阵，具体是：

Ⅰ建立机翼有限元模型，通过Nastran软件对建立的机翼有限元模型进行模态分析，模态分析中取机翼的前9阶弹性模态。所述建立的机翼有限元模型中，机翼下翼面的弦向剖面共有8个网格节点，分别位于距前缘0.0％、5.0％、10.0％、27.0％、45.5％、63.5％、82.0％和100％处。沿机翼展向各网格节点数量与翼肋的数量相同，并与各翼肋的展向位置一致。

Ⅱ建立气动力模型，并将得到的模态分析结果导入到Zaero软件中，得到机翼在给定减缩频率下的广义气动力系数矩阵：

Q(ik)＝F(ik)+iG(ik)(1)

其中Q(ik)为给定减缩频率下的广义气动力系数矩阵，k＝ωb/V为减缩频率，ω为机翼振荡的圆频率，b为机翼参考弦长，V为来流速度。F(ik)和G(ik)分别表示广义气动力系数矩阵的实部和虚部。

对Q(ik)拟合时采用MS法的拟合公式，即

$Q_{ap}\left(\stackrel{\‾}{s}\right)=F_{ap}\left(ik\right)+{\mathrm{iG}}_{ap}\left(ik\right)=A_0+A_1\stackrel{\‾}{s}+A_2{\stackrel{\‾}{s}}^2+D{(\stackrel{\‾}{s}I-R)}^{-1}E\·\stackrel{\‾}{s}---\left(2\right)$

其中s为拉氏变量。表示拟合得到的广义气动力系数矩阵，F_ap(ik)和G_ap(ik)分别表示的实部和虚部。A₀∈R^n×n、A₁∈R^n×n、A₂∈R^n×n、D∈R^n×m、R∈R^m×m、E∈R^m×n均为待定系数矩阵，n为结构模态数，m为气动滞后根数量。I为m阶单位矩阵。R为由气动滞后根组成的对角矩阵，表示为

R＝diag(x)＝diag([x₁x₂…x_m])(3)

其中x表示气动滞后根组成的向量。

所述建立的气动力模型的弦向被均分为15等份，展向被均分为24等份。

步骤2，给定气动滞后根向量的初值x₀，对MS法的拟合公式进行总体拟合并计算该总体拟合的误差f(x₀)。

在对得到的广义气动力系数矩阵进行拟合前，需先给定气动滞后根向量的初值，初始化方法为：当m＝1时，x₀＝[-0.3]；当m＝2时，x₀＝[-0.3-0.5]；当m＝3时，x₀＝[-0.3-0.5-0.7]；当m＝4时，x₀＝[-0.3-0.5-0.7-0.9]。

根据给定的气动滞后根向量的初值x₀得到由气动滞后根组成的对角矩阵R；对MS法的拟合公式进行总体拟合并计算总体拟合误差f(x₀)，具体过程是：

Ⅰ取结构的第r阶模态为关键模态，并令D矩阵的第r行元素为任意非零常数，其中D矩阵的第r行元素代表结构的第r阶模态对应的广义气动力。此时方程(2)变为

$Q_{ap,r}\left(\stackrel{\‾}{s}\right)=A_{0,r}+A_{1,r}\stackrel{\‾}{s}+A_{2,r}{\stackrel{\‾}{s}}^2+D_r{\left(\stackrel{\‾}{s}I-R\right)}^{-1}E\·\stackrel{\‾}{s}---\left(4\right)$

其中下标r均表示矩阵的第r行。

Ⅱ约束方程在s＝0处的有理逼近值等于气动力系数矩阵列表值，在s＝ik_f处的实部有理逼近值等于矩阵的列表值，在s＝ik_g处的虚部有理逼近值等于矩阵的列表值，其中k_f和k_g均为指定的减缩频率值。之后从D矩阵第r行出发拟合E矩阵各列元素，所述E矩阵第j列的加权最小二乘求解式如下

$\underset{l=1}{\overset{L}{\Σ}}\left(P_l^T{W}_{rj}^2{P}_l+Q_l^T{W}_{rj}^2{Q}_l\right)E_j=\underset{l=1}{\overset{L}{\Σ}}\left\{P_l^T{W}_{rj}^2{\stackrel{\‾}{F}}_j\left({\mathrm{ik}}_l\right)+Q_l^T{W}_{rj}^2{\stackrel{\‾}{G}}_j\left({\mathrm{ik}}_l\right)\right\}---\left(5\right)$

其中E_j表示E矩阵第j列，

$\begin{array}{c}{P}_l=k_l^2{D}_r\[{(k_l^{2}I+R^2)}^{-1}-{(k_f^{2}I+R^2)}^{-1}\]\\ Q_l=k_l{D}_r\[{(k_l^{2}I+R^2)}^{-1}-{(k_g^{2}I+R^2)}^{-1}\]R\end{array}---\left(6\right)$

$W_{rj}=\frac{1}{\underset{l}{max}\left\{\left|Q_{rj}\left({\mathrm{ik}}_l\right)\right|,1\right\}}---\left(7\right)$

其中L为减缩频率个数。Q_rj(ik_l)表示Q(ik)第r行第j列元素在减缩频率k_l处的值，W_rj表示Q(ik)的加权矩阵W的第r行第j列元素。和分别表示F(ik_l)和G(ik_l)的第j列。

因当前机翼在颤振点处的减缩频率为0.07，故取k_f＝k_g＝0.05以使颤振点附近的拟合精度更高。

Ⅲ求解出E矩阵后，再逐行用加权最小二乘法求解D矩阵的各元素，第r行除外。D矩阵第i行的加权最小二乘求解式如下

$\underset{l=1}{\overset{L}{\Σ}}\left(P_l^{*T}{W}_i^2{P}_l^{*}+Q_l^{*T}{W}_i^2{Q}_l^{*}\right)D_i^T=\underset{l=1}{\overset{L}{\Σ}}\left(P_l^{*T}{W}_i^2{\hat{F}}_i\left({\mathrm{ik}}_l\right)+Q_l^{*T}{\hat{G}}_i\left({\mathrm{ik}}_l\right)\right)---\left(8\right)$

其中D_i表示D矩阵第i行，

$\begin{array}{c}{P}_l^{*}=k_l^2{E}^T\[{(k_l^{2}I+R^2)}^{-1}-{(k_f^{2}I+R^2)}^{-1}\]\\ Q_l^{*}=k_l{E}^{T}R\[{(k_l^{2}I+R^2)}^{-1}-{(k_g^{2}I+R^2)}^{-1}\]\end{array}---\left(9\right)$

W_i ²表示加权矩阵W第i行各元素的平方值构成的对角矩阵，和分别表示F(ik_l)和G(ik_l)的第i行；

Ⅳ计算对MS法的拟合公式进行总体拟合结果的总误差f(x₀)

$f\left(x_0\right)=\sqrt{\underset{l=1}{\overset{L}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\underset{i=1}{\overset{n}{\Σ}}{W}_{ij}^2\left\{{\[F_{ap,ij}\left(k_l\right)-F_{ij}\left(k_l\right)\]}^2+{\[G_{ap,ij}\left(k_l\right)-G_{ij}\left(k_l\right)\]}^2\right\}}---\left(10\right)$

其中i和j分别表示各矩阵的行和列。

步骤3，开始对滞后根向量进行非线性优化。

给定初始Hessian阵的逆阵H₀，为m阶单位阵；令迭代下标kk＝0；

步骤4，计算第一个迭代点x₀处的梯度值其第i项的算式如下

${(\▿f_0)}_i=\frac{f(x_0+\mathrm{\α}\·e_i)-f(x_0-\mathrm{\α}\·e_i)}{2\mathrm{\α}},(i=1,2,...m)---\left(11\right)$

其中e_i为第i个元素为1的m阶单位向量，α为一个充分小的实数，此处取为0.001。f(x₀+α·e_i)和f(x₀-α·e_i)的计算方法同步骤2中f(x₀)的计算方法。

步骤5，确定搜索方向。

通过公式(12)确定搜索方向

$d_{kk}=-H_{kk}\▿f_{kk}---\left(12\right)$

步骤6，沿d_kk线性搜索步长因子α_kk。具体过程是：

Ⅰ给定参数0<ξ<0.5和0<β<1。取线性搜索时计算目标函数值的最大许可次数为N₁。令mm＝0代表本次循环目标函数值的计算次数，令nn＝0代表本次循环自标量的越界次数；

Ⅱ令x_kk,st＝x_kk+β^mm+nnd_kk，x_kk,st代表从x_kk出发得到的试探点；

Ⅲ判断或是否满足，其中表示x_kk,st的第i项，x ⁱ和分别表示自变量第i项的下界和上界。若都满足则转至步骤Ⅳ；否则令nn+1→nn并转至步骤Ⅱ,其中“→”表示赋值运算；

Ⅳ计算试探点x_kk,st处的目标函数值f(x_kk,st)；

Ⅴ判断mm是否超过N₁。若不超过，则转到步骤Ⅵ；否则取mm为使f(x_kk,st)最小的数，令步长因子α_kk＝β^mm+nn，f_kk+1＝f(x_kk,st)，停止线性搜索；

Ⅵ判断试探点处Armijo不等式的满足情况

$f\left(x_{kk,st}\right)-f\left(x_{kk}\right)\&le;{\mathrm{\&xi;\&beta;}}^{mm+nn}\&dtri;f_{kk}^T{d}_{kk}---\left(13\right)$

$f\left(x_{kk,st}\right)-f\left(x_{kk}\right)\&GreaterEqual;(1-\mathrm{\&xi;}){\mathrm{\&beta;}}^{mm+nn}\&dtri;f_{kk}^T{d}_{kk}---\left(14\right)$

若满足，则令步长因子α_kk＝β^mm+nn，f_kk+1＝f(x_kk,st)，停止线性搜索；否则令mm+1→mm，转至步骤Ⅴ继续搜索计算，直至mm超过N₁或者在试探点处满足Armijo不等式(13)和(14)。

步骤7，求下一个迭代点x_kk+1：

x_kk+1＝x_kk+α_kkd_kk(15)

步骤8，计算点x_kk+1处的梯度值计算方法同步骤4。

步骤9，判断非线性优化过程是否收敛，具体是：

Ⅰ判断||x_kk||>ε₆是否成立，若不成立则转步骤Ⅱ，否则再判断||x_kk+1-x_kk||/||x_kk||≤ε₁是否成立，若成立则停止优化迭代，否则转步骤Ⅲ；

Ⅱ判断||x_kk+1-x_kk||≤ε₂是否成立，若成立则停止优化迭代，否则转步骤Ⅲ；

Ⅲ判断|f_kk|>ε₇是否成立，若不成立则转步骤Ⅳ，否则再判断|f_kk-f_kk+1|/|f_kk|≤ε₃是否成立，若成立则停止优化迭代，否则转步骤Ⅴ；

Ⅳ判断|f_kk-f_kk+1|≤ε₄是否成立，若成立则停止优化迭代，否则转步骤Ⅴ；

Ⅴ判断是否成立，若成立则停止优化迭代，否则转步骤10，继续进行优化迭代。

在上述的收敛准则中，ε₁、ε₂、ε₃、ε₄、ε₅为误差限，取为10^-6；ε₆、ε₇分别用于判断||x_k||和|f_k|的数量级，取为10^-4。

步骤10，计算H_kk+1。具体是：

Ⅰ令自变量变化为s_kk＝x_kk+1-x_kk，梯度差为

Ⅱ通过改进的BFGS校正公式(16)和(17)计算逆Hessian阵H_kk+1：

$H_{kk+1}=\left\{\begin{array}{cc}\left(I-\frac{s_{kk}{\stackrel{\&OverBar;}{y}}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}}\right)H_{kk}\left(I-\frac{{\stackrel{\&OverBar;}{y}}_{kk}{s}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}}\right)+\frac{s_{kk}{s}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}},& {\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}>0\\ H_{kk},& {\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}\&le;0\end{array}\right.---\left(16\right)$

其中

${\stackrel{\&OverBar;}{y}}_{kk}=y_{kk}+\frac{12(f_{kk}-f_{kk+1})+5\&dtri;f_{kk+1}^T{s}_{kk}+(7-{\mathrm{\&alpha;}}_{kk})\&dtri;f_{kk}^T{s}_{kk}}{s_{kk}^T{\mathrm{\&theta;}}_{kk}}{\mathrm{\&theta;}}_{kk}---\left(17\right)$

式中，表示修正的梯度差，θ_kk∈R^m是满足的任意向量，此处取为s_kk。

步骤11，令迭代下标kk→kk+1，转到步骤5继续迭代，直至满足步骤9中的收敛条件。至此，完成了对NACA0012对称翼型的M6机翼的非定常气动力最小状态有理近似的非线性优化算法。

为求解MS拟合式中的各未知矩阵，传统方法是先对整个D矩阵以规定的方式初始化后再对D矩阵和E矩阵迭代求解直至收敛，这种方法计算效率较低。若在传统MS法的基础上对滞后根进行非线性优化，则整个算法就包含了内外两层迭代过程，其计算时间将大大超过其他方法。而本发明在步骤2中则解决了这个问题，先对D矩阵的关键模态行初始化再顺序求解E矩阵的各列和D矩阵的其余行，其优点有两个方面：首先是避免了传统方法中对D矩阵和E矩阵的迭代过程，节省了计算时间，有利于外层非线性优化过程的实施；其次是有效提高了广义气动力矩阵中关键模态行元素的拟合精度，这有利于提高时域颤振分析结果精度和突风响应分析结果精度。另外，本发明在步骤10中采用了改进的BFGS校正公式进行改进，充分利用了当前迭代点处的函数值信息，有效提高了逆Hessian阵的近似精度，这有利于提高外层非线性优化算法的效率。

图6给出了当m＝1时广义气动力矩阵中关键模态行两个元素的拟合结果，图中共比较了五种方法，分别为MS-Dr法、MS法、MS-opt法、Roger法和Roger-opt法。其中“MS-Dr法”为本发明所使用的方法，表明本发明采用MS法的拟合形式，且起始于对D矩阵第r行的初始化。“MS-opt法”表示采用现有技术的非线性算法对MS法的气动滞后根进行优化的方法。“Roger-opt法”表示采用现有技术的非线性算法对Roger法的气动滞后根进行优化的方法。图中“Real”表示被拟合项的离散列表值。

为便于对比，“MS-opt法”和“Roger-opt法”中的非线性优化算法均与本发明中的优化算法框架一致。从图中可以看出，在未进行非线性优化时，Roger方法和MS法都出现了误差较大的情况，而在进行非线性优化后两者的拟合精度都有明显的改善，从而验证了本发明中非线性优化算法的效果。另外，本发明的精度接近优化后的Roger方法，但因其保留了MS法的形式，因此由其出发得到的状态空间模型能够保持MS法模型阶数低的优点。

本发明能够将频域下离散的非定常气动力系数矩阵以MS法的形式延拓至拉氏域，且能有效提高计算效率和关键模态项的精度。

附图说明

图1是本发明的具体实施过程示意图；

图2是机翼的整体有限元模型；

图3是机翼前梁腹板、后梁腹板和翼肋的有限元模型；

图4是机翼有限元模型的弦向剖面；

图5是机翼的气动力模型；

图6是部分元素的拟合结果，其中，6a是Q(ik)矩阵第一行第一列元素的拟合效果，6b是第一行第九列元素的拟合效果；

图7是本发明的流程图。图中：

1.实部；2.虚部；3.MS-Dr法；4.MS法；5.MS-opt法；6.Roger法；7.Roger-opt法；8.Real；9.翼肋；10.前梁腹板；11.蒙皮；12.后梁腹板。

具体实施方式

本实施例是NACA0012对称翼型的M6机翼的非定常气动力最小状态有理近似的非线性优化算法。所述机翼根部完全固支。机翼的参数如下：根弦长为0.8139m，翼尖弦长为0.4573m，展长为1.1963m，前缘后掠角为30°，后缘后掠角为15.8°，无扭转角，蒙皮厚度为0.002m，梁腹板厚度为0.0015m，翼肋厚度为0.0015m，材料为铝合金，弹性模量E＝70GPa，泊松比μ＝0.30，密度为ρ＝2.7×10³kg/m³。

本实施例的具体过程是：

步骤1，建立机翼有限元模型和气动力模型，计算机翼在给定减缩频率下的广义气动力系数矩阵，具体是：

Ⅰ建立机翼有限元模型：在Patran软件中建立机翼有限元模型并用Nastran软件进行模态分析。采用三角形和四边形壳单元进行有限元建模，图2为机翼的整体有限元模型，图3为机翼前梁腹板、后梁腹板和翼肋的有限元模型。机翼由翼肋9、梁腹板和蒙皮11三部分组成，所述的梁腹板包括前梁腹板10和后梁腹板12。所述11个翼肋沿展向均匀分布。所述前梁腹板10位于距离机翼前缘27.0％处，后梁腹板12位于距离机翼前缘63.5％处。机翼下翼面的弦向剖面共有8个网格节点Z_i，i＝1～8。如图4所示，8个网格节点分别记为Z₁、Z₂、Z₃、Z₄、Z₅、Z₆、Z₇、Z₈，分别位于距前缘0.0％、5.0％、10.0％、27.0％、45.5％、63.5％、82.0％和100％处。沿机翼展向各网格节点数量与翼肋的数量相同，并与各翼肋的展向位置一致。通过Nastran软件对建立的机翼有限元模型进行模态分析，模态分析中取机翼的前9阶弹性模态。

Ⅱ建立气动力模型：在Zaero软件中建立机翼气动力模型并计算广义气动力系数矩阵。机翼的气动力模型如图5所示，气动力模型的弦向被均分为15等份，展向被均分为24等份。将通过Nastran软件得到的模态分析结果导入到Zaero中，得到机翼在给定减缩频率下的广义气动力系数矩阵：

Q(ik)＝F(ik)+iG(ik)(1)

其中Q(ik)为给定减缩频率下的广义气动力系数矩阵，k＝ωb/V为减缩频率，ω为机翼振荡的圆频率，b为机翼参考弦长，V为来流速度。F(ik)和G(ik)分别表示广义气动力系数矩阵的实部和虚部。

本实施例中，在获得所述的广义气动力系数矩阵时，设定马赫数为0.7，减缩频率为k＝0.0、0.05、0.1、0.15、0.2、0.25、0.3、0.35、0.4、0.5、0.6、0.7、0.8、0.9、1.0。

对Q(ik)拟合时采用MS法的拟合公式，即

$Q_{ap}\left(\stackrel{\&OverBar;}{s}\right)=F_{ap}\left(ik\right)+{\mathrm{iG}}_{ap}\left(ik\right)=A_0+A_1\stackrel{\&OverBar;}{s}+A_2{\stackrel{\&OverBar;}{s}}^2+D{(\stackrel{\&OverBar;}{s}I-R)}^{-1}E\&CenterDot;\stackrel{\&OverBar;}{s}---\left(2\right)$

其中s为拉氏变量。表示拟合得到的广义气动力系数矩阵，F_ap(ik)和G_ap(ik)分别表示的实部和虚部。A₀∈R^n×n、A₁∈R^n×n、A₂∈R^n×n、D∈R^n×m、R∈R^m×m、E∈R^m×n均为待定系数矩阵，n为结构模态数，m为气动滞后根数量。I为m阶单位矩阵。R为由气动滞后根组成的对角矩阵，可表示为

R＝diag(x)＝diag([x₁x₂…x_m])(3)

其中x表示气动滞后根组成的向量。

步骤2，给定气动滞后根向量的初值x₀，对MS法的拟合公式进行总体拟合并计算该总体拟合的误差f(x₀)。

在对得到的广义气动力系数矩阵进行拟合前，需先给定气动滞后根向量的初值，初始化方法为：当m＝1时，x₀＝[-0.3]；当m＝2时，x₀＝[-0.3-0.5]；当m＝3时，x₀＝[-0.3-0.5-0.7]；当m＝4时，x₀＝[-0.3-0.5-0.7-0.9]。

根据给定的气动滞后根向量的初值x₀得到由气动滞后根组成的对角矩阵R；对MS法的拟合公式进行总体拟合并计算总体拟合误差f(x₀)，具体过程是：

Ⅰ取结构的第r阶模态为关键模态，并令D矩阵的第r行元素为任意非零常数，其中D矩阵的第r行元素代表结构的第r阶模态对应的广义气动力。此时方程(2)变为

$Q_{ap,r}\left(\stackrel{\&OverBar;}{s}\right)=A_{0,r}+A_{1,r}\stackrel{\&OverBar;}{s}+A_{2,r}{\stackrel{\&OverBar;}{s}}^2+D_r{\left(\stackrel{\&OverBar;}{s}I-R\right)}^{-1}E\&CenterDot;\stackrel{\&OverBar;}{s}---\left(4\right)$

其中下标r均表示矩阵的第r行。

因在Zaero软件中用P-K法进行频域颤振分析时发现，颤振模态主要由一阶弯曲模态产生，故取一阶模态为关键模态，即r＝1，因为此模态不仅对结构的稳定性影响最大，其对结构的响应特性也有显著影响。

Ⅱ约束方程在s＝0处的有理逼近值等于气动力系数矩阵列表值，在s＝ik_f处的实部有理逼近值等于矩阵的列表值，在s＝ik_g处的虚部有理逼近值等于矩阵的列表值，其中k_f和k_g均为指定的减缩频率值。之后从D矩阵第r行出发拟合E矩阵各列元素，其中E矩阵第j列的加权最小二乘求解式如下

$\underset{l=1}{\overset{L}{\&Sigma;}}\left(P_l^T{W}_{rj}^2{P}_l+Q_l^T{W}_{rj}^2{Q}_l\right)E_j=\underset{l=1}{\overset{L}{\&Sigma;}}\left\{P_l^T{W}_{rj}^2{\stackrel{\&OverBar;}{F}}_j\left({\mathrm{ik}}_l\right)+Q_l^T{W}_{rj}^2{\stackrel{\&OverBar;}{G}}_j\left({\mathrm{ik}}_l\right)\right\}---\left(5\right)$

其中E_j表示E矩阵第j列，

$\begin{array}{c}{P}_l=k_l^2{D}_r\&lsqb;{(k_l^{2}I+R^2)}^{-1}-{(k_f^{2}I+R^2)}^{-1}\&rsqb;\\ Q_l=k_l{D}_r\&lsqb;{(k_l^{2}I+R^2)}^{-1}-{(k_g^{2}I+R^2)}^{-1}\&rsqb;R\end{array}---\left(6\right)$

$W_{rj}=\frac{1}{\underset{l}{max}\left\{\left|Q_{rj}\left({\mathrm{ik}}_l\right)\right|,1\right\}}---\left(7\right)$

其中L为减缩频率个数。Q_rj(ik_l)表示Q(ik)第r行第j列元素在减缩频率k_l处的值，W_rj表示Q(ik)的加权矩阵W的第r行第j列元素。和分别表示F(ik_l)和G(ik_l)的第j列。

因当前机翼在颤振点处的减缩频率为0.07，故取k_f＝k_g＝0.05以使颤振点附近的拟合精度更高。

Ⅲ求解出E矩阵后，再逐行用加权最小二乘法求解D矩阵的各元素，第r行除外。D矩阵第i行的加权最小二乘求解式如下

$\underset{l=1}{\overset{L}{\&Sigma;}}\left(P_l^{*T}{W}_i^2{P}_l^{*}+Q_l^{*T}{W}_i^2{Q}_l^{*}\right)D_i^T=\underset{l=1}{\overset{L}{\&Sigma;}}\left(P_l^{*T}{W}_i^2{\hat{F}}_i\left({\mathrm{ik}}_l\right)+Q_l^{*T}{W}_i^2{\hat{G}}_i\left({\mathrm{ik}}_l\right)\right)---\left(8\right)$

其中D_i表示D矩阵第i行，

$\begin{array}{c}{P}_l^{*}=k_l^2{E}^T\&lsqb;{\left(k_l^{2}I+R^2\right)}^{-1}-{\left(k_f^{2}I+R^2\right)}^{-1}\&rsqb;\\ Q_l^{*}=k_l{E}^{T}R\&lsqb;{\left(k_l^{2}I+R^2\right)}^{-1}-{\left(k_g^{2}I+R^2\right)}^{-1}\&rsqb;\end{array}---\left(9\right)$

W_i ²表示加权矩阵W第i行各元素的平方值构成的对角矩阵，和分别表示F(ik_l)和G(ik_l)的第i行；

Ⅳ计算对MS法的拟合公式进行总体拟合结果的总误差f(x₀)

$f\left(x_0\right)=\sqrt{\underset{l=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{n}{\&Sigma;}}\underset{i=1}{\overset{n}{\&Sigma;}}{W}_{ij}^2\left\{{\&lsqb;F_{ap,ij}\left(k_l\right)-F_{ij}\left(k_l\right)\&rsqb;}^2+{\&lsqb;G_{ap,ij}\left(k_l\right)-G_{ij}\left(k_l\right)\&rsqb;}^2\right\}}---\left(10\right)$

其中i和j分别表示各矩阵的行和列。

步骤3，开始对滞后根向量进行非线性优化。

给定初始Hessian阵的逆阵H₀，为m阶单位阵；令迭代下标kk＝0；

步骤4，计算第一个迭代点x₀处的梯度值其第i项的算式如下

${(\&dtri;f_0)}_i=\frac{f(x_0+\mathrm{\&alpha;}\&CenterDot;e_i)-f(x_0-\mathrm{\&alpha;}\&CenterDot;e_i)}{2\mathrm{\&alpha;}},(i=1,2,...m)---\left(11\right)$

其中e_i为第i个元素为1的m阶单位向量，α为一个充分小的实数，此处取为0.001。f(x₀+α·e_i)和f(x₀-α·e_i)的计算方法同步骤2中f(x₀)的计算方法。

步骤5，确定搜索方向

通过公式(12)确定搜索方向

$d_{kk}=-H_{kk}\&dtri;f_{kk}---\left(12\right)$

步骤6，沿d_kk线性搜索步长因子α_kk，具体过程是：

Ⅰ给定参数0<ξ<0.5和0<β<1，所述ξ和β分别为0.3和0.7。取线性搜索时计算目标函数值的最大许可次数为N₁，本实施例中，所述N₁为20。令mm＝0代表本次循环目标函数值的计算次数，令nn＝0代表本次循环自标量的越界次数；

Ⅱ令x_kk,st＝x_kk+β^mm+nnd_kk，x_kk,st代表从x_kk出发得到的试探点；

Ⅲ判断或是否满足，其中表示x_kk,st的第i项，x ⁱ和分别表示自变量第i项的下界和上界，在此分别取为-3.0和-0.1。若都满足则转至步骤Ⅳ；否则令nn+1→nn并转至步骤Ⅱ,其中“→”表示赋值运算；

Ⅳ计算试探点x_kk,st处的目标函数值f(x_kk,st)；

Ⅴ判断mm是否超过N₁。若不超过，则转到步骤Ⅵ；否则取mm为使f(x_kk,st)最小的数，令步长因子α_kk＝β^mm+nn，f_kk+1＝f(x_kk,st)，停止线性搜索；

f.判断试探点处Armijo不等式的满足情况

$f\left(x_{kk,st}\right)-f\left(x_{kk}\right)\&le;{\mathrm{\&xi;\&beta;}}^{mm+nn}\&dtri;f_{kk}^T{d}_{kk}---\left(13\right)$

$f\left(x_{kk,st}\right)-f\left(x_{kk}\right)\&GreaterEqual;(1-\mathrm{\&xi;}){\mathrm{\&beta;}}^{mm+nn}\&dtri;f_{kk}^T{d}_{kk}---\left(14\right)$

若满足，则令步长因子α_kk＝β^mm+nn，f_kk+1＝f(x_kk,st)，停止线性搜索；否则令mm+1→mm，转至步骤Ⅴ继续搜索计算，直至mm超过N₁或者在试探点处满足Armijo不等式(13)和(14)。

步骤7，求下一个迭代点x_kk+1：

x_kk+1＝x_kk+α_kkd_kk(15)

步骤8，计算点x_kk+1处的梯度值计算方法同步骤4。

步骤9，判断非线性优化过程是否收敛，具体是：

Ⅰ判断||x_kk||>ε₆是否成立，若不成立则转步骤Ⅱ，否则再判断||x_kk+1-x_kk||/||x_kk||≤ε₁是否成立，若成立则停止优化迭代，否则转步骤Ⅲ；

Ⅱ判断||x_kk+1-x_kk||≤ε₂是否成立，若成立则停止优化迭代，否则转步骤Ⅲ；

Ⅲ判断|f_kk|>ε₇是否成立，若不成立则转步骤Ⅳ，否则再判断|f_kk-f_kk+1|/|f_kk|≤ε₃是否成立，若成立则停止优化迭代，否则转步骤Ⅴ；

Ⅳ判断|f_kk-f_kk+1|≤ε₄是否成立，若成立则停止优化迭代，否则转步骤Ⅴ；

Ⅴ判断是否成立，若成立则停止优化迭代，否则转步骤10，继续进行优化迭代。

在上述的收敛准则中，ε₁、ε₂、ε₃、ε₄、ε₅为误差限，取为10^-6；ε₆、ε₇分别用于判断||x_k||和|f_k|的数量级，取为10^-4。

步骤10，计算H_kk+1。具体是：

Ⅰ令自变量变化为s_kk＝x_kk+1-x_kk，梯度差为

Ⅱ通过改进的BFGS校正公式(16)和(17)计算逆Hessian阵H_kk+1：

$H_{kk+1}=\left\{\begin{array}{cc}\left(I-\frac{s_{kk}{\stackrel{\&OverBar;}{y}}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}}\right)H_{kk}\left(I-\frac{{\stackrel{\&OverBar;}{y}}_{kk}{s}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}}\right)+\frac{s_{kk}{s}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}},& {\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}>0\\ H_{kk},& {\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}\&le;0\end{array}\right.---\left(16\right)$

其中

${\stackrel{\&OverBar;}{y}}_{kk}=y_{kk}+\frac{12\left(f_{kk}-f_{kk+1}\right)+5\&dtri;f_{kk+1}^T{s}_{kk}+\left(7-{\mathrm{\&alpha;}}_{kk}\right)\&dtri;f_{kk}^T{s}_{kk}}{s_{kk}^T{\mathrm{\&theta;}}_{kk}}{\mathrm{\&theta;}}_{kk}---\left(17\right)$

式中，表示修正的梯度差，θ_kk∈R^m是满足的任意向量，此处取为s_kk。

所述改进的BFGS校正公式，即本实施例中的公式(16)和公式(17)2014年在西北工业大学博士论文《大规模结构并行优化方法及其工程应用研究》中公开。

步骤11，令迭代下标kk→kk+1，转到步骤5继续迭代，直至满足步骤9中的收敛条件。至此，完成了对NACA0012对称翼型的M6机翼的非定常气动力最小状态有理近似的非线性优化算法。

图6给出了当m＝1时广义气动力矩阵中关键模态行两个元素的拟合结果，图中共比较了五种方法，分别为MS-Dr法、MS法、MS-opt法、Roger法和Roger-opt法。其中“MS-Dr法”为本发明所使用的方法，表明本发明采用MS法的拟合形式，且起始于对D矩阵第r行的初始化。“MS-opt法”表示采用现有技术的非线性算法对MS法的气动滞后根进行优化的方法。“Roger-opt法”表示采用现有技术的非线性算法对Roger法的气动滞后根进行优化的方法。图中“Real”表示被拟合项的离散列表值。

为便于对比，“MS-opt法”和“Roger-opt法”中的非线性优化算法均与本发明中的优化算法框架一致。从图中可以看出：

1.在采用相同数量滞后根的情况下，Roger方法的拟合精度高于MS法，但在未进行非线性优化时两种方法都会出现误差较大的情况。而在进行非线性优化后，两者的拟合精度都有明显的改善。这验证了本发明中非线性优化算法的效果。

2.本实施例的精度接近优化后的Roger方法。但因其保留了MS法的形式，因此由其出发得到的状态空间模型能够保持MS法模型阶数低的优点。

表1从拟合误差和计算效率两个方面对比了MS-Dr法与MS法和MS-opt法。其中f_r表示广义气动力矩阵Q(ik)的关键模态行的总拟合误差，如式(18)所示，f表示所有元素的总拟合误差，如式(10)所示。t表示算法的总耗时。所有结果均在一个CPU为2.4GHz，内存4G的计算机上运行得到。

$f_r=\sqrt{\underset{l=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{n}{\&Sigma;}}{W}_{rj}^2\left\{{\&lsqb;F_{ap,rj}\left(k_l\right)-F_{rj}\left(k_l\right)\&rsqb;}^2+{\&lsqb;G_{ap,rj}\left(k_l\right)-G_{rj}\left(k_l\right)\&rsqb;}^2\right\}}---\left(18\right)$

从表1中可以看出MS-Dr法对关键模态行的拟合精度和拟合效率明显优于MS法和MS-opt法，且随着气动滞后根数量的增加，MS-Dr法的优势不断增大。当m＝4时，MS-Dr法的f_r仅为MS-opt法的9％，且其计算时间仅为MS-opt法的0.5％。从f的对比可以看出，MS-Dr法的总体误差略大于MS法和MS-opt法，这是提高关键模态行的拟合精度所不能避免的，但可以通过适当增加滞后根数量来弥补。

表1MS-Dr法与MS法、MS-opt法的拟合效果对比

为深入分析MS-Dr法的拟合效果，以下由其拟合结果出发，构造气动弹性系统的状态空间方程，并用根轨迹法进行颤振分析。具体做法参考文献4“Karpel,M.,TimeDomainAeroservoelasticModelingUsingWeightedUnsteadyAerodynamicForces,JournalofGuidance,Control,andDynamics,Vol.13,No.1,1990,pp.30–37.”。同时，由MS-opt法和Roger-opt法拟合结果得到的状态空间方程也被用来进行颤振分析以进行比较，前者状态空间模型的构建方法同参考文献4，后者参考文献5“Tiffany,S.H.,andAdams,W.M.,Non-LinearProgrammingExtensionstoRationalFunctionApproximationMethodsforUnsteadyAerodynamicForces,NASATP-2776,July1988.”。

对当前模型P-K法计算得到的颤振速度为445.8m/s，颤振频率为13.69Hz。下文将以此为参考结果来检验拟合精度。表2给出了这3种方法时域状态空间模型的颤振分析结果，其中V_f表示颤振速度误差，ω_f表示颤振频率误差。从表2中可以看出，虽然气动滞后根数量的增加会提高模型的整体拟合精度，但由Roger-opt法和MS-opt法得到颤振结果却并没有明显改善，反而可能出现随滞后根增加误差明显增大的情况。这是因为整体模型精度的提高并不能保证关键模态相关项精度的改善，而这对模型的稳定性分析和动态响应分析是至关重要的。但MS-Dr法则能一直保持较高的精度，且随着气动滞后根数量的增加误差有明显减小的趋势。

表2MS-D_r法与MS法、MS-opt法的颤振分析误差对比

为进行气动伺服弹性系统的控制律设计、时域仿真和迭代优化设计，系统的运动方程需要从频域模型转换为一阶时域状态空间模型，而其中的关键技术就是将频域的非定常气动力复延拓到拉氏域，本实施例对当前广泛使用的MS法提出了进一步的改进，在明显提高计算效率和关键模态相关项的拟合精度的同时保证了模型的整体精度，能够准确地计算出模型的颤振特征，在气动伺服弹性系统的分析设计中具有重要意义。

Claims (4)

1.一种非定常气动力最小状态有理近似的非线性优化算法，其特征在于，具体过程是：

步骤1，建立机翼有限元模型和气动力模型，计算机翼在给定减缩频率下的广义气动力系数矩阵，具体是：

Ⅰ建立机翼有限元模型，通过Nastran软件对建立的机翼有限元模型进行模态分析，模态分析中取机翼的前9阶弹性模态；

Ⅱ建立气动力模型，并将得到的模态分析结果导入到Zaero软件中，得到机翼在给定减缩频率下的广义气动力系数矩阵：

Q(ik)＝F(ik)+iG(ik)(1)

其中Q(ik)为给定减缩频率下的广义气动力系数矩阵，k＝ωb/V为减缩频率，ω为机翼振荡的圆频率，b为机翼参考弦长，V为来流速度；F(ik)和G(ik)分别表示广义气动力系数矩阵的实部和虚部；

对Q(ik)拟合时采用MS法的拟合公式，即

$Q_{ap}\left(\stackrel{\&OverBar;}{s}\right)=F_{ap}\left(ik\right)+{\mathrm{iG}}_{ap}\left(ik\right)=A_0+A_1\stackrel{\&OverBar;}{s}+A_2{\stackrel{\&OverBar;}{s}}^2+D{(\stackrel{\&OverBar;}{s}I-R)}^{-1}E\&CenterDot;\stackrel{\&OverBar;}{s}---\left(2\right)$

其中s为拉氏变量；表示拟合得到的广义气动力系数矩阵，F_ap(ik)和G_ap(ik)分别表示的实部和虚部；A₀∈R^n×n、A₁∈R^n×n、A₂∈R^n×n、D∈R^n×m、R∈R^m×m、E∈R^m×n均为待定系数矩阵，n为结构模态数，m为气动滞后根数量；I为m阶单位矩阵；R为由气动滞后根组成的对角矩阵，表示为

R＝diag(x)＝diag([x₁x₂…x_m])(3)

其中x表示气动滞后根组成的向量；

步骤2，给定气动滞后根向量的初值x₀，对MS法的拟合公式进行总体拟合并计算该总体拟合的误差f(x₀)；

给定气动滞后根向量的初值x₀；

根据给定的气动滞后根向量的初值x₀得到由气动滞后根组成的对角矩阵R；对MS法的拟合公式进行总体拟合并计算总体拟合误差f(x₀)，具体过程是：

Ⅰ取结构的第r阶模态为关键模态，并令D矩阵的第r行元素为任意非零常数，其中D矩阵的第r行元素代表结构的第r阶模态对应的广义气动力；此时方程(2)变为

$Q_{ap,r}\left(\stackrel{\&OverBar;}{s}\right)=A_{0,r}+A_{1,r}\stackrel{\&OverBar;}{s}+A_{2,r}{\stackrel{\&OverBar;}{s}}^2+D_r{(\stackrel{\&OverBar;}{s}I-R)}^{-1}E\&CenterDot;\stackrel{\&OverBar;}{s}---\left(4\right)$

其中下标r均表示矩阵的第r行；

Ⅱ约束方程在s＝0处的有理逼近值等于气动力系数矩阵列表值，在s＝ik_f处的实部有理逼近值等于矩阵的列表值，在s＝ik_g处的虚部有理逼近值等于矩阵的列表值，其中k_f和k_g均为指定的减缩频率值；之后从D矩阵第r行出发拟合E矩阵各列元素，所述E矩阵第j列的加权最小二乘求解式如下

$\underset{l=1}{\overset{L}{\&Sigma;}}(P_l^T{W}_{rj}^2{P}_l+{Q_l}^T{W}_{rj}^2{Q}_l)E_j=\underset{l=1}{\overset{L}{\&Sigma;}}\left\{{P_l}^T{W}_{rj}^2{\stackrel{\&OverBar;}{F}}_j\left({\mathrm{ik}}_l\right)+{Q_l}^T{W}_{rj}^2{\stackrel{\&OverBar;}{G}}_j\left({\mathrm{ik}}_l\right)\right\}---\left(5\right)$

其中E_j表示E矩阵第j列，

$\begin{array}{c}{P}_l=k_l^2{D}_r\&lsqb;{(k_l^{2}I+R^2)}^{-1}-{(k_f^{2}I+R^2)}^{-1}\&rsqb;\\ Q_l=k_l{D}_r\&lsqb;{(k_l^{2}I+R^2)}^{-1}-{(k_g^{2}I+R^2)}^{-1}\&rsqb;R\end{array}---\left(6\right)$

$W_{rj}=\frac{1}{\underset{l}{\max }\left\{\right|Q_{rj}\left({\mathrm{ik}}_l\right)|,1\}}---\left(7\right)$

其中L为减缩频率个数；Q_rj(ik_l)表示Q(ik)第r行第j列元素在减缩频率k_l处的值，W_rj表示Q(ik)的加权矩阵W的第r行第j列元素；和分别表示F(ik_l)和G(ik_l)的第j列；

因当前机翼在颤振点处的减缩频率为0.07，故取k_f＝k_g＝0.05以使颤振点附近的拟合精度更高；

Ⅲ求解出E矩阵后，再逐行用加权最小二乘法求解D矩阵的各元素，第r行除外；D矩阵第i行的加权最小二乘求解式如下

$\underset{l=1}{\overset{L}{\&Sigma;}}(P_l^{*T}{W}_i^2{P}_l^{*}+Q_l^{*T}{W}_i^2{Q}_l^T)D_i^T=\underset{l=1}{\overset{L}{\&Sigma;}}\left(P_l^{*T}{W}_i^2{\hat{F}}_i\left({\mathrm{ik}}_l\right)+Q_l^{*T}{W}_i^2{\hat{G}}_i\left({\mathrm{ik}}_l\right)\right)---\left(8\right)$

其中D_i表示D矩阵第i行，

$\begin{array}{c}{P}_l^{*}=k_l^2{E}^T\&lsqb;{(k_l^{2}I+R^2)}^{-1}-{(k_f^{2}I+R^2)}^{-1}\&rsqb;\\ Q_l^{*}=k_l{E}^{T}R\&lsqb;{(k_l^{2}I+R^2)}^{-1}-{(k_g^{2}I+R^2)}^{-1}\&rsqb;\end{array}---\left(9\right)$

W_i ²表示加权矩阵W第i行各元素的平方值构成的对角矩阵，和分别表示F(ik_l)和G(ik_l)的第i行；

Ⅳ计算对MS法的拟合公式进行总体拟合结果的总误差f(x₀)

$f\left(x_0\right)=\sqrt{\underset{l=1}{\overset{L}{\&Sigma;}}\underset{j=1}{\overset{n}{\&Sigma;}}\underset{i=1}{\overset{n}{\&Sigma;}}{W}_{ij}^2\left\{{\&lsqb;F_{ap,ij}\left(k_l\right)-F_{ij}\left(k_l\right)\&rsqb;}^2+{\&lsqb;G_{ap,ij}\left(k_l\right)-G_{ij}\left(k_l\right)\&rsqb;}^2\right\}}---\left(10\right)$

其中i和j分别表示各矩阵的行和列；

步骤3，开始对滞后根向量进行非线性优化；

给定初始Hessian阵的逆阵H₀，为m阶单位阵；令迭代下标kk＝0；

步骤4，计算第一个迭代点x₀处的梯度值其第i项的算式如下

${(\&dtri;f_0)}_i=\frac{f(x_0+\mathrm{\&alpha;}\&CenterDot;e_i)-f(x_0-\mathrm{\&alpha;}\&CenterDot;e_i)}{2\mathrm{\&alpha;}},(i=1,2,...m)---\left(11\right)$

其中e_i为第i个元素为1的m阶单位向量，α为一个充分小的实数，此处取为0.001；f(x₀+α·e_i)和f(x₀-α·e_i)的计算方法同步骤2中f(x₀)的计算方法；

步骤5，确定搜索方向；

通过公式(12)确定搜索方向

$d_{kk}=-H_{kk}\&dtri;f_{kk}---\left(12\right)$

步骤6，沿d_kk线性搜索步长因子α_kk；具体过程是：

Ⅰ给定参数0<ξ<0.5和0<β<1；取线性搜索时计算目标函数值的最大许可次数为N₁；令mm＝0代表本次循环目标函数值的计算次数，令nn＝0代表本次循环自标量的越界次数；

Ⅱ令x_kk,st＝x_kk+β^mm+nnd_kk，x_kk,st代表从x_kk出发得到的试探点；

Ⅲ判断 $s_{\mathrm{kk},\mathrm{st}}^i>{\underset{\&OverBar;}{x}}^i,(i=1,...,m)$ 或 $s_{\mathrm{kk},\mathrm{st}}^i<{\stackrel{\&OverBar;}{x}}^i,(i=1,...,m)$ 是否满足，其中表示x_kk,st的第i项，x ⁱ和分别表示自变量第i项的下界和上界；若都满足则转至步骤Ⅳ；否则令nn+1→nn并转至步骤Ⅱ,其中“→”表示赋值运算；

Ⅳ计算试探点x_kk,st处的目标函数值f(x_kk,st)；

Ⅴ判断mm是否超过N₁；若不超过，则转到步骤Ⅵ；否则取mm为使f(x_kk,st)最小的数，令步长因子α_kk＝β^mm+nn，f_kk+1＝f(x_kk,st)，停止线性搜索；

Ⅵ判断试探点处Armijo不等式的满足情况

$f\left(x_{kk,st}\right)-f\left(x_{kk}\right)\&le;{\mathrm{\&xi;\&beta;}}^{mm+nn}\&dtri;f_{kk}^T{d}_{kk}---\left(13\right)$

$f\left(x_{kk,st}\right)-f\left(x_{kk}\right)\&GreaterEqual;(1-\mathrm{\&xi;}){\mathrm{\&beta;}}^{mm+nn}\&dtri;f_{kk}^T{d}_{kk}---\left(14\right)$

若满足，则令步长因子α_kk＝β^mm+nn，f_kk+1＝f(x_kk,st)，停止线性搜索；否则令mm+1→mm，转至步骤Ⅴ继续搜索计算，直至mm超过N₁或者在试探点处满足Armijo不等式(13)和(14)；

步骤7，求下一个迭代点x_kk+1：

x_kk+1＝x_kk+α_kkd_kk(15)

步骤8，计算点x_kk+1处的梯度值计算方法同步骤4；

步骤9，判断非线性优化过程是否收敛，具体是：

Ⅰ判断||x_kk||>ε₆是否成立，若不成立则转步骤Ⅱ，否则再判断||x_kk+1-x_kk||/||x_kk||≤ε₁是否成立，若成立则停止优化迭代，否则转步骤Ⅲ；

Ⅱ判断||x_kk+1-x_kk||≤ε₂是否成立，若成立则停止优化迭代，否则转步骤Ⅲ；

Ⅲ判断|f_kk|>ε₇是否成立，若不成立则转步骤Ⅳ，否则再判断|f_kk-f_kk+1|/|f_kk|≤ε₃是否成立，若成立则停止优化迭代，否则转步骤Ⅴ；

Ⅳ判断|f_kk-f_kk+1|≤ε₄是否成立，若成立则停止优化迭代，否则转步骤Ⅴ；

Ⅴ判断是否成立，若成立则停止优化迭代，否则转步骤10，继续进行优化迭代；

在上述的收敛准则中，ε₁、ε₂、ε₃、ε₄、ε₅为误差限，取为10^-6；ε₆、ε₇分别用于判断||x_k||和|f_k|的数量级，取为10^-4；

步骤10，计算H_kk+1；具体是：

Ⅰ令自变量变化为s_kk＝x_kk+1-x_kk，梯度差为

Ⅱ通过改进的BFGS校正公式(16)和(17)计算逆Hessian阵H_kk+1：

$H_{kk+1}=\left\{\begin{array}{cc}\left(I-\frac{s_{kk}{\stackrel{\&OverBar;}{y}}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}}\right)H_{kk}\left(I-\frac{{\stackrel{\&OverBar;}{y}}_{kk}{s}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}}\right)+\frac{s_{kk}{s}_{kk}^T}{{\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}},& {\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}>0\\ H_{kk},& {\stackrel{\&OverBar;}{y}}_{kk}^T{s}_{kk}\&le;0\end{array}\right.---\left(16\right)$

其中

${\stackrel{\&OverBar;}{y}}_{kk}=y_{kk}+\frac{12(f_{kk}-f_{kk+1})+5\&dtri;f_{kk+1}^T{s}_{kk}+(7-{\mathrm{\&alpha;}}_{kk})\&dtri;f_{kk}^T{s}_{kk}}{s_{kk}^T{\mathrm{\&theta;}}_{kk}}{\mathrm{\&theta;}}_{kk}---\left(17\right)$

式中，表示修正的梯度差，θ_kk∈R^m是满足的任意向量，此处取为s_kk；

步骤11，令迭代下标kk→kk+1，转到步骤5继续迭代，直至满足步骤9中的收敛条件；至此，完成了对NACA0012对称翼型的M6机翼的非定常气动力最小状态有理近似的非线性优化算法。

2.如权利要求1所述一种非定常气动力最小状态有理近似的非线性优化算法，其特征在于，所述建立的机翼有限元模型中，机翼下翼面的弦向剖面共有8个网格节点分别位于距前缘0.0％、5.0％、10.0％、27.0％、45.5％、63.5％、82.0％和100％处；沿机翼展向各网格节点数量与翼肋的数量相同，并与各翼肋的展向位置一致。

3.如权利要求1所述一种非定常气动力最小状态有理近似的非线性优化算法，其特征在于，所述建立的气动力模型的弦向被均分为15等份，展向被均分为24等份。

4.如权利要求1所述一种非定常气动力最小状态有理近似的非线性优化算法，其特征在于，在对得到的广义气动力系数矩阵进行拟合前，需先给定气动滞后根向量的初值，初始化方法为：当m＝1时，x₀＝[-0.3]；当m＝2时，x₀＝[-0.3-0.5]；当m＝3时，x₀＝[-0.3-0.5-0.7]；当m＝4时，x₀＝[-0.3-0.5-0.7-0.9]。