CN103116706A - Configured control optimization method for high-speed aircrafts based on pneumatic nonlinearity and coupling - Google Patents

Abstract

The invention discloses a configured control optimization method for high-speed aircrafts based on pneumatic nonlinearity and coupling and aims to solve the technical problem that the existing configured control optimization method for high-speed aircrafts is poor. According to the technical scheme, the method includes: firstly, defining pneumatic nonlinearity evaluation indicators and pneumatic coupling evaluation indicators; secondly determining configured control optimization target functions; and thirdly establishing a configured control optimization model. Optimization is targeted at the pneumatic nonlinearity evaluation indicators and pneumatic coupling evaluation indicators of aircrafts, and a multi-target configured control optimization algorithm in coordination of ideal point method and simulated annealing algorithm is adopted. Pneumatic nonlinearity and coupling features of the aircrafts are decreased by continuously changing the appearances of the aircrafts, so that a linear model is more similar to a nonlinear model, a decoupling model and a coupling model and the pneumatic nonlinearity and coupling features of the aircrafts are optimal under certain limitations.

Description

Based on pneumatic non-linear with the coupling high-speed aircraft with the control optimization method

Technical field

The present invention relates to a kind of high-speed aircraft with the control optimization method, particularly relate to a kind of based on pneumatic non-linear with the coupling high-speed aircraft with the control optimization method.

Background technology

The optimization of high-speed aircraft aerodynamic configuration combines the aerodynamic performance analysis of design object with best practice, by the profile of continuous change design object, the aeroperformance that makes aircraft is issued to optimum satisfying certain constraint condition.

The profile optimization of hypersonic aircraft is continued to use the Multidisciplinary Optimization method of conventional aircraft, use engineering method as document 1 " Advancements in Multidisciplinary Design Optimization Applied to Hypersonic Vehicles to Achieve Closure (AIAA-2008-2591) " based on Parametric CAD model and pneumatic subject, completed lifting body, the hypersonic multidisciplinary multiple-objection optimization that repeats vehicle (RLV) of Waverider.Document 2 " pneumatic design of quafric curve interface bomb body and optimization (Tang Wei; Zhang Yong; Li Weiji. aerospace journal; 2004; 25 (4): 429-433) " set up the bomb body profile with the quafric curve method, in adopting, the volt Newtonian theory is estimated hypersonic aerodynamic characteristic, has realized having the reentry vehicle layout optimization in quafric curve cross section.Document 3 " based on the Waverider optimal design of nerual network technique (Zhang Fengtao; Cui Kai, Yang Guowei. mechanics journal, 2009; 41 (3): 418-423) " neural network response surface technology is combined with the CFD aerodynamic analysis, realized the aerodynamic configuration optimization of rider layout.

The aerodynamic configuration optimization of the disclosed high-speed aircraft of document is intended to increase lift-drag ratio, the slope of lift curve of aircraft.Aircraft aerodynamic configuration optimization research is still common with aerofoil profile, pays close attention to relatively less to the optimal design of fuselage bomb body.In addition, the optimization of aircraft aerodynamic configuration is normal adopts partial approach based on gradient as the optimal design algorithm, and the method easily is absorbed in local optimum and causes the of overall importance relatively poor of optimization solution.Traditional disposal route is that rule of thumb formula and decision-making technique are compromised to performance, and the method efficient is low, low precision.Development along with computing technique, numerical simulation combines with multidisciplinary optimization (MDO) technology becomes the Main Means of current design, as can process in conjunction with empirical method and Hi-Fi numerical modeling method the initial question of pneumatic/structure Coupling in MDO.

Summary of the invention

In order to overcome existing high-speed aircraft with the poor deficiency of control optimization method effect of optimization, the invention provides a kind of high-speed aircraft based on pneumatic non-linear and coupling with the control optimization method.The method as optimization aim, adopts multiple goal that ideal point method and simulated annealing cooperates mutually with controlling optimized algorithm the pneumatic non-linear evaluation index of aircraft and pneumatic degree of coupling evaluation index.By change of flight device profile constantly to reduce the pneumatic non-linear and coupled characteristic of aircraft, improve the similarity of inearized model and nonlinear model, Decoupled Model and coupling model, the pneumatic non-linear and coupled characteristic that makes aircraft is issued to optimum satisfying certain constraint condition.

The technical solution adopted for the present invention to solve the technical problems is: a kind of based on pneumatic non-linear with the coupling high-speed aircraft with the control optimization method, be characterized in comprising the following steps:

Step 1, be the pneumatic nonlinearity of pitch channel with pneumatic non-linear evaluation index definition

, namely

$\mathrm{DON}{L}_{m_z\left(\mathrm{\α}\right)}=\frac{|C_m\left(\mathrm{\α}\right)-C_{m\_L}\left(\mathrm{\α}\right)|}{|C_m\left(\mathrm{\α}\right)-C_{m\_L}\left(\mathrm{\α}\right)|+\left|C_{m\_L}\left(\mathrm{\α}\right)\right|}---\left(1\right)$

$=\frac{|m_z\left(\mathrm{\α}\right)-m_z\left(0\right)-{m_z}^{\′}\left(0\right)\mathrm{\α}|}{|m_z\left(\mathrm{\α}\right)-m_z\left(0\right)-{m_z}^{\′}\left(0\right)\mathrm{\α}|+|m_z\left(0\right)+{m_z}^{\′}\left(0\right)\mathrm{\α}|}$

In formula, C _m(α) be the stable pitching moment coefficient of aircraft; C _{m_L}(α) be C _mLinear segment (α).

1) C _m(α) computing formula.

C _m(α)=C _{Bomb body}+ K _{α α}C _Aerofoil=C _{M_zhi1}+ C _{M_zhi2}+ C _{M_zhu}+ K _{α α}C _{m_w}(2)

Wherein,

