CN116165896A - Airplane self-adaptive control method based on online frequency domain recursion identification - Google Patents
Airplane self-adaptive control method based on online frequency domain recursion identification Download PDFInfo
- Publication number
- CN116165896A CN116165896A CN202310168326.4A CN202310168326A CN116165896A CN 116165896 A CN116165896 A CN 116165896A CN 202310168326 A CN202310168326 A CN 202310168326A CN 116165896 A CN116165896 A CN 116165896A
- Authority
- CN
- China
- Prior art keywords
- aircraft
- control
- identification
- instruction
- frequency domain
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 99
- 238000012937 correction Methods 0.000 claims abstract description 12
- 230000004044 response Effects 0.000 claims abstract description 10
- 230000003068 static effect Effects 0.000 claims abstract description 6
- 230000000694 effects Effects 0.000 claims description 28
- 230000003044 adaptive effect Effects 0.000 claims description 24
- 238000004088 simulation Methods 0.000 claims description 22
- 230000005284 excitation Effects 0.000 claims description 15
- 230000005484 gravity Effects 0.000 claims description 11
- 238000006243 chemical reaction Methods 0.000 claims description 10
- 230000006870 function Effects 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 8
- 230000001629 suppression Effects 0.000 claims description 8
- 238000012546 transfer Methods 0.000 claims description 8
- 230000009466 transformation Effects 0.000 claims description 8
- 230000008859 change Effects 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 238000005312 nonlinear dynamic Methods 0.000 claims description 2
- 230000010363 phase shift Effects 0.000 claims description 2
- 238000005316 response function Methods 0.000 claims description 2
- 238000000926 separation method Methods 0.000 claims description 2
- 238000012545 processing Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 4
- 238000003775 Density Functional Theory Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 238000013528 artificial neural network Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000005279 excitation period Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
- G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
- G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
- G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
The invention belongs to the technical field of airplane control, and relates to an airplane self-adaptive control method based on online frequency domain recursion identification. The pneumatic derivative can be accurately identified by using a recursive identification method, and the identification error of the static coefficient and the rudder efficiency coefficient is not more than 1%. The improved self-adaptive dynamic inverse control introduces an advanced correction link on the basis of the traditional self-adaptive dynamic inverse control, so that the self-adaptive dynamic inverse controller is faster in response, can realize high-precision control under the condition of inaccurate pneumatic model by combining an on-line identification method, is applied to aircraft pitching attitude control, and has smaller overshoot, shorter rise time and shorter peak time and good control quality compared with the traditional controller and the common self-adaptive dynamic inverse controller.
Description
Technical Field
The invention belongs to the technical field of airplane control, and relates to an airplane self-adaptive control method based on online frequency domain recursion identification.
Background
The control system adopted by most of the existing aircrafts in service is a traditional controller, such as a PID controller, etc., but because the real aircraft pneumatic model has the characteristics of nonlinearity, the control performance of the traditional method cannot fully meet the expectations and the robustness is poor, meanwhile, the secondary problems of long development period, etc. of the control system are brought, and in order to overcome the problem, a great deal of research has been focused on nonlinear control system designs, such as an adaptive control method, an intelligent control method, etc.
Among the published studies, most of the pneumatic recognition methods are studied separately from the adaptive control methods, and there are plentiful results of the studies in both aspects.
The pneumatic identification function comprises correction of a pneumatic model, control reconstruction under the fault condition, real-time adjustment of control parameters and the like, and in some application scenes, the pneumatic model needs to meet the requirement of real-time identification, and along with the increase of accumulated flight data, accurate identification is gradually realized. The method adopted here is dynamic parameter identification, and the method for identifying the aircraft dynamics comprises two main types of methods, namely a time domain identification method and a frequency domain identification method, specifically comprises a regression method represented by a least square method, a plurality of methods are proposed at present based on a least square method, and the common identification methods further comprise a maximum likelihood method, a Kalman filtering method and the like. In terms of the frequency domain method, some practical methods have also been proposed. Furthermore, intelligent methods are also being used in a step-wise fashion in pneumatic recognition, with neural networks typically being used for dynamic recognition.
In the aspect of nonlinear self-adaptive control methods, there are many different technical routes, for example, a series control system combining a neural network or a fuzzy controller with a PID controller is used for realizing real-time adjustment of control gain based on control errors or other observables, the self-adaptive dynamic inverse method (NDI) control is also a self-adaptive control method, and the dynamic inverse method is a control method with stronger adaptability and universality.
The online pneumatic identification and the self-adaptive control method are combined, so that online adjustment is realized, and the online pneumatic identification and self-adaptive control method is a learning control method route proposed in recent years. Some research projects have also performed considerable technical exploration work on this control optimization concept. However, in these studies, online recursive application of the frequency domain has been studied with little effort, and the method of combining frequency domain identification with the improved adaptive dynamic inverse control method has been more involved.
Disclosure of Invention
In order to solve the problems, the invention provides an airplane self-adaptive control method based on online frequency domain recursion identification.
The technical scheme of the invention is as follows:
an aircraft self-adaptive control method based on online frequency domain recursion identification comprises 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 can be 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 in the body coordinate system. The variable superscript indicates the derivative of the variable, as follows.
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 force, F Ay Is a lateral force phi T For engine mounting angle F Ty For engine thrust at O b y b Component of the axis.
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 Respectively, the airplane is seated around the airplane bodyThe product of inertia of each axis is marked.
The angular momentum H equation is:
in the formula ,Hx 、H y 、H z The components of angular momentum in the body coordinate system, respectively.
The movement of the aircraft about the centroid can be expressed as:
in the formula ,LA 、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 over the center of mass, so that the moment is not generated.
(2) Frequency domain online recursion identification method
The frequency domain online recursion identification comprises three parts of multi-sinusoidal control surface excitation, data trend removal and CZT conversion, a specific flow chart is shown in fig. 3, after an aircraft reaches a preset flight state, excitation signals are applied to the control surface, and after sensor data are obtained, aerodynamic force or moment coefficients are calculated; and 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 obtaining an identification result, obtaining a final identification result, and judging whether the identification is converged or not.
The flow-based on-line recursion identification method is suitable for force or moment identification.
The contents of the three parts are described in detail below.
(2.1) multiple sinusoidal control surface excitation:
when the aircraft is in stable flight, the information content of the flight data is small, and the low-information-content flight data cannot support accurate identification of the aerodynamic coefficient of the aircraft, so that an effective maneuvering mode is required to be adopted, and the information content of the flight data is improved to meet the identification requirement.
The method adopted by the invention is a control surface excitation signal with multiple sine superposition inputs, and the control surface excitation signal can be superposed with a control instruction to be used as a final control surface deflection instruction. The form of the multi-sinusoidal excitation input is:
here, the disturbance input applied to the jth control plane is set to be 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 real-time trending method must be adopted, and a fourth-order Baud Wo Sigao pass filter is adopted for trending.
The filtering is performed by adopting a recursive method, which is a key point, because the method can realize the trend removal of the sensor data in real time, and finally achieves the effect of real-time identification.
The expression for the filter can be written as a differential equation of the rational transfer function as follows:
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 ,na 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:
the most commonly adopted method for converting time domain data into frequency domain is DFT method, but there are many disadvantages of DFT method, especially for aircraft dynamics identification, only a certain frequency band with a small range is usually needed to be concerned, and the DFT can introduce interference information, so that the CZT method is adopted here, and the data frequency domain transformation effect of the method is shown in FIG. 4.
Starting point z of the CZT transformation 0 The method can be arbitrarily selected, so that narrow-band high-resolution analysis can be performed on input data from any frequency, and meanwhile, the frequency resolution can be adjusted according to the requirement, so that precise analysis can be performed on any frequency band of interest.
The step of the CZT transformation can be expressed as:
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 0 represents the initial sampling point radius, θ 0 Indicating the phase angle of the initial sampling point,
represents the equal division angle of two adjacent points, W 0 Representing the helix elongation.
Fourier transforms G (k) and H (k) are then calculated for G (n) and H (n), respectively, and Y (k) =g (k) H (k) is calculated.
Then calculate the L-point Fourier inverse transform X (z) of Y (k) k )：
The conventional 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 (ω) can be 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 commands for each variable, using time scale separation of the variables, nonlinear Dynamic Inversion (NDI) is employed in turn to generate commands for faster variables, see fig. 5 for the controller architecture. The dynamic inverse of the outermost loop is to convert the instruction of track inclination angle χ and track deflection angle γ into the instructions 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 χ The gain is controlled for the outer loop.
The track tilt derivative is expressed as:
the two types of the combination are simplified to obtain:
the solution formula can be obtained according to a first-order linear non-homogeneous differential equation:
in a certain time period, χ converges to the track inclination instruction χ c Chi (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.
Bringing the control command back into the kinetic equation, it is possible to obtain:
the solution formula can be obtained according to a first-order linear non-homogeneous differential equation:
the same time domain response function as the track angle form can be obtained, which can reach 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, can be obtained:
in the formula ,
the pitching moment is determined by the current identification model of the aircraft, namely, the current identification result obtained by the method of the step (2) and the pitching moment reversely solved by the sensor instruction; k (K) q For angular velocity control gain, M δ,c For controlling rudder deflection moment command.
Substituting control moment command into
Can be obtained in the kinetic equation of (a):
where M is the torque value obtained by the sensor.
(3.2) improving the adaptive interference suppression controller:
aiming at the problems of large error and strong hysteresis based on conventional dynamic inversion under the action of unsteady aerodynamic force, an improved self-adaptive dynamic inversion control method is introduced, and the control effect of the pneumatic identification result under the condition of error is improved. The controller structure is shown in fig. 6. 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,
the desired dynamic response for the corresponding variable for the lower right corner 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, the situation can be improved through a series lead correction link, and a slow loop expected system can be written as follows:
wherein G(s) is an advance correction device:
ω C omega is the complex plane pole position D For the zero position, the pole-zero right lower corner identifies the variable name representing the corresponding variable. When selecting parameters, ω C <ω D In the case of G(s) being the leading element, when ω C >ω D G(s) is a hysteresis, so ω should be taken here C <ω D 。
After the lead correction link is introduced, the lead correction link needs 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 solving for the nominal angular rate in the slow loop. The conventional 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:
the invention has the beneficial effects that:
the invention discloses a self-adaptive control method based on real-time frequency domain recursion identification, and focuses on the realizability of real-time application and the effect of the control method. The characteristic that the frequency domain aerodynamic parameter identification method is insensitive to noise makes the frequency domain aerodynamic parameter identification method have great value in practical application, and the introduced recursive trend method further reduces the calculation amount requirement. The pneumatic derivative can be accurately identified by using a recursive identification method, and the identification error of the static coefficient and the rudder efficiency coefficient is not more than 1%. The improved dynamic inverse control introduces an advanced correction link on the basis of the traditional dynamic inverse control, so that the dynamic inverse controller responds faster, and can realize high-precision control under the condition of inaccurate pneumatic model by combining an online identification method.
Drawings
FIG. 1 is a flow chart of an improved adaptive control method based on online frequency domain pneumatic recognition;
FIG. 2 is a flow chart of an improved adaptive control algorithm based on online frequency domain pneumatic recognition;
FIG. 3 is a flow chart of a method for identifying pneumatic parameters by online recursion of a frequency domain;
FIG. 4 is a schematic diagram of a chirped z-transform function;
FIG. 5 is a flow chart of an adaptive dynamic inverse control method;
FIG. 6 is a flow chart of a modified adaptive dynamic inverse control method;
FIG. 7 is a control surface excitation signal and component diagram;
FIG. 8 is a diagram of identifying time interval flight conditions;
FIG. 9 is a graph of identifying interval flight data trending results;
FIG. 10 is a graph showing the comparison of the static coefficient identification result with the true value;
FIG. 11 is a graph showing the comparison of the identification result of the rudder efficiency coefficient and the true value;
FIG. 12 is a comparison of PID and dynamic inverse step signal simulated tracking effects;
FIG. 13 is a comparison of PID and dynamic inverse step signal simulation tracking rudder deflection effects;
FIG. 14 is a comparison of PID and dynamic inverse step signal simulation tracking pitch angle rate effects;
FIG. 15 is a comparison of PID, dynamic inversion and modified dynamic inversion step signal simulation tracking effects;
FIG. 16 is a comparison of PID, dynamic inversion, and modified dynamic inversion step signal simulation tracking rudder deflection effects;
FIG. 17 is a comparison of PID, dynamic inverse and modified dynamic inverse step signal simulation tracking pitch rate effects.
Detailed Description
The following describes the embodiments of the present invention further with reference to the drawings and technical schemes.
The improved self-adaptive control method based on-line frequency domain pneumatic identification is shown in fig. 1 and 2.
In the simulation example of the method implementation of the part, the method generally comprises two major parts, namely a frequency domain online recursion identification part and an adaptive control part based on an online identification result. Aiming at frequency domain real-time identification, real-time flight data are acquired through real-time flight simulation, firstly, recurrence trend processing is carried out on the real-time data, analysis is carried out on a trending result, then, the frequency domain recurrence identification program is substituted, the accuracy of the identification result is observed, and the identification precision is determined. In the control simulation part, the control effects of a conventional PID controller and an adaptive dynamic inverse controller are compared, then an improved NDI controller with advanced correction is introduced, the control effects of the three are compared, and the robustness of the improved adaptive dynamic inverse controller is verified. Finally, a full-flow simulation including real-time frequency domain identification and improved self-adaptive dynamic inverse control is provided, and the method can meet the requirement of real-time operation.
(1) Inputting an initial state, giving a target state
TABLE 1 simulation initiation parameters
The target state is to keep the pitch angle 11 deg. flat flight.
(2) Establishing a flight dynamics model
1) Ground coordinate system
Taking a certain point of the ground (initial position of the airplane) as an origin O g ,O g x g The axis is positioned in the horizontal plane along the take-off direction and forwards is positive, O g y g The axis being in the horizontal plane and being in contact with O g x g The axis is vertical, right is positive, O g z g The axis is vertical to the ground and is positive downwards.
2) Machine body coordinate system
Taking the mass center of the airplane as the origin O of a machine body coordinate system b The coordinate system is fixedly connected with the plane, O b x b The axis coincides with the axis of the machine body and is positioned in the longitudinal symmetrical plane of the machine body, and the forward direction is positive, O b y b Perpendicular to the longitudinal symmetry plane, right is positive, O b z b Is positioned in the longitudinal symmetry plane and is O b x b Vertical, positive downward.
3) Velocity coordinate System/airflow coordinate System
Taking the centroid of the aircraft as the origin of coordinates O a ，O a x a Is consistent with the airspeed direction of the aircraft, and forwards is positive, O a z a The axis is vertical and is positioned in the longitudinal symmetrical plane, the downward direction is positive, O a y a Perpendicular to plane O a x a z a And positive to the right.
4) Track coordinate system
The track coordinate system is fixedly connected with the aircraft, and takes the mass center of the aircraft as an origin O s ，O s x s Pointing in the projection direction of the aircraft speed in the plane of longitudinal symmetry of the aircraft, O s z s The axis is in the plane of longitudinal symmetry of the aircraft and is perpendicular to the axis and points downwards the aircraft, O s y s The axis is perpendicular to the plane and directed to the right of the aircraft.
The invention uses several coordinate conversion relations, which are respectively:
1) Conversion matrix from ground coordinate system to machine body coordinate system
2) Conversion matrix from air flow coordinate system to machine body coordinate system
In the formula, theta, phi and phi respectively represent pitch angle, roll angle and yaw angle, alpha represents attack angle, and beta represents sideslip angle.
(1.2) six degree of freedom dynamics modeling of aircraft
The aircraft centroid dynamics equation can be 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 in the body coordinate system. The variable superscript indicates the derivative of the variable, as follows.
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 force, F Ay Is a lateral force phi T For engine mounting angle F Ty For engine thrust at O b y b Component of the axis.
The aircraft inertia matrix I is:
in the inertial matrix, I xx 、I yy 、I zz Three items respectively represent the aircraft around the engine body axis O b x b 、O b y b 、O b z b Moment of inertia of three axes, I xy 、I xz 、……I yz The product of inertia of the aircraft around each axis of the machine body coordinate system is respectively calculated.
The angular momentum H equation is:
in the formula ,Hx 、H y 、H z The three terms are components of angular momentum in the body coordinate system, respectively.
The movement of the aircraft about the centroid can be expressed as:
in the formula ,LA 、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 over the center of mass, so that the moment is not generated.
(3) Frequency domain on-line identification flow simulation example
In this embodiment, the pitch moment coefficient of the aircraft is taken as an example for simulation. The pitch moment coefficient of an aircraft can be reduced by linearization into the following form, where the moment is of primary concern, which can be considered as being related to some amount of flight state and satisfies a transient linear relationship:
in the formula ,Cl 、C m 、C n The right lower corner mark of the three variables represents the derivative of the moment coefficient to the variable, the right lower corner mark is 0 represents the moment coefficient caused by the body configuration, b represents the span length, and c represents the average aerodynamic chord length.
Sometimes, the modeling term on the right side of the equation also comprises some high-order terms or coupling terms of different variables, and proper adjustment is needed to be made for different aircraft configurations to accurately model, and trade-off is made between modeling variable number and dynamic corresponding information quantity, so that the aims of accurately identifying and reducing calculated quantity as much as possible are achieved.
The analysis is mainly performed by taking pitching moment as an example, and the variable X is defined as a vector consisting of data directly acquired from an aircraft sensor, and comprises the following contents:
the recurrence trend of the data can be realized by the trend removal method, and the vector after finishing the trend removal is defined as X d The content of it is expressed as follows:
the trended data can be subjected to recursive Fourier transform by adopting the method to obtain
The expression is as follows:
in the frequency domain, the problem solved by recursion can be expressed in the form of a recursive least squares, expressed as follows:
the value of the identification result can be found by the following expression:
finally, static information filtered in the trend link needs to be recovered:
taking longitudinal attitude control as an example, after identification is finished, control surface instruction solving can be performed based on an identification result, wherein the formula is as follows:
the control surface is applied with an active excitation signal to excite the dynamic characteristics of the aircraft, so that the flight data contains more information, and the identification accuracy is improved.
Where T is the excitation period, and the corresponding image is shown in fig. 7:
the selected identification stage flight state of the invention is shown in fig. 8:
the flying height is about 5000m, the speed is about 150m/s, and the dynamic identification is carried out in a horizontal flying state.
Using a 4 th order bute Wo Sigao pass filter, the filter-related information is expressed as follows:
table 2 filter states
Firstly, the data needs to be subjected to trending treatment, and in order to realize real-time identification, a recursion trending method is adopted, and the comparison of the data before and after treatment is shown in fig. 9.
In the figure, the solid line represents raw data that is not distinguishable when there is a constant component or a component that is linearly related to time. In the above figures, it is evident that there is a significant constant component of both the angle of attack and the elevator angle, and therefore these data must be processed. The other two terms are relatively small in deviation, but also need to be processed.
In the figure, the dotted line is a result of real-time recurrence trend, in order to more intuitively observe the trend removal effect, after the simulation is finished, a more accurate batch processing trend removal method is adopted offline to perform data processing, and all the results are displayed in one picture, and it must be pointed out that the batch processing method does not meet the requirement of real-time operation, because it needs to process all flight data simultaneously, if batch processing is performed on all previous data after each data acquisition, the effect similar to the recurrence method can be achieved, but the calculation amount is very large, and the requirement on the hardware calculation capability of an airplane is very high. As is apparent from the graph, the recursive method requires about 20 seconds to converge the mole with the same accuracy as the batch method, because at the initial time, the flight data is very small, and the constant component contained in the flight data cannot be accurately identified.
According to the embodiment, when the recursive real-time trending method is adopted, after trending processing is started, a period of time is required to wait for the recursive pneumatic identification to be performed, the specific waiting time length is determined through flight data analysis, and the waiting time adopted is 20 seconds.
After the trend removal is completed, frequency domain recurrence identification can be performed, and the identification results of fig. 10 and 11 are:
in the stage of waiting for convergence of the real-time filter, the identification result is artificially set to be 0, in order to ensure that the identification result has larger fluctuation at the beginning of recursion, the first 90 samples are taken for batch processing once, the result is taken as an initial value for subsequent recursion, and in the stage of waiting for data collection of batch processing, the identification result is artificially set to be 0.
In the figure, it can be seen that the initial value error obtained after batch processing is very small, and no deviation phenomenon occurs in the subsequent recursion process, so that the identification is rapid, the convergence is good, and the error is small. To more clearly illustrate the accuracy of the identification, a period of time after the start of the recursion is selected, the relative errors of the parameter identification are averaged, and the calculated results are shown in the following table:
TABLE 3 relative error in pneumatic parameter identification
(4) Self-adaptive dynamic inverse control simulation based on real-time pneumatic parameter identification
Firstly, performing simulation verification on a conventional configuration self-adaptive dynamic inverse controller, and exploring the control effect comparison of the self-adaptive dynamic inverse controller and a PID controller under a step instruction. The pneumatic identification link of the last part is firstly carried out, the conversion of the control method is carried out after the identification is completed by two parts of recursion trend and real-time recursion identification, and finally, the tracking observation tracking effect of the step pitch angle instruction is respectively carried out under the PID and self-adaptive dynamic inverse control modes by giving the step pitch angle instruction. The simulation of this section is performed according to such a flow, but since the previous section has been studied and described in detail with respect to the accuracy of recognition, only the image of the instruction trace section is shown in this chapter.
As is apparent from simulation figures 12, 13 and 14, compared with the PID controller, the adaptive dynamic inverse controller has a faster response speed and a smaller overshoot, and as can be seen from the images of the deflection angle and the pitch angle of the elevator, the adaptive dynamic inverse controller can realize the deflection of the control surface at a larger angle and a faster adjustment response at the same time 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 values, three indexes of rise time, peak time and overshoot are selected to quantitatively compare the control effects of the two values.
Table 4 comparison of control effect of PID controller and adaptive dynamic inverse controller
Quantitative analysis shows that the self-adaptive dynamic inverse has better control effect than a PID controller on three indexes. Controller for controlling a power supply
(5) Improved adaptive dynamic inverse control simulation based on real-time pneumatic parameter identification
And (3) performing simulation verification of the improved self-adaptive dynamic inverse control method under the same condition, and likewise, adopting all links including recursion trend removal, identification and step pitch angle instruction, and only displaying images of a step instruction tracking stage.
FIGS. 15, 16, 17 illustrate the step response tracking effect of improved adaptive dynamic inverse control methods, conventional adaptive dynamic inverse control, and PID control.
As can be seen from the simulation diagram, compared with the PID controller and the self-adaptive dynamic inverse controller, the improved self-adaptive dynamic inverse controller has further improvement in the aspects of response speed and overshoot, and as can be seen from the elevator deflection angle image and the pitch angle speed image, the improved self-adaptive dynamic inverse controller has better response effect of the elevator deflection angle. In order to compare the difference between the two values, three indexes of rise time, peak time and overshoot are selected to quantitatively compare the control effects of the three indexes.
Table 5PID controller, adaptive dynamic inverse controller and improved adaptive dynamic inverse controller control effect comparison
From the above analysis, it can be concluded that: the improved self-adaptive dynamic inverse control has better control effect compared with the conventional self-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 ,Hx 、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 ,LA 、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 ,na 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 0 represents the initial sampling point radius, θ 0 Indicating the phase angle of the initial sampling point,
represents the equal division angle of two adjacent points, W 0 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:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310168326.4A CN116165896B (en) | 2023-02-27 | 2023-02-27 | Airplane self-adaptive control method based on online frequency domain recursion identification |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310168326.4A CN116165896B (en) | 2023-02-27 | 2023-02-27 | Airplane self-adaptive control method based on online frequency domain recursion identification |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116165896A true CN116165896A (en) | 2023-05-26 |
| CN116165896B CN116165896B (en) | 2023-10-20 |
Family
ID=86418030
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310168326.4A Active CN116165896B (en) | 2023-02-27 | 2023-02-27 | Airplane self-adaptive control method based on online frequency domain recursion identification |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116165896B (en) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6341247B1 (en) * | 1999-12-17 | 2002-01-22 | Mcdonell Douglas Corporation | Adaptive method to control and optimize aircraft performance |
| CN103995463A (en) * | 2014-05-30 | 2014-08-20 | 北京敬科海工科技有限公司 | Method for performing position servo driving on electro-hydraulic proportional valve based on hybrid control |
| CN104635492A (en) * | 2014-12-19 | 2015-05-20 | 中国科学院长春光学精密机械与物理研究所 | Parametric adaptive feed-forward control method of guide head stabilizing platform |
| CN110187713A (en) * | 2019-04-12 | 2019-08-30 | 浙江大学 | A kind of longitudinally controlled method of hypersonic aircraft based on aerodynamic parameter on-line identification |
| CN113625555A (en) * | 2021-06-30 | 2021-11-09 | 佛山科学技术学院 | Adaptive inverse control AGV rotation speed control method based on recursive subspace identification |
| CN115048724A (en) * | 2022-06-27 | 2022-09-13 | 南京航空航天大学 | B-type spline-based method for online identification of aerodynamic coefficient of variant aerospace vehicle |
| CN115169002A (en) * | 2022-07-06 | 2022-10-11 | 哈尔滨工业大学 | Time-varying parameter identification method for direct-air composite flight control system |
-
2023
- 2023-02-27 CN CN202310168326.4A patent/CN116165896B/en active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6341247B1 (en) * | 1999-12-17 | 2002-01-22 | Mcdonell Douglas Corporation | Adaptive method to control and optimize aircraft performance |
| CN103995463A (en) * | 2014-05-30 | 2014-08-20 | 北京敬科海工科技有限公司 | Method for performing position servo driving on electro-hydraulic proportional valve based on hybrid control |
| CN104635492A (en) * | 2014-12-19 | 2015-05-20 | 中国科学院长春光学精密机械与物理研究所 | Parametric adaptive feed-forward control method of guide head stabilizing platform |
| CN110187713A (en) * | 2019-04-12 | 2019-08-30 | 浙江大学 | A kind of longitudinally controlled method of hypersonic aircraft based on aerodynamic parameter on-line identification |
| CN113625555A (en) * | 2021-06-30 | 2021-11-09 | 佛山科学技术学院 | Adaptive inverse control AGV rotation speed control method based on recursive subspace identification |
| CN115048724A (en) * | 2022-06-27 | 2022-09-13 | 南京航空航天大学 | B-type spline-based method for online identification of aerodynamic coefficient of variant aerospace vehicle |
| CN115169002A (en) * | 2022-07-06 | 2022-10-11 | 哈尔滨工业大学 | Time-varying parameter identification method for direct-air composite flight control system |
Non-Patent Citations (1)
| Title |
|---|
| 鲁兴举 等: "基于递推傅里叶变换的飞行器参数在线辨识方法", 航空学报, vol. 35, no. 2, pages 532 - 540 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116165896B (en) | 2023-10-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107807663B (en) | Unmanned aerial vehicle formation maintaining control method based on self-adaptive control | |
| CN107357166B (en) | Model-free self-adaptive robust control method of small unmanned helicopter | |
| CN108873929B (en) | Method and system for autonomous landing of fixed-wing aircraft | |
| CN107563044B (en) | Four-rotor unmanned aerial vehicle path tracking control method based on online safety learning | |
| Guan et al. | Fixed-time control for automatic carrier landing with disturbance | |
| Li et al. | Full control of a quadrotor using parameter-scheduled backstepping method: implementation and experimental tests | |
| CN108594837A (en) | Model-free quadrotor drone contrail tracker and method based on PD-SMC and RISE | |
| CN108595756B (en) | Method and device for estimating flight interference of large envelope | |
| CN111367182A (en) | Hypersonic aircraft anti-interference backstepping control method considering input limitation | |
| CN111522352B (en) | Design method of single-parameter active disturbance rejection attitude controller of multi-rotor aircraft | |
| CN112180965A (en) | High-precision overload control method | |
| CN108638068A (en) | A kind of flying robot's Control System Design method carrying redundancy mechanical arm | |
| CN111026160A (en) | Trajectory tracking control method for quad-rotor unmanned aerial vehicle | |
| Huynh et al. | L 1 adaptive control for quadcopter: Design and implementation | |
| CN111142550B (en) | Civil aircraft aided driving control method and system and flight quality evaluation method | |
| CN111399473B (en) | Self-adaptive fault-tolerant control method and system and aerial robot | |
| CN114721266B (en) | Self-adaptive reconstruction control method under condition of structural failure of control surface of airplane | |
| CN115220467A (en) | Flying wing aircraft attitude control method based on neural network incremental dynamic inverse | |
| Yang et al. | Back-stepping robust control for flexible air-breathing hypersonic vehicle via α-filter-based uncertainty and disturbance estimator | |
| Cao et al. | Discrete-time incremental backstepping controller for unmanned aircrafts subject to actuator constraints | |
| Cheng et al. | Hover-to-cruise transition control for high-speed level flight of ducted fan UAV | |
| CN116165896B (en) | Airplane self-adaptive control method based on online frequency domain recursion identification | |
| CN116360255A (en) | Self-adaptive adjusting control method for nonlinear parameterized hypersonic aircraft | |
| CN113093782B (en) | Unmanned aerial vehicle designated performance attitude control method and system | |
| CN111061283A (en) | Air-breathing hypersonic aircraft height control method based on characteristic model |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |