[发明专利]非定常气动力最小状态有理近似的非线性优化算法

权利要求书

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法的拟合公式，即

Qap(s‾)=Fap(ik)+iGap(ik)=A0+A1s‾+A2s‾2+D(s‾I-R)-1E·s‾---(2)]]>

其中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)变为

Qap,r(s‾)=A0,r+A1,rs‾+A2,rs‾2+Dr(s‾I-R)-1E·s‾---(4)]]>

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

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

Σl=1L(PlTWrj2Pl+QlTWrj2Ql)Ej=Σl=1L{PlTWrj2F‾j(ikl)+QlTWrj2G‾j(ikl)}---(5)]]>

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

Pl=kl2Dr[(kl2I+R2)-1-(kf2I+R2)-1]Ql=klDr[(kl2I+R2)-1-(kg2I+R2)-1]R---(6)]]> 1

Wrj=1maxl{|Qrj(ikl)|,1}---(7)]]>

其中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行的加权最小二乘求解式如下

Σl=1L(Pl*TWi2Pl*+Ql*TWi2Ql*)DiT=Σl=1L(Pl*TWi2F^i(ikl)+Ql*TWi2G^i(ikl))---(8)]]>

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

Pl*=kl2ET[(kl2I+R2)-1-(kf2I+R2)-1]Ql*=klETR[(kl2I+R2)-1-(kg2I+R2)-1]---(9)]]>

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

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

f(x0)=Σl=1LΣj=1nΣi=1nWij2{[Fap,ij(kl)-Fij(kl)]2+[Gap,ij(kl)-Gij(kl)]2}---(10)]]>

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

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

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

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

(▿f0)i=f(x0+α·ei)-f(x0-α·ei)2α,(i=1,2,...m)---(11)]]>

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

步骤5，确定搜索方向；

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

dkk=-Hkk▿fkk---(12)]]>

步骤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(xkk,st)-f(xkk)≤ξβmm+nn▿fkkTdkk---(13)]]>

f(xkk,st)-f(xkk)≥(1-ξ)βmm+nn▿fkkTdkk---(14)]]>

若满足，则令步长因子α_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：

Hkk+1=(I-skky‾kkTy‾kkTskk)Hkk(I-y‾kkskkTy‾kkTskk)+skkskkTy‾kkTskk,y‾kkTskk>0Hkk,y‾kkTskk≤0---(16)]]>

其中

y‾kk=ykk+12(fkk-fkk+1)+5▿fkk+1Tskk+(7-αkk)▿fkkTskkskkTθkkθkk---(17)]]>

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

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