$C_{m\_\mathrm{<figure-callout\; id="1"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}1}=\frac{4\sin \mathrm{\&alpha;}\cos \mathrm{\&alpha;si}{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">n</figure-callout>}}^2{\mathrm{\&theta;}}_1}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1\tan {\mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2\tan {\mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}(\frac{1}{2}{x}_{\mathrm{<figure-callout\; id="1"\; label="c\; L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}{L}_{}^{}{}_1^2-\frac{1}{3}{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1^3)$

$C_{m\_\mathrm{<figure-callout\; id="2"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}2}=\frac{4\sin \mathrm{\&alpha;}\cos \mathrm{\&alpha;si}n{\mathrm{\&theta;}}_2\cos {\mathrm{\&theta;}}_2}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1\tan {\mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2\tan {\mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}[L_1{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_1(x_c-L_1-\frac{1}{2}{L}_2)+{{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}^2{\tan \mathrm{\&theta;}}_2(\frac{1}{2}(x_c-L_1)-\frac{1}{3}{L}_2)]$

$C_{m\_\mathrm{zhu}}=\frac{8{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&alpha;}{L}_3[x_c-L_3/2-({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1+L_2)}{3\mathrm{\&pi;}({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1{\tan \mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_2)L_{\mathrm{ref}}}$

In formula: L ₁Be pointed cone length, L ₂Be truncated cone length, L ₃Be shell of column length, θ ₁Be pointed cone cone angle, θ ₂Be the truncated cone cone angle, α is the angle of attack, x _cBe centroid position, L _refBe reference length.

When the flat shape of the wing is two swept-back wing, its stable pitching moment coefficient C _{m_w}For:

$C_{m\_w}=\frac{(\mathrm{\&gamma;}+1){\mathrm{\&alpha;}}^2}{S_{\mathrm{<figure-callout\; id="1"\; label="ref\; L\; ref"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">ref</figure-callout>}}{L}_{\mathrm{ref}}{}_{}}\sqrt{1+{\left[\frac{4}{(\mathrm{\&gamma;}+1)M_{\&infin;}\mathrm{\&alpha;}}\right]}^2}\left\{{\cot \mathrm{\&lambda;}}_I\right[{b_{0I}}^2(x_c-x_s)-\frac{2}{3}{b_{0I}}^3]-$

$\frac{2}{3}\cot {\mathrm{\&lambda;}}_{\mathrm{II}}{({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1)}^3+[\cot {\mathrm{\&lambda;}}_{\mathrm{II}}(x_c-x_s)-b_{0I}({\cot \mathrm{\&lambda;}}_I+{\cot \mathrm{\&lambda;}}_{\mathrm{II}})]{({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1)}^2$

$+{2b}_{0I}{\cot \mathrm{\&lambda;}}_I(b_{0\mathrm{II}}-b_1)(x_c-x_s-b_{0I})$

$+[({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1){\cot \mathrm{\&lambda;}}_{\mathrm{II}}+b_{0I}{\cot \mathrm{\&lambda;}}_I][2(x_c-b_{0I}-{\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-x_s)b_1+{{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_1}^2]\}$

In formula: γ is coefficient of heat insulation, S _refBe area of reference, L _refBe reference length, λ _IBe inner wing leading edge sweep, λ _IIBe outer wing leading edge sweep, b _0IBe interior wing root chord length, b _0IIBe outer wing root chord length, b ₁Be aerofoil tip string, x _sBe the exposed wing installation site distance apart from the bomb body summit.

When the flat shape of the wing is single swept-back wing, its stable pitching moment coefficient C _{m_w}For:

$C_{m\_w}=\frac{M_z}{q_{\&infin;}{S}_{\mathrm{ref}}{L}_{\mathrm{ref}}}=\frac{(\mathrm{\&gamma;}+1){\mathrm{\&alpha;}}^2}{S_{\mathrm{ref}}{L}_{\mathrm{ref}}}{\cot \mathrm{\&lambda;}}_{\mathrm{<figure-callout\; id="1"\; label="I"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">I</figure-callout>}}\sqrt{1+{\left[\frac{4}{(\mathrm{\&gamma;}+1)M_{\&infin;}\mathrm{\&alpha;}}\right]}^2}[{(b_{0I}-b_1)}^2(x_c-x_s)$

$-\frac{2}{3}{(b_{0I}-b_1)}^3+{2b}_1(b_{0I}-b_1)(x_c-x_s-b_{0I}+b_1)-(b_{0I}-b_1){b_1}^2]$

In formula: γ is coefficient of heat insulation, S _refBe area of reference, L _refBe reference length, λ _IBe inner wing leading edge sweep, b _0IBe interior wing root chord length, b ₁Be aerofoil tip string, x _sBe the exposed wing installation site distance apart from the bomb body summit, K _{α α}It is wing body interference coefficient.

2) C _{m_L}(α) computing formula.

C _{m_L}(α)＝C _{m_L_zhi1}+C _{m_L_zhi2}+C _{m_wb_L}（3）

Wherein,

$C_{m\_\mathrm{<figure-callout\; id="1"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}1\_L}=\frac{4\mathrm{\&alpha;}{\sin }^2{\mathrm{\&theta;}}_1}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1{\tan \mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}(\frac{1}{2}{x}_{\mathrm{<figure-callout\; id="1"\; label="c\; L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}{L}_{}^{}{}_1^2-\frac{1}{3}{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1^3)$

$C_{m\_\mathrm{<figure-callout\; id="2"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}2\_L}=\frac{4\mathrm{\&alpha;}{\sin \mathrm{\&theta;}}_2{\cos \mathrm{\&theta;}}_2}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1{\tan \mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}[\frac{1}{2}{L}_1{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_1(2x_c-{2L}_1-L_2)$

$+{{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}^2{\tan \mathrm{\&theta;}}_2(\frac{1}{2}(x_c-L_1)-\frac{1}{3}{L}_2)$

$C_{m\_\mathrm{wb}\_L}=-0.035\mathrm{\&alpha;\&xi;}\frac{(1-{\mathrm{\&eta;}}_t^2)}{L_{\mathrm{ref}}}(x_c-L_{\mathrm{sh}}-R_{\max }{\mathrm{\&eta;}}_t{\cot \mathrm{\&theta;}}_3+\frac{2}{3}\frac{R_{\max }{\cot \mathrm{\&theta;}}_3(1-{\mathrm{\&eta;}}_t^4)}{1-{\mathrm{\&eta;}}_t^3})$

In formula: L ₁Be pointed cone length, L ₂Be truncated cone length, θ ₁Be pointed cone cone angle, θ ₂Be truncated cone cone angle, θ ₃Be the afterbody angle of throat, α is the angle of attack, x _cBe centroid position, L _refBe reference length, η _tBe afterbody shrinkage ratio, R _maxBe the column part radius.

The degree of coupling of the rudder kick of step 2, the pneumatic degree of coupling evaluation index of definition to roll channel

$\frac{m_x^{{\mathrm{\&delta;}}_y}}{m_x^{{\mathrm{\&delta;}}_x}}=p\frac{{2k}_T{S}_{\mathrm{cw}}{y}_r{\left(K_{\mathrm{\&delta;}0}\right)}_{\mathrm{cwr}}{S}_{\mathrm{cwr}}{\left(K_{\mathrm{\&delta;}0}\right)}_W}{\mathrm{SL}{\left(K_{\mathrm{\&delta;}0}\right)}_{\mathrm{Wr}}{S}_{\mathrm{Wr}}\mathrm{\&sigma;}\frac{({\mathrm{\&eta;}}_k+1){(1-\stackrel{\&OverBar;}{D})}^2}{{\mathrm{\&eta;}}_k+1-2\stackrel{\&OverBar;}{D}}[\stackrel{\&OverBar;}{D}+(1-\stackrel{\&OverBar;}{D})f]}---\left(4\right)$

In formula, when aircraft is single vertical fin layout, p=1; When aircraft was twin-finned layout, p was the function about (α, Ma, h).k _TBe air-flow retardation factor, S _cwBe the vertical fin area, S is aircraft feature area, and L is the aircraft characteristic length, y _rBe the distance of the yaw rudder center of area to the bomb body longitudinal axis, S _cwrBe rudder area, S _WrBe vertical fin area, (K _{δ 0}) _cwrBe the interference coefficient between yaw rudder and wing body, (K _{δ 0}) _WBe wing body interference coefficient, (K _{δ 0}) _WrBe the interference coefficient between vertical fin and wing body, η _kBe contraction coefficient,

Be footpath exhibition ratio, f exposes the wing root chord section to the distance of pressing heart other shore with half length, and σ is that the bomb body relative diameter changes the correction factor that rolling moment is exerted an influence and introduces.

Define the rudder kick of pneumatic degree of coupling evaluation index to the degree of coupling of pitch channel

$\frac{m_z^{{\mathrm{\&delta;}}_y}}{m_z^{{\mathrm{\&delta;}}_z}}={\mathrm{qa}}_0(1+a_{1}h)(1+a_2\mathrm{\&alpha;})(1+a_3\mathrm{Ma})(1+a_4{\mathrm{\&eta;}}_{\mathrm{kew}})(1+a_5{\mathrm{\&lambda;}}_{\mathrm{kew}})(1+a_6\cos {\mathrm{\&chi;}}_{\mathrm{cw}})---\left(5\right)$

$(1+a_7{\mathrm{\&lambda;}}_{\mathrm{kcwr}})(1+a_8{y}_{\mathrm{cwr}})(1+a_9{L}_{\mathrm{cw}})(1+a_{10}{\stackrel{\&OverBar;}{D}}_{\mathrm{cw}})(1+a_{11}\frac{S_{\mathrm{cwr}}}{S})(1+a_{12}\frac{S_{\mathrm{cw}}}{S})$

In formula, when aircraft is single vertical fin layout, q=1; When aircraft was twin-finned layout, q was the function about (α, Ma, h).S _cw, η _kcw, λ _kcw, χ _cwBe respectively the area of vertical tail, expose the vertical fin contraction coefficient, expose vertical fin aspect ratio and angle of sweep; S _cwr, λ _Kcwr,

L _cwThe area that represents respectively yaw rudder, the yaw rudder aspect ratio, the footpath exhibition of aircraft vertical fin is compared and is gone out the outer crab rooting string mid point in position to the distance of Vehicle nose.

Define the Jenkel rudder deflection of pneumatic degree of coupling evaluation index to the degree of coupling of jaw channel

$\frac{m_y^{{\mathrm{\&delta;}}_x}}{m_y^{{\mathrm{\&delta;}}_y}}={\mathrm{mb}}_0(1+b_{1}h)(1+{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_2\mathrm{\&alpha;})(1+b_3\mathrm{Ma})(1+b_4{\mathrm{\&eta;}}_{\mathrm{kw}})(1+b_5{\mathrm{\&lambda;}}_{\mathrm{kw}})(1+b_6\cos {\mathrm{\&chi;}}_w)---\left(6\right)$

$(1+b_7{\mathrm{\&lambda;}}_{\mathrm{kwr}})(1+b_8{y}_{\mathrm{wr}})(1+b_9{L}_w)(1+b_{10}{\stackrel{\&OverBar;}{D}}_{\mathrm{cw}})(1+b_{11}\frac{S_{\mathrm{wr}}}{S})(1+b_{12}\frac{S_w}{S})$

In formula, when aircraft is single sweepback layout, m=1; When aircraft was two sweepback layout, m was the function about (α, Ma, h).S _wThe aerofoil area, η _kwExpose the vertical fin contraction coefficient, λ _kwBe the aspect ratio of the wing, χ _wBe vertical fin angle of sweep, S _wrBe Jenkel rudder area, λ _kwrBe the Jenkel rudder aspect ratio,

Be aircraft missile wing footpath exhibition ratio, L _wFor exposing the wing root chord mid point to the distance of Vehicle nose, y _wrBe the distance of the Jenkel rudder center of area to the bomb body longitudinal axis.

Define the Jenkel rudder deflection of pneumatic degree of coupling evaluation index to the degree of coupling of pitch channel

$\frac{m_z^{{\mathrm{\&delta;}}_x}}{m_z^{{\mathrm{\&delta;}}_z}}={\mathrm{nc}}_0(1+c_{1}h)(1+{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}_2\mathrm{\&alpha;})(1+c_3\mathrm{Ma})(1+c_4{\mathrm{\&eta;}}_{\mathrm{kw}})(1+c_5{\mathrm{\&lambda;}}_{\mathrm{kw}})(1+c_6\cos {\mathrm{\&chi;}}_w)---\left(7\right)$

$(1+c_7{\mathrm{\&lambda;}}_{\mathrm{kwr}})(1+c_8{y}_{\mathrm{wr}})(1+c_9{L}_w)(1+c_{10}{\stackrel{\&OverBar;}{D}}_w)(1+c_{11}\frac{S_{\mathrm{wr}}}{S})(1+c_{12}\frac{S_w}{S})$

In formula, when aircraft is single sweepback layout, n=1; When aircraft was two sweepback layout, n was the function about (α, Ma, h).S _wThe aerofoil area, η _kwExpose the vertical fin contraction coefficient, λ _kwBe the aspect ratio of the wing, χ _wBe vertical fin angle of sweep, S _wrBe Jenkel rudder area, λ _kwrBe the Jenkel rudder aspect ratio, Be aircraft missile wing footpath exhibition ratio, L _wFor exposing the wing root chord mid point to the distance of Vehicle nose, y _wrBe the distance of the Jenkel rudder center of area to the bomb body longitudinal axis.

Coefficient a in formula (5), formula (6) and formula (7) ₀～a ₁₂, b ₀～b ₁₂, c ₀～c ₁₂The match value that is based on a large amount of CFD computational datas and utilizes the nonlinear multivariable least-square fitting approach to calculate.

$\left\{\begin{array}{c}{a}_0=0.08263,a_1=6.138e-05,a_2=-0.08484,a_3=-0.03714,a_4=0.02908\\ a_5=0.04677,a_6=0.05188,a_7=0.01160,a_8=0.05631,a_9=0.00876\\ a_{10}=0.06555,a_{11}=0.59681,a_{12}=0.15504\end{array}\right.$

$\left\{\begin{array}{c}{b}_0=1.01895,{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_1=4.422e-05,{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_2=6.5758,b_3=-0.03879,b_4=-0.03635\\ b_5=-0.24995,b_6=-1.19032,b_7=-0.05226,b_8=-0.25926,b_9=-0.10399\\ b_{10}=-0.74664,b_{11}=-2.13691,b_{12}=-0.06801\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="c"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">c</figure-callout>}}_0=-0.0269,{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}_1=-1.2473e-05,{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}_2=-0.6602,c_3=-0.0121,c_4=0.0124\\ c_5=0.0737,c_6=0.4357,c_7=0.0178,c_8=0.0772,c_9=0.0343\\ c_{10}=0.2224,c_{11}=0.6530,c_{12}=0.0233\end{array}\right.$

Step 3, aircraft on overall trajectory arbitrarily unique point i place all have one group of pneumatic non-linear and degree of coupling evaluation index:

$\{{\left({\mathrm{DONL}}_{m_z\left(\mathrm{\&alpha;}\right)}\right)}_i,{(m_x^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z})}_i,{(m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z})}_i,{(m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y})}_i,{(m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z})}_i\}$

Wherein, i=1,2...N, N are the sum of unique point on overall trajectory.

Therefore evaluation index requirement owing to carrying out need to satisfying on overall trajectory with the result that control is optimized non-linear and the degree of coupling at unique point arbitrarily place at overall trajectory multi-characteristic points place needs construct with minor function:

$f_1\left(X\right)=\max \left\{{\left|{\mathrm{DONL}}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\},$ i＝1,2...N

$f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\},$ i＝1,2...N

$f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\},$ i＝1,2...N

$f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\},$ i＝1,2...N

$f_5\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\},$ i＝1,2...N

In formula, X is the vector (θ that contour structures parameter and the flight status parameter by aircraft consists of ₁, θ ₂... h, Ma, α) ^T

The purpose of optimizing with control according to high-speed aircraft and the pneumatic nonlinearity of aircraft, the definition of degree of coupling evaluation index, the objective function that high-speed aircraft overall trajectory multi-characteristic points is optimized with control is:

${\min f}_1\left(X\right)=\max \left\{{\left|\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\}$

$\min f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\}$

$\min f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}$

$\min f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\}$

$\min f_5\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}$

Wherein, i=1,2...N, N are the sum of unique point on overall trajectory.

Step 4, set up high-speed aircraft overall trajectory multi-characteristic points with the control Optimized model:

minf ₁(X)

minf ₂(X)

；（8）

minf ₅(X)

s.t.g _i(X)≤0,i＝1,2,...m

In formula, g _i(X)≤0, i=1,2 ... m is the constraint condition with the control Optimized model.

By ideal point method, multi-objective optimization question is converted to single-object problem, is about to a plurality of objective functions and synthesizes the simple target function.Can change the expression formula of rear single-goal function h (X) according to the computing formula of ideal point method:

$h\left(X\right)={\left[\underset{j=1}{\overset{5}{\mathrm{\&Sigma;}}}\frac{{(f_j\left(x\right)-f_j^0)}^2}{f_j^0}\right]}^{1/2}---\left(9\right)$

Wherein,

Be function f _j(X) ideal value.

Utilize ideal point method to obtain the as follows with the control Optimized model of single-goal function:

min?h(X)（10）

s.t.g _i(X)≤0,i＝1,2,...m

Utilize simulated annealing single-goal function with the feasible zone of control Optimized model in search satisfy the optimization solution of pneumatic nonlinearity and degree of coupling evaluation index requirement.In simulated annealing utilization simulation thermodynamics, the temperature-fall period of classical particle system is found the solution the extreme value of planning problem, and its basic step is:

1. given cooling program parameter and iteration initial solution X thereof ₀And h (X ₀);

2. parametric t=t _kThe time make L _kInferior heuristic search;

According to current solution X _kProduce a random vector Z _kObtain X _kThe new exploration point X ' of neighborhood _k:

X′ _k＝X _k+Z _k

3. produce one at (0,1) upper equally distributed random number η, calculate at given current iteration point X _kAnd temperature t _kThe lower transition probability P corresponding with the Metropolis acceptance criterion has:

$P=\left\{\begin{array}{cc}1,& h\left(X_k^{\&prime;}\right)\&le;h\left(X_k\right)\\ e^{\frac{h\left(X_k\right)-h\left(X_k^{\&prime;}\right)}{t_k}},& h\left(X_k^{\&prime;}\right)>h\left(X_k\right)\end{array}\right.$

If η≤P accepts new explanation X _k=X ' _k, h (X _k)=h (X ' _k), otherwise X _kConstant.

4. sound out point search less than L _kInferior, turn step 2., otherwise turn step 5.;

5. do you satisfy stopping criterion for iteration?

Be, stop, current solution is Approximate Global Optimal Solution; No, turn step 6..

6. according to given temperature damping's function, produce new temperature control parameter t _k+1And the length L of Markov chain _k+1

7. repeating step 2.-6., until find optimum solution.

The invention has the beneficial effects as follows: due to the pneumatic non-linear evaluation index of aircraft and pneumatic degree of coupling evaluation index as optimization aim, adopt multiple goal that ideal point method and simulated annealing cooperates mutually with controlling optimized algorithm.By change of flight device profile constantly to reduce the pneumatic non-linear and coupled characteristic of aircraft, improve the similarity of inearized model and nonlinear model, Decoupled Model and coupling model, the pneumatic non-linear and coupled characteristic that makes aircraft has been issued to optimum satisfying certain constraint condition.

Below in conjunction with drawings and Examples, the present invention is elaborated.

Description of drawings

Fig. 1 is the related aircraft profile vertical view of the inventive method.

Fig. 2 is the related aircraft profile side view of the inventive method.

Fig. 3 is the process flow diagram of the inventive method.

Embodiment

The present invention is based on pneumatic non-linear high-speed aircraft with coupling as follows with control optimization method concrete steps:

With reference to Fig. 1～3.The present invention is directed to wing body assembly aircraft and carry out aerodynamic characteristic with control optimization.The aircraft bomb body is made of pointed cone section, truncated cone section, cylindrical section and contraction afterbody; The aircraft missile wing is divided into single swept-back wing, two swept-wing layout; The aircraft vertical fin is divided into single vertical fin, twin-finned layout.In addition, elevating rudder (Jenkel rudder) and yaw rudder are installed on aircraft missile wing, vertical fin.

One, the definition of pneumatic non-linear evaluation index.

Pneumatic non-linear evaluation index is the pneumatic nonlinearity of pitch channel

, it is defined as:

$\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}=\frac{|C_m\left(\mathrm{\&alpha;}\right)-C_{m\_L}\left(\mathrm{\&alpha;}\right)|}{|C_m\left(\mathrm{\&alpha;}\right)-C_{m\_L}\left(\mathrm{\&alpha;}\right)|+\left|C_{m\_L}\left(\mathrm{\&alpha;}\right)\right|}---\left(1\right)$

$=\frac{|m_z\left(\mathrm{\&alpha;}\right)-m_z\left(0\right)-{m_z}^{\&prime;}\left(0\right)\mathrm{\&alpha;}|}{|m_z\left(\mathrm{\&alpha;}\right)-m_z\left(0\right)-{m_z}^{\&prime;}\left(0\right)\mathrm{\&alpha;}|+|m_z\left(0\right)+{m_z}^{\&prime;}\left(0\right)\mathrm{\&alpha;}|}$

Wherein, C _m(α) be the stable pitching moment coefficient of aircraft; C _{m_L}(α) be C _mLinear segment (α).

1) C _m(α) computing formula.

C _m(α)=C _{Bomb body}+ K _{α α}C _Aerofoil=C _{M_zhi1}+ C _{M_zhi2}+ C _{M_zhu}+ K _{α α}C _{m_w}(2)

Wherein,

$C_{m\_\mathrm{<figure-callout\; id="1"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}1}=\frac{4\sin \mathrm{\&alpha;}\cos \mathrm{\&alpha;si}{n}^2{\mathrm{\&theta;}}_1}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1\tan {\mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2\tan {\mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}(\frac{1}{2}{x}_{\mathrm{<figure-callout\; id="1"\; label="c\; L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}{L}_{}^{}{}_1^2-\frac{1}{3}{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1^3)$

$C_{m\_\mathrm{<figure-callout\; id="2"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}2}=\frac{4\sin \mathrm{\&alpha;}\cos \mathrm{\&alpha;si}n{\mathrm{\&theta;}}_2\cos {\mathrm{\&theta;}}_2}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1\tan {\mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2\tan {\mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}[L_1{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_1(x_c-L_1-\frac{1}{2}{L}_2)+{{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}^2{\tan \mathrm{\&theta;}}_2(\frac{1}{2}(x_c-L_1)-\frac{1}{3}{L}_2)]$

$C_{m\_\mathrm{zhu}}=\frac{8{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2\mathrm{\&alpha;}{L}_3[x_c-L_3/2-({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1+L_2)}{3\mathrm{\&pi;}({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1{\tan \mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_2)L_{\mathrm{ref}}}$

In formula: L ₁Be pointed cone length, L ₂Be truncated cone length, L ₃Be shell of column length, θ ₁Be pointed cone cone angle, θ ₂Be the truncated cone cone angle, α is the angle of attack, x _cBe centroid position, L _refBe reference length.

When the flat shape of the wing is two swept-back wing, its stable pitching moment coefficient C _{m_w}For:

$C_{m\_w}=\frac{(\mathrm{\&gamma;}+1){\mathrm{\&alpha;}}^2}{S_{\mathrm{<figure-callout\; id="1"\; label="ref\; L\; ref"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">ref</figure-callout>}}{L}_{\mathrm{ref}}{}_{}}\sqrt{1+{\left[\frac{4}{(\mathrm{\&gamma;}+1)M_{\&infin;}\mathrm{\&alpha;}}\right]}^2}\left\{{\cot \mathrm{\&lambda;}}_I\right[{b_{0I}}^2(x_c-x_s)-\frac{2}{3}{b_{0I}}^3]-$

$\frac{2}{3}\cot {\mathrm{\&lambda;}}_{\mathrm{II}}{({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1)}^3+[\cot {\mathrm{\&lambda;}}_{\mathrm{II}}(x_c-x_s)-b_{0I}({\cot \mathrm{\&lambda;}}_I+{\cot \mathrm{\&lambda;}}_{\mathrm{II}})]{({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1)}^2$

$+{2b}_{0I}{\cot \mathrm{\&lambda;}}_I({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1)(x_c-x_s-b_{0I})$

$+[({\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-b_1){\cot \mathrm{\&lambda;}}_{\mathrm{II}}+b_{0I}{\cot \mathrm{\&lambda;}}_I][2(x_c-b_{0I}-{\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}-x_s){\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_1+{{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_1}^2]\}$

In formula: γ is coefficient of heat insulation, S _refBe area of reference, L _refBe reference length, λ _IBe inner wing leading edge sweep, λ _IIBe outer wing leading edge sweep, b _0IBe interior wing root chord length, b _0IIBe outer wing root chord length, b ₁Be aerofoil tip string, x _sBe the exposed wing installation site distance apart from the bomb body summit.

When the flat shape of the wing is single swept-back wing, its stable pitching moment coefficient C _{m_w}For:

$C_{m\_w}=\frac{M_z}{q_{\&infin;}{S}_{\mathrm{ref}}{L}_{\mathrm{ref}}}=\frac{(\mathrm{\&gamma;}+1){\mathrm{\&alpha;}}^2}{S_{\mathrm{ref}}{L}_{\mathrm{ref}}}{\cot \mathrm{\&lambda;}}_I\sqrt{1+{\left[\frac{4}{(\mathrm{\&gamma;}+1)M_{\&infin;}\mathrm{\&alpha;}}\right]}^2}[{(b_{0I}-b_1)}^2(x_c-x_s)$

$-\frac{2}{3}{(b_{0I}-b_1)}^3+{2b}_1(b_{0I}-b_1)(x_c-x_s-b_{0I}+b_1)-(b_{0I}-b_1){{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_1}^2]$

In formula: γ is coefficient of heat insulation, S _refBe area of reference, L _refBe reference length, λ _IBe inner wing leading edge sweep, b _0IBe interior wing root chord length, b ₁Be aerofoil tip string, x _sBe the exposed wing installation site distance apart from the bomb body summit.

K _{α α}Be use for reference " missile aerodynamics " (Miao Ruisheng occupies virtuous inscription, Wu Jiasheng. Beijing: National Defense Industry Press, 2006) in wing body interference coefficient that the research of wing body interference effect is introduced.

2) C _{m_L}(α) computing formula.

C _{m_L}(α)＝C _{m_Lzhi1}+C _{m_L_zhi2}+C _{m_wb_L}（3）

Wherein,

$C_{m\_\mathrm{<figure-callout\; id="1"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}1\_L}=\frac{4\mathrm{\&alpha;}{\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">sin</figure-callout>}}^2{\mathrm{\&theta;}}_1}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1{\tan \mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}(\frac{1}{2}{x}_{\mathrm{<figure-callout\; id="1"\; label="c\; L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}{L}_{}^{}{}_1^2-\frac{1}{3}{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1^3)$

$C_{m\_\mathrm{<figure-callout\; id="2"\; label="zhi"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">zhi</figure-callout>}2\_L}=\frac{4\mathrm{\&alpha;}{\sin \mathrm{\&theta;}}_2{\cos \mathrm{\&theta;}}_2}{{({\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1{\tan \mathrm{\&theta;}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}[\frac{1}{2}{L}_1{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2{\tan \mathrm{\&theta;}}_1(2x_c-{2L}_1-L_2)$

$+{{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}^2{\tan \mathrm{\&theta;}}_2(\frac{1}{2}(x_c-L_1)-\frac{1}{3}{L}_2)$

$C_{m\_\mathrm{wb}\_L}=-0.035\mathrm{\&alpha;\&xi;}\frac{(1-{\mathrm{\&eta;}}_t^2)}{L_{\mathrm{ref}}}(x_c-L_{\mathrm{sh}}-R_{\max }{\mathrm{\&eta;}}_t{\cot \mathrm{\&theta;}}_3+\frac{2}{3}\frac{R_{\max }{\cot \mathrm{\&theta;}}_3(1-{\mathrm{\&eta;}}_t^4)}{1-{\mathrm{\&eta;}}_t^3})$

In formula: L ₁Be pointed cone length, L ₂Be truncated cone length, θ ₁Be pointed cone cone angle, θ ₂Be truncated cone cone angle, θ ₃Be the afterbody angle of throat, α is the angle of attack, x _cBe centroid position, L _refBe reference length, η _tBe afterbody shrinkage ratio, R _maxBe the column part radius.

Two, the definition of pneumatic degree of coupling evaluation index.

Pneumatic degree of coupling evaluation index comprises four: the degree of coupling of rudder kick to roll channel

The degree of coupling of rudder kick to pitch channel

The degree of coupling of Jenkel rudder deflection to jaw channel

The degree of coupling of Jenkel rudder deflection to pitch channel

Its definition is respectively:

$\frac{m_x^{{\mathrm{\&delta;}}_y}}{m_x^{{\mathrm{\&delta;}}_x}}=p\frac{{2k}_T{S}_{\mathrm{cw}}{y}_r{\left(K_{\mathrm{\&delta;}0}\right)}_{\mathrm{cwr}}{S}_{\mathrm{cwr}}{\left(K_{\mathrm{\&delta;}0}\right)}_W}{\mathrm{SL}{\left(K_{\mathrm{\&delta;}0}\right)}_{\mathrm{Wr}}{S}_{\mathrm{Wr}}\mathrm{\&sigma;}\frac{({\mathrm{\&eta;}}_k+1){(1-\stackrel{\&OverBar;}{D})}^2}{{\mathrm{\&eta;}}_k+1-2\stackrel{\&OverBar;}{D}}[\stackrel{\&OverBar;}{D}+(1-\stackrel{\&OverBar;}{D})f]}---\left(4\right)$

When aircraft is single vertical fin layout: p=1; When aircraft was twin-finned layout: p was the function about (α, Ma, h).In formula: k _TBe air-flow retardation factor, S _cwBe the vertical fin area, S is aircraft feature area, and L is the aircraft characteristic length, y _rBe the distance of the yaw rudder center of area to the bomb body longitudinal axis, S _cwrBe rudder area, S _WrBe vertical fin area, (K _{δ 0}) _cwrBe the interference coefficient between yaw rudder and wing body, (K _{δ 0}) _WBe wing body interference coefficient, (K _{δ 0}) _WrBe the interference coefficient between vertical fin and wing body, η _kBe contraction coefficient, Be footpath exhibition ratio, f exposes the wing root chord section to the distance of pressing heart other shore (be similar to and get 0.5) with half length, and σ is that the bomb body relative diameter changes the correction factor that rolling moment is exerted an influence and introduces.

$\frac{m_z^{{\mathrm{\&delta;}}_y}}{m_z^{{\mathrm{\&delta;}}_z}}={\mathrm{qa}}_0(1+a_{1}h)(1+a_2\mathrm{\&alpha;})(1+a_3\mathrm{Ma})(1+a_4{\mathrm{\&eta;}}_{\mathrm{kew}})(1+a_5{\mathrm{\&lambda;}}_{\mathrm{kew}})(1+a_6\cos {\mathrm{\&chi;}}_{\mathrm{cw}})---\left(5\right)$

$(1+a_7{\mathrm{\&lambda;}}_{\mathrm{kcwr}})(1+a_8{y}_{\mathrm{cwr}})(1+a_9{L}_{\mathrm{cw}})(1+a_{10}{\stackrel{\&OverBar;}{D}}_{\mathrm{cw}})(1+a_{11}\frac{S_{\mathrm{cwr}}}{S})(1+a_{12}\frac{S_{\mathrm{cw}}}{S})$

When aircraft is single vertical fin layout: q=1; When aircraft was twin-finned layout: q was the function about (α, Ma, h).In formula: S _cw, η _kcw, λ _kcw, χ _cwBe respectively the area of vertical tail, expose the vertical fin contraction coefficient, expose vertical fin aspect ratio and angle of sweep; S _cwr, λ _Kcwr,

L _cwThe area that represents respectively yaw rudder, the yaw rudder aspect ratio, the footpath exhibition of aircraft vertical fin is compared and is gone out the outer crab rooting string mid point in position to the distance of Vehicle nose.

$\frac{m_y^{{\mathrm{\&delta;}}_x}}{m_y^{{\mathrm{\&delta;}}_y}}={\mathrm{mb}}_0(1+b_{1}h)(1+{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_2\mathrm{\&alpha;})(1+b_3\mathrm{Ma})(1+b_4{\mathrm{\&eta;}}_{\mathrm{kw}})(1+b_5{\mathrm{\&lambda;}}_{\mathrm{kw}})(1+b_6\cos {\mathrm{\&chi;}}_w)---\left(6\right)$

$(1+b_7{\mathrm{\&lambda;}}_{\mathrm{kwr}})(1+b_8{y}_{\mathrm{wr}})(1+b_9{L}_w)(1+b_{10}{\stackrel{\&OverBar;}{D}}_{\mathrm{cw}})(1+b_{11}\frac{S_{\mathrm{wr}}}{S})(1+b_{12}\frac{S_w}{S})$

When aircraft is single sweepback layout: m=1; When aircraft was two sweepback layout: m was the function about (α, Ma, h).In formula: S _wThe aerofoil area, η _kwExpose the vertical fin contraction coefficient, λ _kwBe the aspect ratio of the wing, χ _wBe vertical fin angle of sweep, S _wrBe Jenkel rudder area, λ _kwrBe the Jenkel rudder aspect ratio,

Be aircraft missile wing footpath exhibition ratio, L _wFor exposing the wing root chord mid point to the distance of Vehicle nose, y _wrBe the distance of the Jenkel rudder center of area to the bomb body longitudinal axis.

$\frac{m_z^{{\mathrm{\&delta;}}_x}}{m_z^{{\mathrm{\&delta;}}_z}}={\mathrm{nc}}_0(1+c_{1}h)(1+{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}_2\mathrm{\&alpha;})(1+c_3\mathrm{Ma})(1+c_4{\mathrm{\&eta;}}_{\mathrm{kw}})(1+c_5{\mathrm{\&lambda;}}_{\mathrm{kw}})(1+c_6\cos {\mathrm{\&chi;}}_w)---\left(7\right)$

$(1+c_7{\mathrm{\&lambda;}}_{\mathrm{kwr}})(1+c_8{y}_{\mathrm{wr}})(1+c_9{L}_w)(1+c_{10}{\stackrel{\&OverBar;}{D}}_w)(1+c_{11}\frac{S_{\mathrm{wr}}}{S})(1+c_{12}\frac{S_w}{S})$

When aircraft is single sweepback layout: n=1; When aircraft was two sweepback layout: n was the function about (α, Ma, h).In formula: S _wThe aerofoil area, η _kwExpose the vertical fin contraction coefficient, λ _kwBe the aspect ratio of the wing, χ _wBe vertical fin angle of sweep, S _wrBe Jenkel rudder area, λ _kwrBe the Jenkel rudder aspect ratio,

Be aircraft missile wing footpath exhibition ratio, L _wFor exposing the wing root chord mid point to the distance of Vehicle nose, y _wrBe the distance of the Jenkel rudder center of area to the bomb body longitudinal axis.

Coefficient a in formula (5), formula (6) and formula (7) ₀～a ₁₂, b ₀～b ₁₂, c ₀～c ₁₂The match value that is based on a large amount of CFD computational datas and utilizes the nonlinear multivariable least-square fitting approach to calculate.

$\left\{\begin{array}{c}{a}_0=0.08263,a_1=6.138e-05,a_2=-0.08484,a_3=-0.03714,a_4=0.02908\\ a_5=0.04677,a_6=0.05188,a_7=0.01160,a_8=0.05631,a_9=0.00876\\ a_{10}=0.06555,a_{11}=0.59681,a_{12}=0.15504\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_0=1.01895,{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_1=4.422e-05,{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">b</figure-callout>}}_2=6.5758,b_3=-0.03879,b_4=-0.03635\\ b_5=-0.24995,b_6=-1.19032,b_7=-0.05226,b_8=-0.25926,b_9=-0.10399\\ b_{10}=-0.74664,b_{11}=-2.13691,b_{12}=-0.06801\end{array}\right.$

$\left\{\begin{array}{c}{c}_0=-0.0269,{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}_1=-1.2473e-05,{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">c</figure-callout>}}_2=-0.6602,c_3=-0.0121,c_4=0.0124\\ c_5=0.0737,c_6=0.4357,c_7=0.0178,c_8=0.0772,c_9=0.0343\\ c_{10}=0.2224,c_{11}=0.6530,c_{12}=0.0233\end{array}\right.$

Three, determining with control optimization aim function.

The present invention is directed to high-speed aircraft and optimize with control at the aerodynamic configuration at overall trajectory multi-characteristic points place, and aircraft on overall trajectory arbitrarily unique point i place all have one group of pneumatic non-linear and degree of coupling evaluation index:

$\{{\left({\mathrm{DONL}}_{m_z\left(\mathrm{\&alpha;}\right)}\right)}_i,{(m_x^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z})}_i,{(m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z})}_i,{(m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y})}_i,{(m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z})}_i\}$

Wherein, i=1,2...N, N are the sum of unique point on overall trajectory.

Therefore evaluation index requirement owing to carrying out need to satisfying on overall trajectory with the result that control is optimized non-linear and the degree of coupling at unique point arbitrarily place at overall trajectory multi-characteristic points place needs construct with minor function:

$f_1\left(X\right)=\max \left\{{\left|{\mathrm{DONL}}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\},$ i＝1,2...N

$f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\},$ i＝1,2...N

$f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\},$ i＝1,2...N

$f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\},$ i＝1,2...N

$f_5\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\},$ i＝1,2...N

Wherein, X is the vector (θ that contour structures parameter and flight status parameter by aircraft consist of ₁, θ ₂... h, Ma, α) ^T

The objective function that high-speed aircraft overall trajectory multi-characteristic points is optimized with control comprises pneumatic non-linear and coupled characteristic evaluation index, and wherein pneumatic degree of coupling evaluation index by four sub-index constitutes, therefore the invention belongs to the multi-objective nonlinear optimization problem.

The purpose that high-speed aircraft overall trajectory multi-characteristic points is optimized with control is that profile by change of flight device constantly is to reduce aircraft in the pneumatic non-linear and coupled characteristic evaluation index at the place of unique point i arbitrarily.The purpose of optimizing with control according to high-speed aircraft and the pneumatic nonlinearity of aircraft, the definition of degree of coupling evaluation index, the objective function that high-speed aircraft overall trajectory multi-characteristic points is optimized with control is:

${\min f}_1\left(X\right)=\max \left\{{\left|\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\}$

$\min f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\}$

$\min f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}$

$\min f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\}$

$\min f_5\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}$

Wherein, i=1,2...N, N are the sum of unique point on overall trajectory.

Four, with the foundation of controlling Optimized model.

High-speed aircraft overall trajectory multi-characteristic points is as follows with the control Optimized model:

minf ₁(X)

minf ₂(X)

；（8）

minf ₅(X)

s.t.g _i(X)≤0,i＝1,2,...m

Wherein, g _i(X)≤0, i=1,2 ... m is the constraint condition with the control Optimized model, and the Main Function of this constraint condition is that the aerodynamic arrangement of the high-speed aircraft before and after optimizing with control is consistent.

, introduce a high-speed aircraft overall trajectory multi-characteristic points based on " bipyramid+two sweepback+single vertical fin " layout and optimize example with the control Optimized model for the described high-speed aircraft overall trajectory of formula (8) multi-characteristic points.

1) high-speed aircraft contour structures parameter.

The example aircraft contour structures parameter of " bipyramid+two sweepback+single vertical fin " layout is:

θ ₁＝9° θ ₂＝4.5° θ ₃＝45° x _c＝4.5m L ₁＝3m L ₂＝2.5m L ₃＝0.3m L ₄＝0.2m

b _0I＝2.8m b _0II＝2m χ _I＝60° χ _II＝30° b _wr＝0.15m l _wr＝0.5m

χ _cw＝5.7° b _cws＝0.3m b _cwl＝0.25m b _cwr＝0.15m l _cwr＝0.3m

After determining the layout and elementary structure parameter thereof of aircraft, can calculate by given formal parameter value the value of following intermediate variable in order to simplify calculating:

η _t＝0.7023 R _max＝0.6719m S _ref＝1.4183m ² b ₁＝1.9635m l＝2.5R _max

2) choose based on the unique point of overall trajectory.

This example is chosen four unique points altogether on overall trajectory, the flight status parameter at each unique point place is:

Unique point 1: height h=61965m, Mach number Ma=21.74, angle of attack=15 °

Unique point 2: height h=41155m, Mach number Ma=20.95, angle of attack=15 °

Unique point 3: height h=20000m, Mach number Ma=11.14, angle of attack=10 °

Unique point 4: height h=0m, Mach number Ma=6.76, angle of attack=2 °

3) implementation step of optimizing with control.

At first, high-speed aircraft overall trajectory multi-characteristic points is set with the feasible zone of control Optimized model.Constraint condition in formula (8) is set to:

$s.t.=\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="Z"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">Z</figure-callout>}}_0\&CenterDot;(1-0.2)<Z\\ Z<{\mathrm{<figure-callout\; id="0"\; label="Z"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">Z</figure-callout>}}_0\&CenterDot;(1+0.2)\\ b_{0I}+{\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}<{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2+L_3\\ b_{0I}<2.5\&CenterDot;R_{\max }\&CenterDot;{\tan \mathrm{\&chi;}}_I\\ 2.5\&CenterDot;R_{\max }\&CenterDot;{\tan \mathrm{\&chi;}}_I\&CenterDot;{\tan \mathrm{\&chi;}}_{\mathrm{II}}<b_{0I}\&CenterDot;{\tan \mathrm{\&chi;}}_I+{\mathrm{<figure-callout\; id="0"\; label="b"\; filenames="HDA00002853860200011.png"\; state="\{\{state\}\}">b</figure-callout>}}_{0\mathrm{II}}\&CenterDot;{\tan \mathrm{\&chi;}}_I\\ x_c<{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_1+{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">L</figure-callout>}}_2+L_3\\ I_{\mathrm{wr}}<{1.5\&CenterDot;R_{\max }}_{}\\ b_{\mathrm{cwr}}<b_{\mathrm{cwl}}\\ b_{\mathrm{cwl}}<b_{\mathrm{cws}}\\ l_{\mathrm{cwr}}\&CenterDot;{\tan \mathrm{\&chi;}}_{\mathrm{cw}}<b_{\mathrm{cws}}-b_{\mathrm{cwl}}\end{array}\right.$

Wherein, Z is the vector (θ that the contour structures parameter by aircraft consists of ₁, θ ₂..., b _cwr, l _cwr) ^T, Z ₀Expression aircraft original shape parameter vector.

The present invention's reference " non-linear Optimal calculation method " (Zhang Guangcheng. Beijing: Higher Education Publishing House, 2005) theoretical with multiple-objection optimization commonly used at present for the optimum theory of nonlinear problem in, adopt that ideal point method and simulated annealing cooperates mutually with controlling optimized algorithm.

At first, by ideal point method, multi-objective optimization question is converted to single-object problem, is about to a plurality of objective functions and synthesizes the simple target function.Can change the expression formula of rear single-goal function h (X) according to the computing formula of ideal point method:

$h\left(X\right)={\left[\underset{j=1}{\overset{5}{\mathrm{\&Sigma;}}}\frac{{(f_j\left(x\right)-f_j^0)}^2}{f_j^0}\right]}^{1/2}---\left(9\right)$

Wherein, Be function f _j(X) ideal value.

Utilize formula (9) with a plurality of objective function f ₁～f ₅Convert simple target function h to.Wherein, the parameter f in formula (9) ₁ ⁰～f ₅ ⁰Value as follows:

f ₁ ⁰＝1e-4,f ₂ ⁰＝1e-4,f ₃ ⁰＝1e-4,f ₄ ⁰＝1e-2,f ₅ ⁰＝1e-2

According to the ultimate principle of simulated annealing, wherein parameter is arranged.This example adopts permanent initial temperature, fixing temperature decline length, fixing iterative steps and zero degree stop criterion, and its design parameter value is as follows:

Initial temperature T ₀=500 °; T=5 ° of temperature decline length Δ;

Final temperature T=0 °; Iterative steps L=1000

Operation use that VC writes based on ideal point method and simulated annealing with controlling optimizer.Obtain the basic formal parameter value of the optimization solution of objective function h (X) in feasible zone and corresponding high-speed aircraft after the program end of run.

Utilize ideal point method to obtain the as follows with the control Optimized model of single-goal function:

min?h(X)（10）

s.t.g _i(X)≤0,i＝1,2,...m

Secondly, utilize simulated annealing single-goal function with the feasible zone of control Optimized model in search satisfy the optimization solution of pneumatic nonlinearity and the requirement of degree of coupling evaluation index.In simulated annealing utilization simulation thermodynamics, the temperature-fall period of classical particle system is found the solution the extreme value of planning problem, and its basic step is:

1. given cooling program parameter and iteration initial solution X thereof ₀And h (X ₀);

2. parametric t=t _kThe time make L _kInferior heuristic search;

According to current solution X _kProduce a random vector Z _kObtain X _kThe new exploration point X ' of neighborhood _k:

X′ _k＝X _k+Z _k

3. produce one at (0,1) upper equally distributed random number η, calculate at given current iteration point X _kAnd temperature t _kThe lower transition probability P corresponding with the Metropolis acceptance criterion has:

$P=\left\{\begin{array}{cc}1,& h\left(X_k^{\&prime;}\right)\&le;h\left(X_k\right)\\ e^{\frac{h\left(X_k\right)-h\left(X_k^{\&prime;}\right)}{t_k}},& h\left(X_k^{\&prime;}\right)>h\left(X_k\right)\end{array}\right.$

If η≤P accepts new explanation X _k=X ' _k, h (X _k)=h (X ' _k), otherwise X _kConstant.

4. sound out point search less than L _kInferior, turn step 2., otherwise turn step 5.;

5. do you satisfy stopping criterion for iteration?

Be, stop, current solution is Approximate Global Optimal Solution; No, turn step 6..

6. according to given temperature damping's function, produce new temperature control parameter t _k+1And the length L of Markov chain _k+1

7. repeating step 2.-6., until find optimum solution.

Five, interpretation of result.

This example carry out the overall trajectory multi-characteristic points optimize with control before the evaluation index value of pneumatic non-linear and coupled characteristic of high-speed aircraft be:

$\left\{\begin{array}{c}{f}_1\left(X\right)=\max \left\{{\left|\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\}=0.56005\\ f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\}=0.076547\\ f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}=0.073565\\ f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\}=2.100998\\ f_5\left(X\right)=\max \left\{{|{\max }_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}=0.344559\end{array}\right.$

This example carry out the overall trajectory multi-characteristic points optimize with control after the evaluation index value of pneumatic non-linear and coupled characteristic of high-speed aircraft be:

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="f"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">f</figure-callout>}}_1^{*}\left(X\right)=\max \left\{{\left|\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\}=0.037411\\ {\mathrm{<figure-callout\; id="2"\; label="f"\; filenames="HDA00002853860200012.png"\; state="\{\{state\}\}">f</figure-callout>}}_2^{*}\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\}=0.051086\\ f_3^{*}\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}=0.072101\\ f_4^{*}\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\}=0.000849\\ f_5^{*}\left(X\right)=\max \left\{{|{\max }_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}=0.44836\end{array}\right.$

Aircraft contour structures parameter after optimization is:

θ ₁＝9.738° θ ₂＝4.905° θ ₃＝47.52° L ₁＝3.456m L ₂＝2.845m L ₃＝0.2712m L ₄＝0.2316m

x _c＝5.319m b _0I＝2.4416m b _0II＝1.76m χ _I＝67.32° χ _II＝28.5° b _wr＝0.1476m l _wr＝0.5m

χ _cw＝6.1788° b _cws＝0.3414m b _cwl＝0.2935m b _cwr＝0.1599m l _cwr＝0.2412m

Pneumatic non-linear and evaluation index value coupled characteristic before and after high-speed aircraft is optimized with control is compared, and can draw the following conclusions:

A) optimizing rear the pneumatic non-linear of high-speed aircraft with control significantly reduces;

B) optimize the coupled characteristic of rear high-speed aircraft with control: rudder kick reduces the degree of coupling of lift-over, pitch channel, and Jenkel rudder deflection significantly reduces the degree of coupling of jaw channel, and Jenkel rudder deflection increases the degree of coupling of pitch channel;

C) in general, the pneumatic non-linear and coupled characteristic of the high-speed aircraft after optimizing with control be improved significantly.

Claims (1)

1. One kind based on pneumatic non-linear with the coupling high-speed aircraft with the control optimization method, it is characterized in that comprising the following steps:

   Step 1, be the pneumatic nonlinearity of pitch channel with pneumatic non-linear evaluation index definition , namely

   $\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}=\frac{|C_m\left(\mathrm{\&alpha;}\right)-C_{m\_L}\left(\mathrm{\&alpha;}\right)|}{|C_m\left(\mathrm{\&alpha;}\right)-C_{m\_L}\left(\mathrm{\&alpha;}\right)|+\left|C_{m\_L}\left(\mathrm{\&alpha;}\right)\right|}---\left(1\right)$

   $=\frac{|m_z\left(\mathrm{\&alpha;}\right)-m_z\left(0\right)-{m_z}^{\&prime;}\left(0\right)\mathrm{\&alpha;}|}{|m_z\left(\mathrm{\&alpha;}\right)-m_z\left(0\right)-{m_z}^{\&prime;}\left(0\right)\mathrm{\&alpha;}|+|m_z\left(0\right)+{m_z}^{\&prime;}\left(0\right)\mathrm{\&alpha;}|}$

   In formula, C _m(α) be the stable pitching moment coefficient of aircraft; C _{m_L}(α) be C _mLinear segment (α);

   1) C _m(α) computing formula;

   C _m(α)=C _{Bomb body}+ K _{α α}C _Aerofoil=C _{M_zhi1}+ C _{M_zhi2}+ C _{M_zhu}+ K _{α α}C _{m_w}(2)

   Wherein,

   $C_{m\_\mathrm{zhi}1}=\frac{4\sin \mathrm{\&alpha;}\cos \mathrm{\&alpha;si}{n}^2{\mathrm{\&theta;}}_1}{{(L_1\tan {\mathrm{\&theta;}}_1+L_2\tan {\mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}(\frac{1}{2}{x}_c{L}_1^2-\frac{1}{3}{L}_1^3)$

   $C_{m\_\mathrm{zhi}2}=\frac{4\sin \mathrm{\&alpha;}\cos \mathrm{\&alpha;si}n{\mathrm{\&theta;}}_2\cos {\mathrm{\&theta;}}_2}{{(L_1\tan {\mathrm{\&theta;}}_1+L_2\tan {\mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}[L_1{L}_2{\tan \mathrm{\&theta;}}_1(x_c-L_1-\frac{1}{2}{L}_2)+{L_2}^2{\tan \mathrm{\&theta;}}_2(\frac{1}{2}(x_c-L_1)-\frac{1}{3}{L}_2)]$

   $C_{m\_\mathrm{zhu}}=\frac{8{\sin }^2\mathrm{\&alpha;}{L}_3[x_c-L_3/2-(L_1+L_2)}{3\mathrm{\&pi;}(L_1{\tan \mathrm{\&theta;}}_1+L_2{\tan \mathrm{\&theta;}}_2)L_{\mathrm{ref}}}$

   In formula: L ₁Be pointed cone length, L ₂Be truncated cone length, L ₃Be shell of column length, θ ₁Be pointed cone cone angle, θ ₂Be the truncated cone cone angle, α is the angle of attack, x _cBe centroid position, L _refBe reference length;

   When the flat shape of the wing is two swept-back wing, its stable pitching moment coefficient C _{m_w}For:

   $C_{m\_w}=\frac{(\mathrm{\&gamma;}+1){\mathrm{\&alpha;}}^2}{S_{\mathrm{ref}}{L}_{\mathrm{ref}}}\sqrt{1+{\left[\frac{4}{(\mathrm{\&gamma;}+1)M_{\&infin;}\mathrm{\&alpha;}}\right]}^2}\left\{{\cot \mathrm{\&lambda;}}_I\right[{b_{0I}}^2(x_c-x_s)-\frac{2}{3}{b_{0I}}^3]-$

   $\frac{2}{3}\cot {\mathrm{\&lambda;}}_{\mathrm{II}}{(b_{0\mathrm{II}}-b_1)}^3+[\cot {\mathrm{\&lambda;}}_{\mathrm{II}}(x_c-x_s)-b_{0I}({\cot \mathrm{\&lambda;}}_I+{\cot \mathrm{\&lambda;}}_{\mathrm{II}})]{(b_{0\mathrm{II}}-b_1)}^2$

   $+{2b}_{0I}{\cot \mathrm{\&lambda;}}_I(b_{0\mathrm{II}}-b_1)(x_c-x_s-b_{0I})$

   $+[(b_{0\mathrm{II}}-b_1){\cot \mathrm{\&lambda;}}_{\mathrm{II}}+b_{0I}{\cot \mathrm{\&lambda;}}_I][2(x_c-b_{0I}-b_{0\mathrm{II}}-x_s)b_1+{b_1}^2]\}$

   In formula: γ is coefficient of heat insulation, S _refBe area of reference, L _refBe reference length, λ _IBe inner wing leading edge sweep, λ _IIBe outer wing leading edge sweep, b _0IBe interior wing root chord length, b _0IIBe outer wing root chord length, b ₁Be aerofoil tip string, x _sBe the exposed wing installation site distance apart from the bomb body summit;

   When the flat shape of the wing is single swept-back wing, its stable pitching moment coefficient C _{m_w}For:

   $C_{m\_w}=\frac{M_z}{q_{\&infin;}{S}_{\mathrm{ref}}{L}_{\mathrm{ref}}}=\frac{(\mathrm{\&gamma;}+1){\mathrm{\&alpha;}}^2}{S_{\mathrm{ref}}{L}_{\mathrm{ref}}}{\cot \mathrm{\&lambda;}}_I\sqrt{1+{\left[\frac{4}{(\mathrm{\&gamma;}+1)M_{\&infin;}\mathrm{\&alpha;}}\right]}^2}[{(b_{0I}-b_1)}^2(x_c-x_s)$

   $-\frac{2}{3}{(b_{0I}-b_1)}^3+{2b}_1(b_{0I}-b_1)(x_c-x_s-b_{0I}+b_1)-(b_{0I}-b_1){b_1}^2]$

   In formula: γ is coefficient of heat insulation, S _refBe area of reference, L _refBe reference length, λ _IBe inner wing leading edge sweep, b _0IBe interior wing root chord length, b ₁Be aerofoil tip string, x _sBe the exposed wing installation site distance apart from the bomb body summit, K _{α α}It is wing body interference coefficient;

   2) C _{m_L}(α) computing formula;

   C _{m_L}(α)＝C _{m_L_zhi1}+C _{m_L_zhi2}+C _{m_wb_L}（3）

   Wherein,

   $C_{m\_\mathrm{zhi}1\_L}=\frac{4\mathrm{\&alpha;}{\sin }^2{\mathrm{\&theta;}}_1}{{(L_1{\tan \mathrm{\&theta;}}_1+L_2{\tan \mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}(\frac{1}{2}{x}_c{L}_1^2-\frac{1}{3}{L}_1^3)$

   $C_{m\_\mathrm{zhi}2\_L}=\frac{4\mathrm{\&alpha;}{\sin \mathrm{\&theta;}}_2{\cos \mathrm{\&theta;}}_2}{{(L_1{\tan \mathrm{\&theta;}}_1+L_2{\tan \mathrm{\&theta;}}_2)}^2{L}_{\mathrm{ref}}}[\frac{1}{2}{L}_1{L}_2{\tan \mathrm{\&theta;}}_1(2x_c-{2L}_1-L_2)$

   $+{L_2}^2{\tan \mathrm{\&theta;}}_2(\frac{1}{2}(x_c-L_1)-\frac{1}{3}{L}_2)$

   $C_{m\_\mathrm{wb}\_L}=-0.035\mathrm{\&alpha;\&xi;}\frac{(1-{\mathrm{\&eta;}}_t^2)}{L_{\mathrm{ref}}}(x_c-L_{\mathrm{sh}}-R_{\max }{\mathrm{\&eta;}}_t{\cot \mathrm{\&theta;}}_3+\frac{2}{3}\frac{R_{\max }{\cot \mathrm{\&theta;}}_3(1-{\mathrm{\&eta;}}_t^4)}{1-{\mathrm{\&eta;}}_t^3})$

   In formula: L ₁Be pointed cone length, L ₂Be truncated cone length, θ ₁Be pointed cone cone angle, θ ₂Be truncated cone cone angle, θ ₃Be the afterbody angle of throat, α is the angle of attack, x _cBe centroid position, L _refBe reference length, η _tBe afterbody shrinkage ratio, R _maxBe the column part radius;

   The degree of coupling of the rudder kick of step 2, the pneumatic degree of coupling evaluation index of definition to roll channel

   

   $\frac{m_x^{{\mathrm{\&delta;}}_y}}{m_x^{{\mathrm{\&delta;}}_x}}=p\frac{{2k}_T{S}_{\mathrm{cw}}{y}_r{\left(K_{\mathrm{\&delta;}0}\right)}_{\mathrm{cwr}}{S}_{\mathrm{cwr}}{\left(K_{\mathrm{\&delta;}0}\right)}_W}{\mathrm{SL}{\left(K_{\mathrm{\&delta;}0}\right)}_{\mathrm{Wr}}{S}_{\mathrm{Wr}}\mathrm{\&sigma;}\frac{({\mathrm{\&eta;}}_k+1){(1-\stackrel{\&OverBar;}{D})}^2}{{\mathrm{\&eta;}}_k+1-2\stackrel{\&OverBar;}{D}}[\stackrel{\&OverBar;}{D}+(1-\stackrel{\&OverBar;}{D})f]}---\left(4\right)$

   In formula, when aircraft is single vertical fin layout, p=1; When aircraft was twin-finned layout, p was the function about (α, Ma, h); k _TBe air-flow retardation factor, S _cwBe the vertical fin area, S is aircraft feature area, and L is the aircraft characteristic length, y _rBe the distance of the yaw rudder center of area to the bomb body longitudinal axis, S _cwrBe rudder area, S _WrBe vertical fin area, (K _{δ 0}) _cwrBe the interference coefficient between yaw rudder and wing body, (K _{δ 0}) _WBe wing body interference coefficient, (K _{δ 0}) _WrBe the interference coefficient between vertical fin and wing body, η _kBe contraction coefficient,

   

   Be footpath exhibition ratio, f exposes the wing root chord section to the distance of pressing heart other shore with half length, and σ is that the bomb body relative diameter changes the correction factor that rolling moment is exerted an influence and introduces;

   Define the rudder kick of pneumatic degree of coupling evaluation index to the degree of coupling of pitch channel

   

   $\frac{m_z^{{\mathrm{\&delta;}}_y}}{m_z^{{\mathrm{\&delta;}}_z}}={\mathrm{qa}}_0(1+a_{1}h)(1+a_2\mathrm{\&alpha;})(1+a_3\mathrm{Ma})(1+a_4{\mathrm{\&eta;}}_{\mathrm{kew}})(1+a_5{\mathrm{\&lambda;}}_{\mathrm{kew}})(1+a_6\cos {\mathrm{\&chi;}}_{\mathrm{cw}})---\left(5\right)$

   $(1+a_7{\mathrm{\&lambda;}}_{\mathrm{kcwr}})(1+a_8{y}_{\mathrm{cwr}})(1+a_9{L}_{\mathrm{cw}})(1+a_{10}{\stackrel{\&OverBar;}{D}}_{\mathrm{cw}})(1+a_{11}\frac{S_{\mathrm{cwr}}}{S})(1+a_{12}\frac{S_{\mathrm{cw}}}{S})$

   In formula, when aircraft is single vertical fin layout, q=1; When aircraft was twin-finned layout, q was the function about (α, Ma, h); S _cw, η _kcw, λ _kcw, χ _cwBe respectively the area of vertical tail, expose the vertical fin contraction coefficient, expose vertical fin aspect ratio and angle of sweep; S _cwr, λ _Kcwr,

   

   L _cwThe area that represents respectively yaw rudder, the yaw rudder aspect ratio, the footpath exhibition of aircraft vertical fin is compared and is gone out the outer crab rooting string mid point in position to the distance of Vehicle nose;

   Define the Jenkel rudder deflection of pneumatic degree of coupling evaluation index to the degree of coupling of jaw channel

   

   $\frac{m_y^{{\mathrm{\&delta;}}_x}}{m_y^{{\mathrm{\&delta;}}_y}}={\mathrm{mb}}_0(1+b_{1}h)(1+b_2\mathrm{\&alpha;})(1+b_3\mathrm{Ma})(1+b_4{\mathrm{\&eta;}}_{\mathrm{kw}})(1+b_5{\mathrm{\&lambda;}}_{\mathrm{kw}})(1+b_6\cos {\mathrm{\&chi;}}_w)---\left(6\right)$

   $(1+b_7{\mathrm{\&lambda;}}_{\mathrm{kwr}})(1+b_8{y}_{\mathrm{wr}})(1+b_9{L}_w)(1+b_{10}{\stackrel{\&OverBar;}{D}}_{\mathrm{cw}})(1+b_{11}\frac{S_{\mathrm{wr}}}{S})(1+b_{12}\frac{S_w}{S})$

   In formula, when aircraft is single sweepback layout, m=1; When aircraft was two sweepback layout, m was the function about (α, Ma, h); S _wThe aerofoil area, η _kwExpose the vertical fin contraction coefficient, λ _kwBe the aspect ratio of the wing, χ _wBe vertical fin angle of sweep, S _wrBe Jenkel rudder area, λ _kwrBe the Jenkel rudder aspect ratio,

   

   Be aircraft missile wing footpath exhibition ratio, L _wFor exposing the wing root chord mid point to the distance of Vehicle nose, y _wrBe the distance of the Jenkel rudder center of area to the bomb body longitudinal axis;

   Define the Jenkel rudder deflection of pneumatic degree of coupling evaluation index to the degree of coupling of pitch channel

   

   $\frac{m_z^{{\mathrm{\&delta;}}_x}}{m_z^{{\mathrm{\&delta;}}_z}}={\mathrm{nc}}_0(1+c_{1}h)(1+c_2\mathrm{\&alpha;})(1+c_3\mathrm{Ma})(1+c_4{\mathrm{\&eta;}}_{\mathrm{kw}})(1+c_5{\mathrm{\&lambda;}}_{\mathrm{kw}})(1+c_6\cos {\mathrm{\&chi;}}_w)---\left(7\right)$

   $(1+c_7{\mathrm{\&lambda;}}_{\mathrm{kwr}})(1+c_8{y}_{\mathrm{wr}})(1+c_9{L}_w)(1+c_{10}{\stackrel{\&OverBar;}{D}}_w)(1+c_{11}\frac{S_{\mathrm{wr}}}{S})(1+c_{12}\frac{S_w}{S})$

   In formula, when aircraft is single sweepback layout, n=1; When aircraft was two sweepback layout, n was the function about (α, Ma, h); S _wThe aerofoil area, η _kwExpose the vertical fin contraction coefficient, λ _kwBe the aspect ratio of the wing, χ _wBe vertical fin angle of sweep, S _wrBe Jenkel rudder area, λ _kwrBe the Jenkel rudder aspect ratio,

   

   Be aircraft missile wing footpath exhibition ratio, L _wFor exposing the wing root chord mid point to the distance of Vehicle nose, y _wrBe the distance of the Jenkel rudder center of area to the bomb body longitudinal axis;

   Coefficient a in formula (5), formula (6) and formula (7) ₀～a ₁₂, b ₀～b ₁₂, c ₀～c ₁₂The match value that is based on a large amount of CFD computational datas and utilizes the nonlinear multivariable least-square fitting approach to calculate;

   $\left\{\begin{array}{c}{a}_0=0.08263,a_1=6.138e-05,a_2=-0.08484,a_3=-0.03714,a_4=0.02908\\ a_5=0.04677,a_6=0.05188,a_7=0.01160,a_8=0.05631,a_9=0.00876\\ a_{10}=0.06555,a_{11}=0.59681,a_{12}=0.15504\end{array}\right.$

   $\left\{\begin{array}{c}{b}_0=1.01895,b_1=4.422e-05,b_2=6.5758,b_3=-0.03879,b_4=-0.03635\\ b_5=-0.24995,b_6=-1.19032,b_7=-0.05226,b_8=-0.25926,b_9=-0.10399\\ b_{10}=-0.74664,b_{11}=-2.13691,b_{12}=-0.06801\end{array}\right.$

   $\left\{\begin{array}{c}{c}_0=-0.0269,c_1=-1.2473e-05,c_2=-0.6602,c_3=-0.0121,c_4=0.0124\\ c_5=0.0737,c_6=0.4357,c_7=0.0178,c_8=0.0772,c_9=0.0343\\ c_{10}=0.2224,c_{11}=0.6530,c_{12}=0.0233\end{array}\right.$

   Step 3, aircraft on overall trajectory arbitrarily unique point i place all have one group of pneumatic non-linear and degree of coupling evaluation index:

   $\{{\left({\mathrm{DONL}}_{m_z\left(\mathrm{\&alpha;}\right)}\right)}_i,{(m_x^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z})}_i,{(m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z})}_i,{(m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y})}_i,{(m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z})}_i\}$

   Wherein, i=1,2...N, N are the sum of unique point on overall trajectory;

   Therefore evaluation index requirement owing to carrying out need to satisfying on overall trajectory with the result that control is optimized non-linear and the degree of coupling at unique point arbitrarily place at overall trajectory multi-characteristic points place needs construct with minor function:

   $f_1\left(X\right)=\max \left\{{\left|{\mathrm{DONL}}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\},$ i＝1,2...N

   $f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\},$ i＝1,2...N

   $f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\},$ i＝1,2...N

   $f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\},$ i＝1,2...N

   $f_5\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\},$ i＝1,2...N

   In formula, X is the vector (θ that contour structures parameter and the flight status parameter by aircraft consists of ₁, θ ₂... h, Ma, α) ^T

   The purpose of optimizing with control according to high-speed aircraft and the pneumatic nonlinearity of aircraft, the definition of degree of coupling evaluation index, the objective function that high-speed aircraft overall trajectory multi-characteristic points is optimized with control is:

   ${\min f}_1\left(X\right)=\max \left\{{\left|\mathrm{DON}{L}_{m_z\left(\mathrm{\&alpha;}\right)}\right|}_i\right\}$

   $\min f_2\left(X\right)=\max \left\{{|m_x^{{\mathrm{\&delta;}}_y}/m_x^{{\mathrm{\&delta;}}_x}|}_i\right\}$

   $\min f_3\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_y}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}$

   $\min f_4\left(X\right)=\max \left\{{|m_y^{{\mathrm{\&delta;}}_x}/m_y^{{\mathrm{\&delta;}}_y}|}_i\right\}$

   $\min f_5\left(X\right)=\max \left\{{|m_z^{{\mathrm{\&delta;}}_x}/m_z^{{\mathrm{\&delta;}}_z}|}_i\right\}$

   Wherein, i=1,2...N, N are the sum of unique point on overall trajectory;

   Step 4, set up high-speed aircraft overall trajectory multi-characteristic points with the control Optimized model:

   min?f ₁(X)

   min?f ₂(X)

   ；（8）

   min?f ₅(X)

   s.t.g _i(X)≤0,i＝1,2,...m

   In formula, g _i(X)≤0, i=1,2 ... m is the constraint condition with the control Optimized model;

   By ideal point method, multi-objective optimization question is converted to single-object problem, is about to a plurality of objective functions and synthesizes the simple target function; Can change the expression formula of rear single-goal function h (X) according to the computing formula of ideal point method:

   $h\left(X\right)={\left[\underset{j=1}{\overset{5}{\mathrm{\&Sigma;}}}\frac{{(f_j\left(x\right)-f_j^0)}^2}{f_j^0}\right]}^{1/2}---\left(9\right)$

   Wherein,

   

   Be function f _j(X) ideal value;

   Utilize ideal point method to obtain the as follows with the control Optimized model of single-goal function:

   min?h(X)（10）

   s.t.g _i(X)≤0,i＝1,2,...m

   Utilize simulated annealing single-goal function with the feasible zone of control Optimized model in search satisfy the optimization solution of pneumatic nonlinearity and degree of coupling evaluation index requirement; In simulated annealing utilization simulation thermodynamics, the temperature-fall period of classical particle system is found the solution the extreme value of planning problem, and its basic step is:

   1. given cooling program parameter and iteration initial solution X thereof ₀And h (X ₀);

   2. parametric t=t _kThe time make L _kInferior heuristic search;

   According to current solution X _kProduce a random vector Z _kObtain X _kThe new exploration point X ' of neighborhood _k:

   X′ _k＝X _k+Z _k

   3. produce one at (0,1) upper equally distributed random number η, calculate at given current iteration point X _kAnd temperature t _kThe lower transition probability P corresponding with the Metropolis acceptance criterion has:

   $P=\left\{\begin{array}{cc}1,& h\left(X_k^{\&prime;}\right)\&le;h\left(X_k\right)\\ e^{\frac{h\left(X_k\right)-h\left(X_k^{\&prime;}\right)}{t_k}},& h\left(X_k^{\&prime;}\right)>h\left(X_k\right)\end{array}\right.$

   If η≤P accepts new explanation X _k=X ' _k, h (X _k)=h (X ' _k), otherwise X _kConstant;

   4. sound out point search less than L _kInferior, turn step 2., otherwise turn step 5.;

   5. do you satisfy stopping criterion for iteration?

   Be, stop, current solution is Approximate Global Optimal Solution; No, turn step 6.;

   6. according to given temperature damping's function, produce new temperature control parameter t _k+1And the length L of Markov chain _k+17. repeating step 2.-6., until find optimum solution.

Images