CN103809450A - Multi-loop aircraft model cluster flutter restraining composite root locus multistage PID (Proportion Integration Differentiation) robust controller design method - Google Patents

Abstract

The invention provides a multi-loop aircraft model cluster flutter restraining composite root locus multistage PID (Proportion Integration Differentiation) robust controller design method. The method comprises the following steps: directly determining to obtain a model cluster constructed by amplitude frequency and phase frequency features in a full envelope by means of frequency sweeping flight test under the condition that different heights and Mach numbers are given; giving a closed loop pole distribution limitation index under corresponding root locus description through multi-loop model equivalence according to an amplitude frequency margin in the flight envelope and a phase margin military standard requirement, determining the stage number and parameter value of a multistage PID robust controller by adding the multistage PID robust controller, the closed-loop pole distribution limitation index in the full envelope of an air craft and a model identification method in system identification; designing a low-altitude flight aircraft which is accordant with the full envelope, can be used for restraining flutter, has low overshot and is stable on the basis of the concept of closed-loop pole distribution limitation under root locus description.

Description

The multistage PID robust Controller Design of the compound root locus of multiloop dummy vehicle bunch Flutter Suppression method

Technical field

The present invention relates to a kind of controller of aircraft method for designing, particularly the multistage PID robust Controller Design of the compound root locus of multiloop dummy vehicle bunch Flutter Suppression method, belongs to the category such as observation and control technology and flight mechanics.

Background technology

The control of aircraft landing process plays an important role to flight safety; Because flying speed in aircraft landing process changes greatly, even also can face strong nonlinearity problem according to longitudinal model; On the other hand, there is the phenomenons such as saturated, dead band in the control vane of aircraft; Consider from flight safety, when hedgehopping (as take off/land), controller must guarantee that system has certain stability margin, non-overshoot and stationarity, like this, just make hedgehopping controller design very complicated, can not directly apply mechanically existing control theory and carry out the design of aircraft control.

In the design of modern practical flight controller, a small part adopts state-space method to design, and great majority still adopt the classical frequency domain method take PID as representative and carry out controller design against Nyquist Array Method as the modern frequency method of representative.Modern control theory is take state-space method as feature, take analytical Calculation as Main Means, to realize performance index as optimum modern control theory, then have and developed method for optimally controlling, model reference control method, self-adaptation control method, dynamic inversion control method, feedback linearization method, directly nonlinear optimization control, variable-gain control method, neural network control method, fuzzy control method, a series of controller design methods such as robust control method and several different methods combination control, the scientific paper of delivering is ten hundreds of, for example Ghasemi A in 2011 has designed reentry vehicle (the Ghasemi A of Adaptive Fuzzy Sliding Mode Control, Moradi M, Menhaj M B.Adaptive Fuzzy Sliding Mode Control Design for a Low-Lift Reentry Vehicle[J] .Journal of Aerospace Engineering, 2011, 25 (2): 210-216), Babaei A R in 2013 is that non-minimum phase and Nonlinear Flight device have designed fuzzy sliding mode tracking control robot pilot (Babaei A R, Mortazavi M, Moradi M H.Fuzzy sliding mode autopilot design for nonminimum phase and nonlinear UAV[J] .Journal of Intelligent and Fuzzy Systems, 2013, 24 (3): 499-509), a lot of research only rests on the Utopian simulation study stage, and there are three problems in this design: (1), owing to cannot carrying out the extreme low-altitude handling and stability experiment of aircraft, is difficult to obtain the mathematical model of accurate controlled device, (2) stability margin stipulating for army's mark etc. is evaluated the important performance indexes of flight control system, and state-space method far can be expressed with obvious form unlike classical frequency method, (3) too complicated, the constraint of not considering working control device and state of flight of controller architecture, the controller of design physically can not be realized.

