CN116165896A - An Aircraft Adaptive Control Method Based on Online Frequency Domain Recursive Identification - Google Patents

Abstract

The invention belongs to the technical field of aircraft control, and relates to an aircraft adaptive control method based on online frequency-domain recursive identification, including an online recursive frequency-domain aerodynamic parameter identification method and an improved adaptive interference suppression control method. The aerodynamic derivative can be accurately identified by using the recursive identification method, and the identification error of static stability coefficient and rudder effect coefficient does not exceed 1%. The improved adaptive dynamic inverse control introduces the advanced correction link on the basis of the traditional adaptive dynamic inverse control, so that the adaptive dynamic inverse controller responds faster. Combined with the online identification method, it can realize high-precision control under the condition of inaccurate aerodynamic model. Applied to aircraft pitch attitude control, compared with traditional controllers and ordinary adaptive dynamic inverse controllers, the overshoot is smaller, the rise time and peak time are shorter, and it has good control quality.

Description

An adaptive control method for aircraft based on online frequency domain recursive identification

Technical Field

The invention belongs to the technical field of aircraft control and relates to an aircraft adaptive control method based on online frequency domain recursive identification.

Background Art

The control systems used by the vast majority of aircraft in service today are traditional controllers, such as PID controllers. However, due to the nonlinear characteristics of real aircraft aerodynamic models, the control performance of traditional methods cannot fully meet expectations and has poor robustness. It also brings secondary problems such as a long control system development cycle. In order to overcome this problem, in recent years, a large amount of research has focused on nonlinear control system design, such as adaptive control methods, intelligent control methods, etc.

Among the published studies, most of them studied aerodynamic identification methods and adaptive control methods separately, and achieved fruitful research results in both aspects.

The role of aerodynamic identification includes the correction of aerodynamic models, control reconstruction in fault conditions, real-time adjustment of control parameters, etc. In some application scenarios, the aerodynamic model must meet the requirements of real-time identification. As the accumulated flight data increases, accurate identification is gradually achieved. Accurately obtaining the real-time dynamic information of the aircraft is the prerequisite for nonlinear adaptive control adjustment. The method used here is dynamic parameter identification. The methods of aircraft dynamic identification include two categories: time domain identification and frequency domain identification. Specifically, they include regression methods represented by the least squares method. At present, many methods have been proposed based on the least squares method. In addition, commonly used identification methods also include the maximum likelihood method and the Kalman filter method. In terms of frequency domain methods, some practical methods have also been proposed. In addition, intelligent methods are gradually being used in aerodynamic identification, among which the more representative one is the use of neural networks for dynamic identification.

In terms of nonlinear adaptive control methods, there are many different technical routes, such as a series control system that combines a neural network or fuzzy controller with a PID controller to achieve real-time adjustment of the control gain based on the control error or other observables. Adaptive dynamic inversion (NDI) control is also an adaptive control method. The dynamic inversion method is a control method with strong adaptability and versatility. Since this method can better cope with the complicated decoupling work between the channels of nonlinear objects and there is no need to repeatedly adjust the gains of each loop, it is also widely used in aircraft nonlinear control. In the long development process, dynamic inversion control has produced many different forms by adjusting the controller structure or combining with other controllers. Adaptive dynamic inversion based on L1 is a representative direction.

Combining online aerodynamic identification with adaptive control methods to achieve online regulation is a learning control method route proposed in recent years. Some research projects have also conducted considerable technical exploration work on this control optimization idea. However, in these studies, there are few results in the online recursive application research in the frequency domain, and even fewer people have been involved in the method of combining frequency domain identification with improved adaptive dynamic inverse control methods.

Summary of the invention

In order to solve the above problems, the present invention provides an aircraft adaptive control method based on online frequency domain recursive identification.

The technical solution of the present invention:

An aircraft adaptive control method based on online frequency domain recursive identification includes an online recursive frequency domain aerodynamic parameter identification method and an improved adaptive interference suppression control method. The details are as follows:

(1) Aircraft six-degree-of-freedom dynamics modeling

The aircraft center of mass dynamics equation can be expressed as:

Where m is the mass of the aircraft, V = [V _x V _y V _z ] represents the component of the aircraft velocity in the body coordinate system, Ω = [pq r] represents the component of the aircraft angular velocity in the body coordinate system, and [F _x F _y F _z ] represents the component of the total external force on the aircraft in the body coordinate system. The dot above a variable represents the derivative of the variable, and the same applies below.

Combined with the aircraft force analysis, the aircraft's center of mass dynamics equation is finally expressed as follows:

In the formula, G is the gravity acting on the aircraft, D is the drag, L is the lift, F _Ay is the lateral force, φ _T is the engine installation angle, and F _Ty is the component of the engine thrust on the O _b y _b axis.

The aircraft inertia matrix I is:

In the inertia matrix, _Ixx , _Iyy , and _Izz represent the moments of inertia _of the aircraft around the three axes of the fuselage _Obxb , _{Obyb ,} and _{Obzb ,} respectively, and _Ixy , _Ixz , ..., _Iyz are the products of inertia of the aircraft around each axis of the fuselage coordinate system.

The equation for angular momentum H is:

Where H _x , _Hy , and H _z are the components of angular momentum in the body coordinate system.

The motion of an aircraft around its center of mass can be expressed as:

Where _LA , _MA and _NA are the components of aerodynamic torque in the body coordinate system, _LT , _MT and _NT are the components of torque generated by thrust in the body coordinate system, and gravity does not generate torque because it passes through the center of mass.

(2) Frequency domain online recursive identification method

The frequency domain online recursive identification includes three parts: multi-sinusoidal control surface excitation, data detrending and CZT transformation. The specific flow chart is shown in Figure 3. After the aircraft reaches the predetermined flight state, the excitation signal is applied to the control surface, and the aerodynamic force or torque coefficient is calculated after the sensor data is obtained. The acquired data and the calculated aerodynamic coefficient are detrended, and the data that is linearly related to time is eliminated, and the dynamic component is retained. Then, the data is converted from the time domain to the frequency domain band of interest through CZT transformation, and recursive least squares identification is performed. After obtaining the identification result, the static information is restored to obtain the final identification result, and it is judged whether the identification converges.

This streamlined online recursive identification method is suitable for force or torque identification.

The following is a detailed introduction to the contents of the three parts.

(2.1) Multi-sine rudder excitation:

When the aircraft is in a stable flight, the flight data contains very little information. This low-information flight data cannot support accurate identification of the aircraft's aerodynamic coefficients. Therefore, an effective maneuvering method must be adopted to increase the amount of information contained in the flight data to meet the identification requirements.

The method adopted by the present invention is a multi-sine superposition input rudder excitation signal, which can be superimposed with the control command as the final control rudder deflection command. The multi-sine excitation input is in the form of:

Here, the disturbance input applied to the jth control surface is assumed to be u _j , which is the sum of sinusoidal harmonics φ _k with a single phase shift. Where m is the total number of available harmonically related frequencies, T is the time length of the excitation. _{A k} is the amplitude of the kth sinusoidal component, and t is the time vector.

(2.2) Data detrending:

In order to meet the real-time control of the aircraft, a real-time detrending method must be adopted. Here, a fourth-order Butterworth high-pass filter is used for detrending.

A recursive method is used here for filtering, which is a critical point because this method can realize real-time detrending of sensor data and ultimately achieve the effect of real-time identification.

The filter expression can be written in the form of a difference equation of the following rational transfer function:

a(1)y(n)=b(1)x(n)+b(2)x(n-1)+…+b(n _b +1)x(nn _b )

-a(2)y(n-1)-…-a(n _a +1)y(nn _a )

Where n _a is the order of the feedback filter, n _b is the order of the feedforward filter, a(i) represents the denominator coefficient of the rational transfer function vector, b(i) represents the numerator coefficient of the rational transfer function vector, x(i) is the input data, and y(i) is the filtered data.

(2.3) CZT transformation:

When converting time domain data to frequency domain, the most commonly used method is the DFT method, but the DFT method has many shortcomings, especially for aircraft dynamics identification, usually only a small frequency band needs to be concerned, and DFT will introduce interference information. Therefore, the CZT method is used here. The data frequency domain transformation effect of this method is shown in Figure 4.

The starting point z ₀ of the CZT transformation can be selected arbitrarily, so the input data can be analyzed with narrowband high resolution starting from any frequency. At the same time, the frequency resolution can be adjusted as needed, so that any frequency band of interest can be precisely analyzed.

The steps of CZT transformation can be expressed as:

First, form two L-point sequences g(n) and h(n):

Wherein, L satisfies L≥N+M-1, and L＝ ^2M , M is an integer, and A and W are expressed as:

表示相邻两点等分角，W₀表示螺旋线伸展率。A ₀ represents the radius of the initial sampling point, θ ₀ represents the phase angle of the initial sampling point,

W 0 represents the angle divided equally between two adjacent points, and W ₀ represents the helix stretch rate.

Then, the Fourier transforms G(k) and H(k) of g(n) and h(n) are calculated respectively, and Y(k)=G(k)H(k) is calculated.

Then calculate the L-point inverse Fourier transform X(z _k ) of Y(k):

To achieve the effect of real-time identification in the conventional CZT transformation process, it is necessary to use recursive finite Fourier transform at the same time. The next moment _Xi (ω) can be expressed as:

X _i (ω)=X _i-1 (ω)+x(i)e ^-jωiΔt

in

e ^-jωiΔt =e ^-jωΔt e ^-jω(i-1)Δt

(3) Adaptive interference suppression control

(3.1) Adaptive disturbance rejection controller:

In order to track the command of each variable, the time scale separation of the variables is used to sequentially use nonlinear dynamic inversion (NDI) to generate commands for faster variables. The controller structure is shown in Figure 5. The dynamic inversion of the outermost fast loop converts the commands of the trajectory inclination angle χ and the trajectory deflection angle γ into θ and φ commands. For linear tracking, the required control derivative is proportional to the error between the variable and its command. According to the dynamic equations, the roll angle command generated by the inverse process is:

Where, _χc and _γc are the instructions for the trajectory inclination and trajectory deflection, _θc and _φc are the instructions for θ and φ, the lower right subscript of the variable represents the change instruction value, the same below, V is the total velocity of gravity, g is the local acceleration of gravity, and _Kχ is the outer loop control gain.

The trajectory inclination derivative is expressed as:

The above two equations can be simplified to:

According to the first-order linear non-homogeneous differential equation solution formula, we can get:

Within a certain period of time, χ converges to the trajectory inclination angle command χ _c , where χ(0) is the initial value of the trajectory inclination angle.

θ、β的指令求解p、q、r的指令，根据动力学方程：则有如下指令形式：The dynamic inversion of the inner slow loop is to solve the next step of dynamic inversion instructions.

The instructions of θ and β solve the instructions of p, q and r. According to the dynamic equation: the instructions are in the following form:

Among them, K _φ , K _β , K _θ are the angular velocity loop control gains.

Bringing the above control instructions back to the dynamic equation, we can get:

According to the first-order linear non-homogeneous differential equation solution formula, we can get:

The time domain response function with the same form as the track angle can be obtained, and it can converge within a certain period of time.

Solving the inverse of the dynamic equations for the inner loop that tracks these angular rate commands, yielding the required angular acceleration proportional to the angular rate error, yields:

是由飞机当前的辨识模型确定，即通过步骤(2)的方法获得的当前辨识结果和传感器指令反向求解的俯仰力矩；K_q为角速度控制增益，M_δ,c为控制舵偏力矩指令。In the formula,

It is determined by the current identification model of the aircraft, that is, the pitch moment obtained by inversely solving the current identification result and the sensor command obtained by the method of step (2); K _q is the angular velocity control gain, and M _δ,c is the control rudder deflection moment command.

的动力学方程中可得：Substitute the control torque command into

From the kinetic equation, we can get:

Where M is the torque value obtained by the sensor.

(3.2) Improved adaptive disturbance rejection controller:

In order to solve the problem of large error and strong hysteresis based on conventional dynamic inversion under the action of unsteady aerodynamic force, an improved adaptive dynamic inversion control method is introduced to improve the control effect when there are errors in the aerodynamic identification results. The controller structure is shown in Figure 6. The details are as follows:

Considering that the change of airflow angle [α β μ] is slower than the change of angular velocity [p q r], where μ is the track roll angle, the system is divided into a slow-changing subsystem and a fast-changing subsystem. In the slow loop:

为慢回路中右下角标相应变量的期望动态响应。In the formula, ω _αd , ω _βd , ω _μd are the bandwidth frequencies corresponding to [α β μ],

is the expected dynamic response of the corresponding variable in the lower right corner in the slow loop.

Due to the uncertainty of the aerodynamic model, the conventional dynamic inversion method has a slow response speed, which leads to command tracking errors. This situation can be improved by connecting the advance correction link in series. Then the slow loop expectation system can be written as:

Among them, G(s) is the advance correction device:

ω _C is the position of the pole in the complex plane, ω _D is the position of the zero point, and the variable name marked in the lower right corner of the zero pole represents the corresponding variable. When selecting parameters, when ω _C <ω _D , G(s) is the leading link, and when ω _C >ω _D , G(s) is the lagging link, so here we should take ω _C <ω _D.

After the lead correction link is introduced, it needs to be implemented in discrete control simulation, so the transfer function of the lead link is discretized.

The advance link is added when solving the nominal angular rate in the slow loop. The conventional slow loop angular rate solution is expressed as:

q _c =K _θ (θ _c -θ)

After adding the lookahead phase, the expression becomes:

为改进的期望角速度。in,

is the desired angular velocity for improvement.

Written in time domain:

Using approximate differences instead of differentials, we can get the nominal angular rate after adding the lead link:

Beneficial effects of the present invention:

The present invention focuses on the study of an adaptive control method based on real-time frequency domain recursive identification, focusing on the feasibility of its real-time application and the effect of the control method. The frequency domain aerodynamic parameter identification method is insensitive to noise, which makes it of great value in practical applications. The introduced recursive detrending method further reduces the computational requirements. The aerodynamic derivatives can be accurately identified using the recursive identification method, and the identification errors of the static stability coefficient and the rudder efficiency coefficient do not exceed 1%. The improved dynamic inverse control introduces an advance correction link on the basis of the traditional dynamic inverse control, so that the dynamic inverse controller responds faster. Combined with the online identification method, high-precision control can be achieved when the aerodynamic model is inaccurate. It is applied to aircraft pitch attitude control. Compared with traditional controllers and ordinary adaptive dynamic inverse controllers, the overshoot is smaller, the rise time and peak time are shorter, and it has good control quality.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a flow chart of an improved adaptive control method based on online frequency domain aerodynamic identification;

FIG2 is a flow chart of an improved adaptive control algorithm based on online frequency domain aerodynamic identification;

FIG3 is a flow chart of a method for online recursive aerodynamic parameter identification in the frequency domain;

FIG4 is a schematic diagram of a linear frequency modulation z-transform function;

FIG5 is a flow chart of an adaptive dynamic inverse control method;

FIG6 is a flow chart of an improved adaptive dynamic inverse control method;

FIG7 is a diagram of the rudder excitation signal and its components;

FIG8 is a diagram of the flight status of the identification time interval;

Figure 9 is the result of detrending the flight data in the identification interval;

Figure 10 is a comparison diagram of the static stability coefficient identification result and the true value;

Figure 11 is a comparison chart of the rudder efficiency coefficient identification result and the true value;

Figure 12 is a comparison of the simulation tracking effects of PID and dynamic inverse step signals;

Figure 13 is a comparison of the rudder angle tracking effects of PID and dynamic inverse step signal simulation;

FIG14 is a comparison of the effects of PID and dynamic inverse step signal simulation tracking pitch angular velocity;

Figure 15 is a comparison of the simulation tracking effects of PID, dynamic inversion and improved dynamic inversion step signals;

Figure 16 is a comparison of the rudder angle tracking effects of PID, dynamic inversion and improved dynamic inversion step signal simulation;

FIG17 is a comparison of the effects of PID, dynamic inverse and improved dynamic inverse step signal simulation tracking pitch angular velocity.

DETAILED DESCRIPTION

The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.

The improved adaptive control method based on online frequency domain aerodynamic identification of the present invention is shown in FIG1 and FIG2 .

In the simulation example of the method implementation in this part, it generally includes two major parts, namely the frequency domain online recursive identification part and the adaptive control part based on the online identification results. For real-time identification in the frequency domain, real-time flight data is obtained through real-time flight simulation. First, the real-time data is recursively detrended, and the results of detrending are analyzed. Then, it is substituted into the frequency domain recursive identification program to observe the accuracy of the identification results and determine the identification accuracy. In the control simulation part, the control effects of the conventional PID controller and the adaptive dynamic inverse controller are compared, and then the improved NDI controller with lead correction is introduced to compare the control effects of the three, and the robustness of the improved adaptive dynamic inverse controller is verified. Finally, a full-process simulation including real-time frequency domain identification and improved adaptive dynamic inverse control is given to prove that the method of the present invention can meet the requirements of real-time operation.

(1) Input the initial state and give the target state

Table 1 Simulation initial parameters

The target state is to maintain level flight at a pitch angle of 11°.

(2) Establishing a flight dynamics model

1) Ground coordinate system

Take a point on the ground (the initial position of the aircraft) as the origin O _g , the O _g x _g axis is in the horizontal plane along the take-off direction, forward is positive, the O _g y _g axis is in the horizontal plane and perpendicular to the O _g x _g axis, right is positive, and the O _g z _g axis is perpendicular to the ground, downward is positive.

2) Body coordinate system

Take the center of mass of the aircraft as the origin O _b of the fuselage coordinate system. The coordinate system is fixed to the aircraft. The O _b x _b axis coincides with the axis of the fuselage and is in the longitudinal symmetry plane of the fuselage. The forward direction is positive. O _b y _b is perpendicular to the longitudinal symmetry plane and is positive to the right. O _b z _b is located in the longitudinal symmetry plane and is perpendicular to O _b x _b and is positive downward.

3) Velocity coordinate system/airflow coordinate system

Take the center of mass of the aircraft as the coordinate origin _Oa . _Oaxa is consistent with the direction of the aircraft's airspeed, and is positive forward. _{The Oaza axis is} vertical and located in the longitudinal symmetry plane, and is positive downward. _{Oaya is} perpendicular to the _{plane Oaxaza} , and is positive to the right.

4) Track coordinate system

The track coordinate system is fixed to the aircraft, with the aircraft's center of mass as the origin _Os , _Osxs pointing to the projection direction of the aircraft _'s velocity in the aircraft's longitudinal symmetry plane, _{Oszs axis} in the aircraft's longitudinal symmetry plane, perpendicular to the axis and pointing down the aircraft, and _{Osys axis} perpendicular to the plane and pointing to the right of the aircraft.

The present invention uses several coordinate transformation relationships, which are:

1) Transformation matrix from ground coordinate system to body coordinate system

2) Transformation matrix from airflow coordinate system to body coordinate system

Where θ, φ, ψ represent the pitch angle, roll angle, and yaw angle respectively, α represents the angle of attack, and β represents the sideslip angle.

(1.2) Aircraft six-degree-of-freedom dynamics modeling

The aircraft center of mass dynamics equation can be expressed as:

Where m is the mass of the aircraft, V = [V _x V _y V _z ] represents the component of the aircraft velocity in the body coordinate system, Ω = [pq r] represents the component of the aircraft angular velocity in the body coordinate system, and [F _x F _y F _z ] represents the component of the total external force on the aircraft in the body coordinate system. The dot above a variable represents the derivative of the variable, and the same applies below.

Combined with the aircraft force analysis, the aircraft's center of mass dynamics equation is finally expressed as follows:

In the formula, G is the gravity acting on the aircraft, D is the drag, L is the lift, F _Ay is the lateral force, φ _T is the engine installation angle, and F _Ty is the component of the engine thrust on the O _b y _b axis.

The aircraft inertia matrix I is:

In the inertia matrix, the three items _Ixx , _Iyy , and _Izz represent the moments of inertia of the aircraft around the _three axes of the fuselage _{Obxb , Obyb ,} and _{Obzb ,} respectively, and _Ixy , _Ixz ,... _Iyz are the products of inertia of the aircraft around the axes of the fuselage coordinate system.

The equation for angular momentum H is:

Wherein, H _x , _Hy , and H _z are the components of angular momentum in the body coordinate system.

The motion of an aircraft around its center of mass can be expressed as:

Where _LA , _MA and _NA are the components of aerodynamic torque in the body coordinate system, _LT , _MT and _NT are the components of torque generated by thrust in the body coordinate system, and gravity does not generate torque because it passes through the center of mass.

(3) Frequency domain online identification process simulation example

This embodiment uses the pitch moment coefficient of an aircraft as an example for simulation. The pitch moment coefficient of an aircraft can be simplified into the following form by a linearization method. Here, we mainly focus on the moment, which can be regarded as being related to some flight state quantities and satisfying an instantaneous linear relationship:

Where C _l , C _m , and C _n are the pitch, yaw, and roll moment coefficients, respectively. The lower right subscripts of the three variables indicate that a certain variable represents the derivative of the moment coefficient with respect to the variable. The lower right subscript 0 represents the moment coefficient caused by the body configuration, b represents the wingspan, and c represents the average aerodynamic chord length.

Sometimes the modeling terms on the right side of the equation also include some higher-order terms or coupling terms of different variables. Appropriate adjustments need to be made for aircraft of different configurations to accurately model the aircraft, and a trade-off should be made between the number of modeling variables and the amount of dynamic response information, so as to achieve the purpose of accurate identification while minimizing the amount of calculation.

Here we mainly take the pitch moment as an example for analysis, and define the variable X as a vector composed of data directly obtained from the aircraft sensor, which contains the following content:

The detrending method introduced above can be used to achieve recursive detrending of data. The vector after detrending is defined as X _d , and its content is expressed as follows:

其表达式为：The detrended data can be recursively Fourier transformed using the method described above to obtain

Its expression is:

In the frequency domain, the recursive solution problem can be expressed in the form of recursive least squares, as follows:

The numerical value of the identification result can be obtained by the following expression:

Finally, it is necessary to restore the static information filtered out in the detrending process:

Taking longitudinal attitude control as an example, after the identification is completed, the control surface command can be solved based on the identification results. The formula is:

Applying active excitation signals to the control surfaces can stimulate the dynamic characteristics of the aircraft, so that the flight data contains more information and the accuracy of identification is improved. Since the present invention only focuses on the longitudinal pitch force identification, the excitation is only applied to the elevator. The excitation expression is:

Where T is the excitation period, and the corresponding image is shown in Figure 7:

The flight state of the identification phase selected by the present invention is shown in FIG8 :

The flight altitude is about 5000m, the speed is about 150m/s, and dynamic identification is performed in horizontal flight state.

A 4th-order Butterworth high-pass filter is used, and the filter-related information is expressed as follows:

Table 2 Filter status

First, the data needs to be detrended. In order to achieve real-time identification, a recursive detrending method is used here. The comparison of the data before and after processing is shown in Figure 9.

In the figure, the solid line represents the unprocessed raw data. When the data has a constant component or a component that is linearly related to time, it will not be able to be identified. In the above figure, it can be clearly seen that the angle of attack and the elevator angle have obvious constant components, so these data must be processed. Relatively speaking, the other two items have smaller deviations, but they also need to be processed.

The dotted line in the figure is the result of real-time recursive detrending. In order to observe the detrending effect more intuitively, after the simulation, a more accurate batch detrending method is used offline to process the data, and all the results are displayed in a picture. It must be pointed out that this batch processing method does not meet the requirements of real-time operation, because it requires all flight data to be processed simultaneously. If all previous data are batch processed after each data acquisition, the effect similar to the recursive method can be achieved, but this will cause a very large amount of calculation and put forward high requirements on the hardware computing power of the aircraft. The result of offline batch processing is represented by the blue dotted line. It can be clearly seen from the figure that the recursive method needs about 20 seconds to converge to the same accuracy as the batch processing method. This is because at the initial moment, there are very few flight data and it is impossible to accurately identify the constant components contained in the flight data.

It can be seen from this embodiment that when the recursive real-time detrending method is used, it is necessary to wait for a period of time after starting the detrending process before carrying out the recursive aerodynamic identification. The specific waiting time length is determined by flight data analysis. The waiting time used here is 20 seconds.

After detrending, recursive identification in the frequency domain can be performed. The identification results are shown in Figures 10 and 11:

While waiting for the real-time filter to converge, the identification result is manually set to 0. In order to ensure that the identification result does not fluctuate greatly at the beginning of the recursion, the first 90 samples are taken for batch processing, and the result is used as the initial value for subsequent recursion. Similarly, while waiting for the data collection stage of batch processing, the identification result is manually set to 0.

In the figure, we can see that the error of the initial value obtained after batch processing is already very small, and there is no deviation in the subsequent recursive process, which shows that the identification is fast, the convergence is good, and the error is small. In order to more clearly illustrate the identification accuracy, a period of time after the start of recursion is selected to average the relative error of parameter identification, and the calculation results are shown in the following table:

Table 3 Relative error of aerodynamic parameter identification

(4) Adaptive dynamic inversion control simulation based on real-time aerodynamic parameter identification

First, the conventional configuration adaptive dynamic inverse controller is simulated and verified, and the control effect comparison between the adaptive dynamic inverse controller and the PID controller under the step command is explored. Here, the aerodynamic identification link in the previous part is first carried out, including the recursive detrending and real-time recursive identification. The control method is converted only after the identification is completed. Finally, by giving a step pitch angle command, the step pitch angle command is tracked and observed in the PID and adaptive dynamic inverse control modes respectively. The simulation of this part is carried out according to this process, but since the accuracy of the identification has been explored and explained in detail in the previous part, this chapter only shows the image of the command tracking part.

It can be clearly seen from simulation figures 12, 13, and 14 that the adaptive dynamic inverse controller has a faster response speed and a smaller overshoot than the PID controller. From the elevator deflection angle image and the pitch angle velocity image, it can be seen that the adaptive dynamic inverse controller can simultaneously achieve a larger angle of rudder deflection and a faster adjustment response when controlling the aircraft to track the step pitch angle command, thereby achieving a better control effect. In order to compare the difference between the two numerically, the three indicators of rise time, peak time, and overshoot are selected to quantitatively compare the control effects of the two.

Table 4 Comparison of control effects between PID controller and adaptive dynamic inverse controller

Quantitative analysis shows that the adaptive dynamic inversion controller has better control effects than the PID controller in terms of the three indicators.

(5) Improved adaptive dynamic inversion control simulation based on real-time aerodynamic parameter identification

The simulation verification of the improved adaptive dynamic inverse control method was carried out under the same conditions. Similarly, all the links including recursive detrending, identification and step pitch angle command were adopted, and only the image of the step command tracking stage was displayed.

Figures 15, 16 and 17 show the step response tracking effects of the improved adaptive dynamic inverse control method, conventional adaptive dynamic inverse control and PID control.

It can be seen from the simulation diagram that compared with the PID controller and the adaptive dynamic inverse controller, the improved adaptive dynamic inverse controller has further improved the response speed and overshoot. Similarly, from the elevator deflection angle image and the pitch angular velocity image, it can be seen that the improved adaptive dynamic inverse controller has a better response effect of the elevator deflection angle. In order to compare the difference between the two numerically, the three indicators of rise time, peak time and overshoot are selected to quantitatively compare the control effects of the three.

Table 5 Comparison of control effects of PID controller, adaptive dynamic inverse controller and improved adaptive dynamic inverse controller

Through the above analysis, it can be concluded that the improved adaptive dynamic inverse control of the present invention has a better control effect than the conventional adaptive dynamic inverse control.

Claims (1)

1. The aircraft self-adaptive control method based on the online frequency domain recursion identification is characterized by comprising an online recursion frequency domain pneumatic parameter identification method and an improved self-adaptive interference suppression control method; the method comprises the following steps:

(1) Six-degree-of-freedom dynamics modeling of aircraft

The aircraft centroid dynamics equation is expressed as:

wherein m is aircraft mass, V= [ V ] _x V _y V _z ]Representing the component of the aircraft speed in the body coordinate system, Ω= [ pqr ]]Representing the component of the angular velocity of an aircraft in the body coordinate system, [ F ] _x F _y F _z ]Representing the components of the combined external force applied by the aircraft under the coordinate system of the aircraft body; the variable superscript indicates the derivative of the variable;

in connection with aircraft stress analysis, the mass dynamics equation of the aircraft is finally expressed as follows:

wherein G is the gravity of the aircraft, D is the resistance, L is the lift,

is a lateral force phi _T For engine mounting angle->

For engine thrust at O _b y _b A component of the shaft;

the aircraft inertia matrix I is:

in the inertial matrix, I _xx 、I _yy 、I _zz Respectively represent the axis O of the aircraft around the aircraft body _b x _b 、O _b y _b 、O _b z _b Moment of inertia of three axes, I _xy 、I _xz 、……I _yz The inertia products of the aircraft around each axis of the machine body coordinate system are respectively calculated;

the angular momentum H equation is:

in the formula ,H_x 、H _y 、H _z The components of the angular momentum under the machine body coordinate system are respectively;

the movement of the aircraft around the centroid is expressed as:

in the formula ,L_A 、M _A 、N _A Components of aerodynamic moment under the machine body coordinate system, L _T 、M _T 、N _T The components of the moment generated by the thrust under the machine body coordinate system are respectively, and the gravity is not generated due to the fact that the gravity passes through the center of mass;

(2) Frequency domain online recursion identification method

The frequency domain online recursion identification comprises three parts of contents of multi-sine control surface excitation, data trend removal and CZT conversion; after the aircraft reaches a preset flight state, applying an excitation signal to a control surface, acquiring sensor data, and calculating aerodynamic force or moment coefficient; carrying out data trending on the acquired data and the calculated pneumatic coefficient, removing the data which has a linear relation with time, retaining a dynamic component, then converting the data from a time domain into a frequency domain frequency band of interest through CZT conversion, carrying out recursive least square identification, recovering static information after an identification result is obtained, obtaining a final identification result, and judging whether the identification is converged or not; the three parts are specifically as follows:

(2.1) multiple sinusoidal control surface excitation:

a control surface excitation signal which is input by multi-sine superposition is adopted, and the signal is superposed with a control instruction to be used as a final control surface deflection instruction; the form of the multi-sinusoidal excitation input is:

the disturbance input applied to the j-th control plane is u _j Is a sine wave harmonic phi with a single phase shift _k Sum up; where m is the total number of available harmonically related frequencies and T is the length of time of excitation; a is that _k Is the amplitude of the kth sine wave component, and t is the time vector;

(2.2) data trending:

in order to meet the real-time control of the aircraft, a fourth-order Baud Wo Sigao pass filter is adopted for trending;

the expression of the filter is written as a differential equation of the rational transfer function:

a(1)y(n)＝b(1)x(n)+b(2)x(n-1)+…+b(n _b +1)x(n-n _b )-a(2)y(n-1)-…-a(n _a +1)y(n-n _a )

in the formula ,n_a Is the order of the feedback filter, n _b Is the order of the feedforward filter, a (i) represents the denominator coefficient of the rational transfer function vector, b (i) represents the numerator coefficient of the rational transfer function vector, x (i) is the input data, and y (i) is the filtered data;

(2.3) CZT conversion:

when converting the time domain data into the frequency domain, adopting a CZT method; the CZT transformation steps are as follows:

first two L-point sequences g (n) and h (n) are constructed:

wherein L satisfies L.gtoreq.N+M-1, and L=2 ^M M is an integer, A, W is expressed as:

A ₀ represents the initial sampling point radius, θ ₀ Indicating the phase angle of the initial sampling point,

represents the equal division angle of two adjacent points, W ₀ Representing the helix stretching rate;

then fourier transforms G (k) and H (k) of G (n) and H (n), respectively, are calculated, and Y (k) =g (k) H (k) is calculated;

then calculate the L-point Fourier inverse transform X (z) of Y (k) _k )：

The CZT transformation process needs to adopt recursive finite Fourier transformation at the same time to realize the effect of real-time identification, and the latter moment X _i (ω) is expressed as:

X _i (ω)＝X _i-1 (ω)+x(i)e ^-jωiΔt

wherein

e ^-jωiΔt ＝e ^-jωΔt e ^-jω(i-1)Δt

(3) Adaptive interference suppression control

(3.1) an adaptive interference suppression controller:

to track the commands for each variable, nonlinear dynamic inversion is employed in turn to generate commands for faster variables using time scale separation of the variables; the dynamic inversion of the outermost loop is to convert the instruction of the track inclination angle χ and the track deflection angle γ into an instruction of θ and φ; for linear tracking, the required control derivative is proportional to the error between the variable and its command, and the roll angle command generated by the inverse process is, according to the system of kinetic equations:

in the formula ,χ_c and γ_c Is the instruction of track inclination angle and track deflection angle theta _c and φ_c For the theta and phi instructions, the right lower corner mark of the variable represents the instruction value of the change amount, and the following is that V is the gravity combined speed, g is the local gravity acceleration, K _χ Controlling gain for the outer loop;

the track tilt derivative is expressed as:

the two types of the combination are simplified to obtain:

the method is obtained by solving a formula according to a first-order linear non-homogeneous differential equation:

in a certain time period, χ converges to the track inclination instruction χ _c X (0) is the initial value of the track inclination angle;

the dynamic inversion of the inner loop slow loop is to solve the dynamic inversion instruction of the next step according to

Instructions of theta and beta solve instructions of p, q and r, and according to a dynamics equation: then there is the following instruction form:

wherein ,K_φ 、K _β 、K _θ Gain is controlled for the angular velocity loop;

and (3) bringing the control command back to a dynamics equation to obtain:

the method is obtained by solving a formula according to a first-order linear non-homogeneous differential equation:

obtaining a time domain response function which is the same as the track angular form and can achieve convergence in a certain time;

solving the inverse of the kinetic equation that tracks these angular rate command inner loops, produces the required angular acceleration proportional to the angular rate error, yielding:

in the formula ,

is made of flyingDetermining a current identification model of the machine, namely determining a current identification result obtained by the method of the step (2) and a pitching moment reversely solved by a sensor instruction; k (K) _q For angular velocity control gain, M _δ,c A rudder deflection moment command is controlled;

substituting control moment command into

Can be obtained in the kinetic equation of (a):

wherein M is a moment value obtained by a sensor;

(3.2) improving the adaptive interference suppression controller:

an improved self-adaptive dynamic inverse control method is introduced to improve the control effect under the condition that the pneumatic identification result has errors; the method comprises the following steps:

considering that the change in the air flow angle [ αβμ ] is slow relative to the change in the angular velocity [ pqr ], where μ is the track roll angle, the system is thus divided into a slow-varying subsystem and a fast-varying subsystem, in a slow loop:

in the formula ,ω_αd 、ω _βd 、ω _μd Respectively [ alpha beta mu ]]The frequency of the corresponding bandwidth is chosen to be the frequency of the corresponding bandwidth,

expected dynamic response for the corresponding variable for the lower right corner mark in the slow loop;

because of uncertainty of the pneumatic model, the response speed of the conventional dynamic inverse method is low, so that an error of command tracking is brought, and the situation is improved through a series lead correction link, so that a slow loop expects the system to write as follows:

wherein G(s) is an advance correction device:

ω _C omega is the complex plane pole position _D The zero position is the zero position, and the name of the variable is identified by the right lower angle of the zero pole and corresponds to the variable; when selecting parameters, ω _C <ω _D In the case of G(s) being the leading element, when ω _C >ω _D G(s) is a hysteresis, and therefore ω is taken _C <ω _D ；

After the lead correction link is introduced, the lead correction link is required to be realized in discrete control simulation, so that the transfer function of the lead link is subjected to relevant discretization;

the lead link is added when the nominal angular rate is solved in the slow loop; the slow loop angular rate solver is expressed as:

q _c ＝K _θ (θ _c -θ)

after the advance link is added, the expression becomes:

wherein ,

for improved desired angular velocity;

written in the form of the time domain:

the approximate difference is used for replacing the differential, and the nominal angular rate after the advance link is added is obtained by arrangement:

Images