The scholar Rosenbrock of Britain systematically, study in a creative way in the design that how frequency domain method is generalized to multi-variable system and gone, utilize matrix diagonal dominance concept, Multivariable is converted into the design problem of the single-variable system of the classical approach that can know with people, in succession there is Mayne sequence return difference method later, MacFarlane System with Characteristic Locus Method, the methods such as Owens dyadic expansion, common feature is many input more than one outputs, the design of serious associated multi-variable system between loop, turn to the design problem of a series of single-variable systems, and then can select a certain classical approach (frequency response method of Nyquist and Bode, the root-locus technique of Evans etc.) design of completion system, above-mentioned these methods retain and have inherited the advantage of classic graphic-arts technique, do not require accurate especially mathematical model, easily meet the restriction in engineering.Particularly, in the time that employing has the conversational computer-aided design system of people's one machine of graphic display terminal to realize, can give full play to deviser's experience and wisdom, design and both meet quality requirements, be again controller physically attainable, simple in structure; (tall and big far away, sieve becomes, Shen Hui, Hu Dewen, Flexible Satellite Attitude Decoupling Controller Design Using Multiple Variable Frequency Domain Method, aerospace journal, 2007, Vol.28 (2), pp442-447 multivariate frequency method have been carried out improving research both at home and abroad; Xiong Ke, Xia Zhixun, Guo Zhenyun, the hypersonic cruise vehicle multivariable frequency domain approach of banked turn Decoupling design, plays arrow and guidance journal, 2011, Vol.31 (3), pp25-28) still, when this method for designing can taking into account system uncertain problem, conservative property is excessive, under aircraft control vane limited case, can not obtain rational design result; Particularly, in the time of aircraft generation flutter, designed control system is likely difficult to the stability of the system that guarantees.

In sum, current control method can't change at dummy vehicle, design and can suppress that flutter, non-overshoot are little, low-latitude flying controller stably according to the stability margin index in full flight envelope.

Summary of the invention

The technological deficiency that can not design the stability margin index meeting in full flight envelope at aircraft in the situation that full flight envelope inner model changes greatly and can suppress little, the steady low-latitude flying controller of overshoot of flutter in order to overcome existing method, the invention provides the multistage PID robust Controller Design of the compound root locus of a kind of multiloop dummy vehicle bunch Flutter Suppression method, the method directly determines by frequency sweep flight test the model cluster that the amplitude-frequency that obtains in full envelope curve and phase-frequency characteristic form under given differing heights, Mach number condition; According to the amplitude-frequency nargin in flight envelope and the mark requirement of phase margin army, by multiloop model equivalence, provide the Distribution of Closed Loop Poles restriction index under corresponding root locus description, by adding the identification Method in multistage PID controller the restriction index of the Distribution of Closed Loop Poles in the full envelope curve of aircraft and System Discrimination to determine multistage PID robust controller sum of series parameter value; From root locus describe concept that Distribution of Closed Loop Poles limits design meet full flight envelope can suppress that flutter, overshoot are little, low-latitude flying controller stably.

The technical solution adopted for the present invention to solve the technical problems: the multistage PID robust Controller Design of the compound root locus of a kind of multiloop dummy vehicle bunch Flutter Suppression method, is characterized in comprising the following steps:

Under step 1, given differing heights, Mach number by frequency sweep flight test directly by allowing amplitude-frequency and phase-frequency characteristic in the full envelope curve of flight to form primary control surface in the full envelope curve of aircraft and the model cluster of flying height, and can cross over flight envelope and obtain the flutter frequency of aircraft, obtain open-loop transfer function model cluster matrix between corresponding aircraft primary control surface and flying height and be:

Wherein, G is m × m square formation, and m>1 is positive integer, the independent variable that s is Laplace transformation, h is aircraft altitude, and M is Mach number, and Δ is uncertain vector, P is m × m single mode square formation, and D is m × m polynomial expression diagonal matrix, and Q is m × m single mode square formation

for polynomial expression, n>1 is positive integer;

Choose

satisfy condition:

and

Wherein, G _efor m × m square formation, P _efor m × m single mode square formation, D _efor m × m polynomial expression diagonal matrix, d _i,Efor D _ei _erow, I _ecolumn element,

for the I of D _erow, I _ecolumn element, I _e=1,2 ..., m, Q _efor m × m single mode square formation,

for polynomial expression, arg is phase angle mathematic sign;

The controller of aircraft multiloop system is made as:

Wherein, G _cA(s) be m × m square formation, G _a0(s)=diag[G _{c, 1}(s), G _{c, 2}(s) ..., G _c,m(s)] be m × m diagonal matrix; (s) be G _a0(s) I _erow, I _ecolumn element, I _e=1,2 ..., m;

Step 2, controller

(s), I _e=1,2 ..., the design process of m is as follows:

(1) order

the form of embodying is:

$G_{I_E}\left(s\right)=\frac{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)}{B(h,M,s)}[1+\mathrm{\Δk}\left(s\right)]$

Flutter frequency is: ω _aSE(h, M);

Wherein

A(h,M,s)=s ^m+a _m-1(h,M)s ^m-1+a _m-2(h,M)s ^m-2+…+a ₁(h,M)s+a ₀(h,M)、

B (h, M, s)=s ⁿ+ b _n-1(h, M) s ^n-1+ b _n-2(h, M) s ^n-2+ ... + b ₁(h, M) s+b ₀(h, M) is polynomial expression, and s is the variable after laplace transform conventional in transport function, h, and M is respectively flying height and Mach number, and σ (h, M) is the time delay of pitch channel, and K (h, M) is with h, the gain that M changes, a _l(h, M), l=0,1,2 ..., m-1 be in polynomial expression A (h, M, s) with h, M change coefficient bunch, b _i(h, M), i=0,1,2 ..., n-1 be in polynomial expression B (h, M, s) with h, M change coefficient bunch, the indeterminate that △ k (s) is model;

(2) transport function of the multistage PID controller of candidate is:

$G_{c,I_E}\left(s\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Π}}}[K_p\left(i\right)+k_I\left(i\right)/s+k_D\left(i\right)\·s]}{\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1}$

In formula, k _cfor constant gain to be determined, N is integer, represents the progression of multistage PID controller to be determined, k _p(i), k _i(i), k _d(i) i=1,2 ..., N is constant to be determined, k _aSEfor Flutter Suppression gain;

Add after multistage PID controller, the open-loop transfer function of whole system is:

$G_{I_E}\left(s\right)G_{c,I_E}\left(s\right)=\frac{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)}{B(h,M,s)}[1+\mathrm{\Δk}\left(s\right)]\frac{\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_p\left(i\right){+k}_I\left(i\right)/s+k_D\left(i\right)\·s]}{\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1}$

Corresponding root locus equation is:

$e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Λk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_p\left(i\right)\·s+k_1\left(i\right)+k_D\·s^2]+s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1)=0;$

(3) establish s=σ+j ω, wherein: the real part that σ is s, the imaginary part that ω is s, j is the imaginary part of symbol; The stability margin index of system is set as: σ≤-

2,

wherein,

for non-zero real, ξ gives fixed number; Set up according to flight test or wind tunnel test the lagging phase angle that model indeterminate causes

radian, amplitude the stability margin index of system is adjusted into:

with

wherein, △ _mand △ _abe whole real number;

Like this, the stability margin index of system can be converted into: according to

$\left\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)\right[1+\mathrm{\Δk}\left(s\right)\left]\underset{i=1}{\overset{N}{\mathrm{\Π}}}\right[k_p\left(i\right)\·s+k_I\left(i\right)+k_D\left(i\right)\·s^2]+s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}=0$

Or $\left\{\begin{array}{c}\mathrm{Re}\left\{\right\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Δk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_P\left(i\right)\·s+k_I\left(i\right){+k}_D\left(i\right)\·s^2]+\\ s\·B(h,M,s)\·(\frac{k_{\mathrm{SAE}}}{{\mathrm{\ω}}_{\mathrm{SAE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}\\ \mathrm{Im}\left\{\right\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Δk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_P\left(i\right)\·s+k_I\left(i\right)+k_I\left(i\right)+k_D\left(i\right)\·s^2]+\\ s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}\}=0\end{array}\right.$ The root locus obtaining must meet

with

according under this index and maximum likelihood criterion or the common constraint of other criterion, can determine according to the maximum likelihood method in system model Structure Identification or discrimination method progression N, the constant k of multistage PID controller _p(i), k _i(i), k _d(i) i=1,2 ..., N and Flutter Suppression gain k _aSE.

The invention has the beneficial effects as follows: the concept of the Distribution of Closed Loop Poles restriction from root locus is described, by adding multistage PID controller, in full flight envelope according to meeting that given Distribution of Closed Loop Poles restriction requires and identification Method is determined the parameter of multistage PID robust controller, design meet full flight envelope can suppress that flutter, overshoot are little, low-latitude flying controller stably.

Below in conjunction with embodiment, the present invention is elaborated.

Embodiment

Under step 1, given differing heights, Mach number, use Linear chirp

(f ₀for initial frequency, f ₁for cutoff frequency, r=(f ₁-f ₀)/T, T is the frequency sweep time) or logarithm swept-frequency signal f (t)=A (t) sin{2 π f ₀/ r[exp (rt)-1] } (f ₀for initial frequency, f ₁for cutoff frequency, r=ln (f ₁/ f ₀)/T, T is the frequency sweep time) aircraft is encouraged, amplitude-frequency and phase-frequency characteristic in the full envelope curve that can directly obtain allowing to fly, and can cross over flight envelope and obtain the flutter frequency of aircraft, obtain open-loop transfer function model cluster matrix between corresponding aircraft primary control surface and flying height and be:

Wherein, G is m × m square formation, and m>1 is positive integer, the independent variable that s is Laplace transformation, h is aircraft altitude, and M is Mach number, is uncertain vector, and P is m × m single mode square formation, D is m × m polynomial expression diagonal matrix, and Q is m × m single mode square formation

for polynomial expression, n>1 is positive integer;

Choose satisfy condition:

and

Wherein, G _efor m × m square formation, P _efor m × m single mode square formation, D _efor m × m polynomial expression diagonal matrix, di, E is D _ei _erow, I _ecolumn element, for the I of D _erow, I _ecolumn element, I _e=1,2 ..., m, Q _efor m × m single mode square formation, for polynomial expression, arg is phase angle mathematic sign;

The controller of aircraft multiloop system is made as:

Wherein, G _cA(s) be m × m square formation, G _a0(s)=diag[G _{c, 1}(s), G _{c, 2}(s) ..., G _c,m(s)] be m × m diagonal matrix;

(s) be G _a0(s) I _erow, I _ecolumn element, I _e=1,2 ..., m;

Step 2, controller

(s), I _e=1,2 ..., the design process of m is as follows:

(1) order

the form of embodying is:

$G_{I_E}\left(s\right)=\frac{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)}{B(h,M,s)}[1+\mathrm{\Δk}\left(s\right)]$

With flutter frequency ω _aSE(h, M);

Wherein

A(h,M,s)=s ^m+a _m-1(h,M)s ^m-1+a _m-2(h,M)s ^m-2+…+a ₁(h,M)s+a ₀(h,M)、

B (h, M, s)=s ⁿ+ b _n-1(h, M) s ^n-1+ b _n-2(h, M) s ^n-2+ ... + b ₁(h, M) s+b ₀(h, M) is polynomial expression, and s is the variable after laplace transform conventional in transport function, h, and M is respectively flying height and Mach number, and σ (h, M) is the time delay of pitch channel, and K (h, M) is with h, the gain that M changes, a _l(h, M), l=0,1,2 ..., m-1 be in polynomial expression A (h, M, s) with h, M change coefficient bunch, b _i(h, M), i=0,1,2 ..., n-1 be in polynomial expression B (h, M, s) with h, M change coefficient bunch, the indeterminate that △ k (s) is model;

(2) transport function of the multistage PID controller of candidate is:

$G_{c,I_E}\left(s\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Π}}}[K_p\left(i\right)+k_I\left(i\right)/s+k_D\left(i\right)\·s]}{\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1}$

In formula, k _cfor constant gain to be determined, N is integer, represents the progression of multistage PID controller to be determined, k _p(i), k _i(i), k _d(i) i=1,2 ..., N is constant to be determined, k _aSEfor Flutter Suppression gain;

Add after multistage PID controller, the open-loop transfer function of whole system is:

$G_{I_E}\left(s\right)G_{c,I_E}\left(s\right)=\frac{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)}{B(h,M,s)}[1+\mathrm{\Δk}\left(s\right)]\frac{\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_p\left(i\right){+k}_I\left(i\right)/s+k_D\left(i\right)\·s]}{\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1}$

Corresponding root locus equation is:

$e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Λk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_p\left(i\right)\·s+k_1\left(i\right)+k_D\·s^2]+s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1)=0;$ s?B(h,M,s)?(kASEs+1)=0

(3) establish s=σ+j ω, wherein: the real part that σ is s, the imaginary part that ω is s, j is the imaginary part of symbol; The stability margin index of system is set as: σ≤-

,

wherein,

for non-zero real, ξ gives fixed number; Set up according to flight test or wind tunnel test the lagging phase angle that model indeterminate causes

radian, amplitude

the stability margin index of system is adjusted into:

with

wherein, △ _mand △ _abe whole real number;

Like this, the stability margin index of system can be converted into: according to

$\left\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)\right[1+\mathrm{\Δk}\left(s\right)\left]\underset{i=1}{\overset{N}{\mathrm{\Π}}}\right[k_p\left(i\right)\·s+k_I\left(i\right)+k_D\left(i\right)\·s^2]+s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}=0$

Or $\left\{\begin{array}{c}\mathrm{Re}\left\{\right\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Δk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_P\left(i\right)\·s+k_I\left(i\right){+k}_D\left(i\right)\·s^2]+\\ s\·B(h,M,s)\·(\frac{k_{\mathrm{SAE}}}{{\mathrm{\ω}}_{\mathrm{SAE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}\\ \mathrm{Im}\left\{\right\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Δk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_P\left(i\right)\·s+k_I\left(i\right)+k_I\left(i\right)+k_D\left(i\right)\·s^2]+\\ s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}\}=0\end{array}\right.$ The root locus obtaining must meet with

according under this index and maximum likelihood criterion or the common constraint of other criterion, can determine according to the maximum likelihood method in system model Structure Identification or discrimination method progression N, the constant k of multistage PID controller _p(i), k _i(i), k _d(i) i=1,2, ", N and Flutter Suppression gain k _aSE.

Claims (1)

1. the multistage PID robust Controller Design of the compound root locus of a multiloop dummy vehicle bunch Flutter Suppression method, is characterized in comprising the following steps:

Under step 1, given differing heights, Mach number by frequency sweep flight test directly by allowing amplitude-frequency and phase-frequency characteristic in the full envelope curve of flight to form primary control surface in the full envelope curve of aircraft and the model cluster of flying height, and can cross over flight envelope and obtain the flutter frequency of aircraft, obtain open-loop transfer function model cluster matrix between corresponding aircraft primary control surface and flying height and be:

Wherein, G is m × m square formation, and m>1 is positive integer, the independent variable that s is Laplace transformation, h is aircraft altitude, and M is Mach number, and Δ is uncertain vector, and P is m × m single mode square formation, D is m × m polynomial expression diagonal matrix, and Q is m × m single mode square formation for polynomial expression, n>1 is positive integer;

Choose

satisfy condition:

and

Wherein, G _efor m × m square formation, P _efor m × m single mode square formation, D _efor m × m polynomial expression diagonal matrix, d _i,Efor D _ei _erow, I _ecolumn element,

for the I of D _erow, I _ecolumn element, I _e=1,2 ..., m, Q _efor m × m single mode square formation, for polynomial expression, arg is phase angle mathematic sign;

The controller of aircraft multiloop system is made as:

Wherein, G _cA(s) be m × m square formation, G _a0(s)=diag[G _{c, 1}(s), G _{c, 2}(s) ..., G _c,m(s)] be m × m diagonal matrix;

(s) be G _a0(s) I _erow, I _ecolumn element, I _e=1,2 ..., m;

Step 2, controller

(s), I _e=1,2 ..., the design process of m is as follows:

(1) order

(s)= (s, h, M, Δ)/

(s, h, M, Δ), the form of embodying is:

$G_{I_E}\left(s\right)=\frac{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)}{B(h,M,s)}[1+\mathrm{\Δk}\left(s\right)]$

Flutter frequency is: ω _aSE(h, M);

Wherein

A(h,M,s)=s ^m+a _m-1(h,M)s ^m-1+a _m-2(h,M)s ^m-2+…+a ₁(h,M)s+a ₀(h,M)、

B (h, M, s)=s ⁿ+ b _n-1(h, M) s ^n-1+ b _n-2(h, M) s ^n-2+ ... + b ₁(h, M) s+b ₀(h, M) is polynomial expression, and s is the variable after laplace transform conventional in transport function, h, and M is respectively flying height and Mach number, and σ (h, M) is the time delay of pitch channel, and K (h, M) is with h, the gain that M changes, a _l(h, M), l=0,1,2 ..., m-1 be in polynomial expression A (h, M, s) with h, M change coefficient bunch, b _i(h, M), i=0,1,2 ..., n-1 be in polynomial expression B (h, M, s) with h, M change coefficient bunch, the indeterminate that △ k (s) is model;

(2) transport function of the multistage PID controller of candidate is:

$G_{c,I_E}\left(s\right)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Π}}}[K_p\left(i\right)+k_I\left(i\right)/s+k_D\left(i\right)\·s]}{\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1}$

In formula, k _cfor constant gain to be determined, N is integer, represents the progression of multistage PID controller to be determined, k _p(i), k _i(i), k _d(i) i=1,2 ..., N is constant to be determined, k _aSEfor Flutter Suppression gain;

Add after multistage PID controller, the open-loop transfer function of whole system is:

$G_{I_E}\left(s\right)G_{c,I_E}\left(s\right)=\frac{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)}{B(h,M,s)}[1+\mathrm{\Δk}\left(s\right)]\frac{\underset{i=1}{\overset{N}{\∏}}[k_p\left(i\right){+k}_I\left(i\right)/s+k_D\left(i\right)\·s]}{\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1}$

Corresponding root locus equation is:

$e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Λk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_p\left(i\right)\·s+k_1\left(i\right)+k_D\·s^2]+s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1)=0;$

(3) establish s=σ+j ω, wherein: the real part that σ is s, the imaginary part that ω is s, j is the imaginary part of symbol; The stability margin index of system is set as: σ≤-

,

wherein, for non-zero real, ξ gives fixed number; Set up according to flight test or wind tunnel test the lagging phase angle that model indeterminate causes radian, amplitude

the stability margin index of system is adjusted into:

with wherein, △ _mand △ _abe whole real number;

Like this, the stability margin index of system can be converted into: according to

$\left\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)\right[1+\mathrm{\Δk}\left(s\right)\left]\underset{i=1}{\overset{N}{\mathrm{\Π}}}\right[k_p\left(i\right)\·s+k_I\left(i\right)+k_D\left(i\right)\·s^2]+s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}=0$

Or $\left\{\begin{array}{c}\mathrm{Re}\left\{\right\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Δk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_P\left(i\right)\·s+k_I\left(i\right){+k}_D\left(i\right)\·s^2]+\\ s\·B(h,M,s)\·(\frac{k_{\mathrm{SAE}}}{{\mathrm{\ω}}_{\mathrm{SAE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}\\ \mathrm{Im}\left\{\right\{e^{-\mathrm{\σ}(h,M)s}K(h,M)A(h,M,s)[1+\mathrm{\Δk}\left(s\right)]\underset{i=1}{\overset{N}{\mathrm{\Π}}}[k_P\left(i\right)\·s+k_I\left(i\right)+k_I\left(i\right)+k_D\left(i\right)\·s^2]+\\ s\·B(h,M,s)\·(\frac{k_{\mathrm{ASE}}}{{\mathrm{\ω}}_{\mathrm{ASE}}}s+1){\}}_{s=\mathrm{\σ}+\mathrm{j\ω}}\}=0\end{array}\right.$ The root locus obtaining must meet

with

according under this index and maximum likelihood criterion or the common constraint of other criterion, can determine according to the maximum likelihood method in system model Structure Identification or discrimination method progression N, the constant k of multistage PID controller _p(i), k _i(i), k _d(i) i=1,2 ..., N and Flutter Suppression gain k _aSE.