Time-varying Parameter Identification Method for Direct Gas Composite Flight Control System

Abstract

Time-varying parameter identification method for direct gas composite flight control system, which belongs to the field of aircraft flight parameter identification. The present invention solves the problem that the constant aerodynamic parameters obtained by the small perturbation linearization method cannot meet the design requirements of the aircraft. The main technical solution adopted by the method of the present invention is: establishing the body coordinate system of the aircraft and the aircraft velocity coordinate system, establishing the pitch channel attitude dynamics equation based on the coordinate system, and designing a time-varying parameter filter based on the pitch channel attitude dynamics equation; then online identifying the aerodynamic parameters based on the designed filter system. The present invention transforms the time-varying parameter estimation into non-time-varying parameter estimation by polynomial fitting according to the Weierstrass approximation theorem. The method of the present invention can be applied to aircraft flight parameter identification.

Description

Technical Field

[0001] The present invention belongs to the field of aircraft flight parameter identification, and particularly relates to a method for identifying time-varying parameters of a direct-gas composite flight control system. Background Art

[0002] The identification of aircraft model parameters is an important basis and foundation for the design and development of new models of aircraft. Establishing an accurate aircraft dynamics model is an important part of aircraft control system design, flight quality evaluation, aircraft adaptive and self-learning control, and ground simulation. With the improvement of requirements and the development of aircraft aerodynamic design, the maneuverability of aircraft is constantly improving, and large-angle-of-attack and high-overload maneuvering actions are often made. Especially for aircraft adopting a direct-gas composite control strategy, when direct force is involved, it will have a greater impact on the aircraft's motion state, and the resulting non-linear characteristics of the aircraft's aerodynamic characteristics become more obvious. The constant aerodynamic parameters obtained by the small perturbation linearization method can no longer meet the design requirements of the aircraft. How to effectively identify them within the entire range of its time-varying aerodynamic parameters has become an urgent problem that flight controller designers need to solve. Summary of the Invention

[0003] The purpose of the present invention is to propose a method for identifying time-varying parameters of a direct-gas composite flight control system to solve the problem that the constant aerodynamic parameters obtained by the small perturbation linearization method cannot meet the design requirements of the aircraft.

[0004] The technical solution adopted by the present invention to solve the above technical problems is as follows:

[0005] A method for identifying time-varying parameters of a direct-gas composite flight control system, the method specifically includes the following steps:

[0006] Step 1: Design of a time-varying parameter filter system for direct-gas composite flight control

[0007] Step 1.1: Establish the aircraft body coordinate system ox1y1z1 and the aircraft velocity coordinate system ox3y3z3;

[0008] Step 1.2: Based on the coordinate systems in Step 1.1, establish the pitch channel attitude dynamics equation; the pitch channel attitude dynamics equation is:

[0009]

[0010] where ω z is the pitch rate, is the first derivative of ω z , δ z is the rudder deflection angle, F z is the direct force, α is the angle of attack, and the aerodynamic parameters a1, a2, and a3 are respectively defined as: Direct force parameter l z is defined as: l z =-l / J z where J z is the moment of inertia of the aircraft about the coordinate axis oz3, is the partial derivative of the pitching moment with respect to the pitching angular rate, is the partial derivative of the pitching moment with respect to the angle of attack, is the partial derivative of the pitching moment with respect to the rudder deflection angle, and l is the distance from the lateral thrust engine to the aircraft's center of mass;

[0011] Step 1-3: Design a time-varying parameter filter based on the pitching channel attitude dynamics equation in Step 1-2;

[0012] For the aerodynamic parameter a2, approximate a2(α) using a first-order polynomial in α, i.e.:

[0013] a2(α)≈a 20 +k 21 α (3)

[0014] where a 20 and k 21 are constants;

[0015] Substitute Equation (3) into the pitching channel attitude dynamics equation (1) to obtain:

[0016]

[0017] The state variables x1, x2, x3, x4, x5, and x6 of the filter system are defined as: x1 = ω z , x2 = a1, x3 = a 20 , x4 = k 21 , x5 = l z , x6 = a3;

[0018] Then the state vector X of the filter system is:

[0019] X = [x1 x2 x3 x4 x5 x6] T (5)

[0020] Define the state equation of the filter system as:

[0021]

[0022] where is the first derivative of X, and f(X) is the state transition function of the filter system;

[0023]

[0024] The measurement equation h of the filter system is:

[0025] h = [1 0 0 0 0 0] (8)

[0026] Step 2: Based on the filter system designed in Step 1, perform online identification on a 20 and k 21 ;

[0027] Similarly, identify the pneumatic parameters a1 and a3.

[0028] Furthermore, in Step 1, the specific process of establishing the aircraft body coordinate system ox1y1z1 and the aircraft velocity coordinate system ox3y3z3 is as follows:

[0029] Taking the center of mass of the aircraft as the origin o, the central symmetry axis of the aircraft as the ox1 axis, the aircraft velocity direction as the ox3 axis, the oy1 axis and the oy3 axis are both in the vertical plane of the aircraft, and the relationship between the oy1 axis and the oy3 axis is determined by the angle of attack and the yaw angle. The oz1 axis satisfies the right-hand rule with the ox1 axis and the oy1 axis, and the oz3 axis satisfies the right-hand rule with the ox3 axis and the oy3 axis.

[0030] Furthermore, the specific process of Step 2 is as follows:

[0031] Discretize Equation (7) to obtain a discrete system:

[0032]

[0033] where f(X k ) is the state transition function of the filter system at the k-th step, T is the sampling interval, and X k is the state vector of the filter system at the k-th step;

[0034] Then, based on the discrete system, perform online identification on the parameters a 20 and k 21 ;

[0035] Furthermore, the discretization of Equation (7) uses the difference method.

[0036] Furthermore, the online identification of the parameters a 20 and k 21 based on the discrete system uses the EKF algorithm.

[0037] Even further, the process of the EKF algorithm is as follows:

[0038] Step 1: Predict the state and the error covariance matrix:

[0039]

[0040]

[0041] Among them, is the one-step prediction result of the state, is the estimated result of the filter system state transition function at the (k - 1)-th step, P k|k-1 is the one-step prediction result of the error covariance matrix, P k-1|k-1 is the estimated result of the error covariance matrix at the (k - 1)-th step, F k is the Jacobian matrix at the k-th step, and the superscript T represents the transpose of the matrix, Q k is the process noise at the k-th step;

[0042] Step 2: Calculate the Kalman gain:

[0043]

[0044]

[0045]

[0046] Among them, z k is the actual measurement value at the k-th step, is the predicted measurement value calculated using is the residual, is the state transition matrix at the k-th step, R k is the measurement noise covariance matrix at the k-th step, S k is the intermediate variable matrix, and the superscript -1 represents the inverse of the matrix, K k is the gain matrix at the k-th step;

[0047] Step 3: Update the state and the error covariance matrix:

[0048]

[0049] P k|k =(I - K k H k )P k|k-1 (16)

[0050] Among them, is the state estimation result at the k-th step, P k|k is the estimated result of the error covariance matrix at the k-th step, and I is the identity matrix.

[0051] The beneficial effects of the present invention are:

[0052] The present invention proposes a new method for filter modeling for systems with time-varying parameters. According to the Weierstrass approximation theorem, through polynomial fitting, the problem of time-varying parameter estimation is transformed into the problem of non-time-varying parameter estimation. And direct force action is introduced in the simulation. Through EKF joint filtering, the time-varying aerodynamic parameters and direct force parameters are identified online. It can be seen from the simulation results that the method of the present invention can effectively estimate the time-varying aerodynamic parameters and direct force parameters, and can meet the requirements of aircraft design. Description of the Drawings

[0053] Figure 1 is a schematic diagram of the relationship between the aircraft velocity coordinate system and the aircraft body coordinate system;

[0054] In the figure, M z is the pitching moment, and β is the sideslip angle;

[0055] Figure 2 is a diagram of the change of the aircraft angle of attack;

[0056] Figure 3 is a diagram of the ignition of the lateral thrust engine;

[0057] Figure 4 is the comparison diagram of the estimated value, true value and measured value of the pitch rate ω z ;

[0058] Figure 5 is the comparison diagram of the estimated value and true value of a2;

[0059] Figure 6 is the estimation result diagram of the parameter a 20 ;

[0060] Figure 7 is the estimation result diagram of the parameter k 21 ;

[0061] Figure 8 is the estimation result diagram of a1;

[0062] Figure 9 is the estimation result diagram of a3;

[0063] Figure 10 is the estimation result diagram of l z ;

[0064] Figure 11 is a diagram of the change of the rudder deflection angle;

[0065] Figure 12 is a diagram of the change of the pitch angle and the pitch angle command. Detailed Implementation Manner

[0066] Specific implementation method 1. The time-varying parameter identification method for the direct-gas composite flight control system described in this implementation method specifically includes the following steps:

[0067] Step 1. Design of the time-varying parameter filter system for direct-gas composite flight control

[0068] Step 1-1. Establish the body coordinate system ox1y1z1 of the aircraft and the velocity coordinate system ox3y3z3 of the aircraft;

[0069] Step 1-2. Based on the coordinate system in Step 1-1, establish the attitude dynamics equation for the pitch channel; considering the aerodynamic nonlinear characteristics, the attitude dynamics equation for the pitch channel is:

[0070]

[0071] where, ω z is the pitch rate. The comparison of the estimated value, true value, and measured value of the pitch rate is as Figure 4 shown, is the first derivative of ω z , δ z is the rudder deflection angle. The change of the rudder deflection angle is as Figure 11 shown, F z is the direct force, α is the angle of attack, and the aerodynamic parameters a1, a2, and a3 are defined as follows: The direct force parameter l z is defined as: l z =-l / J z , J z is the moment of inertia of the aircraft about the coordinate axis oz3, is the partial derivative of the pitch moment with respect to the pitch rate, is the partial derivative of the pitch moment with respect to the angle of attack, is the partial derivative of the pitch moment with respect to the rudder deflection angle, and l is the distance from the side thrust engine to the center of mass of the aircraft;

[0072] Step 1-3. Design a time-varying parameter filter based on the pitch channel attitude dynamics equation in Step 1-2;

[0073] Taking the aerodynamic parameter a2 of the angle of attack as an example, considering that a2 is a function of the angle of attack α, that is, a2 = a2(α). According to the first law of the Weierstrass approximation theorem (Stone–Weierstrass Theorem), for any continuous function, a polynomial with the highest term degree not greater than n can definitely be found so that the difference between it and the continuous function is zero. That is, a continuous function on a closed interval can be uniformly approximated by a polynomial series. Representing this function by a polynomial, it can be written as:

[0074] a2 = a2(α) = a 20+k1α + k2α 2 +k3α 3 +... (2)

[0075] where a 20 , k1, k2,... are constants;

[0076] For the aerodynamic parameter a2, when the angle of attack varies within the critical angle of attack range, a2(α) can be approximated by a first-order polynomial of α, i.e.:

[0077] a2(α) ≈ a 20 +k 21 α (3)

[0078] where a 20 and k 21 are constants;

[0079] Substituting Equation (3) into the pitch channel attitude dynamics equation (1) gives:

[0080]

[0081] The state variables x1, x2, x3, x4, x5, and x6 of the filter system are defined as: x1 = ω z , x2 = a1, x3 = a 20 , x4 = k 21 , x5 = l z , x6 = a3;

[0082] Then the state vector X of the filter system is:

[0083] X = [x1 x2 x3 x4 x5 x6] T (5)

[0084] Define the state equation of the filter system as:

[0085]

[0086] where is the first derivative of X, and f(X) is the state transition function of the filter system;

[0087]

[0088] The measurement equation h of the filter system is:

[0089] h = [1 0 0 0 0 0] (8)

[0090] Step 2: Perform online identification of a 20 and k 21 based on the filter system designed in Step 1;

[0091] Similarly, the aerodynamic parameters a1 and a3 are identified (the identification method is the same as that of the aerodynamic parameter a2).

[0092] For a direct gas composite control aircraft, the flight state of the aircraft will change violently due to the ignition of the lateral thrust engine. Its aerodynamic parameters are no longer constant as in traditional aircraft, but time-varying. The problem of identifying time-varying aerodynamic parameters has always been the main obstacle to designing an aircraft controller. The present invention aims at the problem of estimating time-varying parameters of a direct gas composite flight control system, and proposes a new engineering solution algorithm. By polynomial fitting, the time-varying parameters are fitted into several polynomials composed of constant parameters for identification, so as to transform the problem of identifying time-varying parameters into the problem of identifying constant parameters and reduce the implementation difficulty. The aircraft of the present invention can be a missile.

[0093] Specific Embodiment 2: Combine Figure 1 This embodiment is described. The difference between this embodiment and Specific Embodiment 1 is that in step 11, the specific process of establishing the aircraft body coordinate system ox1y1z1 and the aircraft velocity coordinate system ox3y3z3 is as follows:

[0094] Taking the centroid of the aircraft as the origin o, taking the central symmetry axis of the aircraft as the ox1 axis, taking the aircraft velocity direction as the ox3 axis, both the oy1 axis and the oy3 axis are in the vertical plane of the aircraft, and the relationship between the oy1 axis and the oy3 axis is determined by the angle of attack and the yaw angle. The oz1 axis satisfies the right-hand rule with the ox1 axis and the oy1 axis, and the oz3 axis satisfies the right-hand rule with the ox3 axis and the oy3 axis.

[0095] Other steps and parameters are the same as those in Specific Embodiment 1.

[0096] Specific Embodiment 3: The difference between this embodiment and Specific Embodiment 1 or 2 is that the specific process of step 2 is as follows:

[0097] Discretize Equation (7) to obtain a discrete system:

[0098]

[0099] where, f(X k ) is the filter system state transition function at the k-th step, T is the sampling interval, and X k is the filter system state vector at the k-th step;

[0100] Then, based on the discrete system, the parameters a 20 and k 21 are identified online.

[0101] Other steps and parameters are the same as those in Specific Embodiment 1 or 2.

[0102] Specific Embodiment 4: Different from one of Specific Embodiments 1 to 3, the discretization of formula (7) is carried out by the difference method.

[0103] Other steps and parameters are the same as those in one of Specific Embodiments 1 to 3.

[0104] Specific Embodiment 5: Different from one of Specific Embodiments 1 to 4, for the online identification of parameters a 20 and k 21 based on the discrete system, the EKF algorithm is adopted.

[0105] Other steps and parameters are the same as those in one of Specific Embodiments 1 to 4.

[0106] Specific Embodiment 6: Different from one of Specific Embodiments 1 to 5, the process of the EKF algorithm is as follows:

[0107] Step 1: Predict the state and the error covariance matrix:

[0108]

[0109]

[0110] Among them, is the one-step prediction result of the state, is the estimation result of the filter system state transition function at the (k - 1)-th step, P k|k-1 is the one-step prediction result of the error covariance matrix, P k-1|k-1 is the estimation result of the error covariance matrix at the (k - 1)-th step, F k is the Jacobian matrix at the k-th step, the superscript T represents the transpose of the matrix, Q k is the process noise at the k-th step;

[0111] Step 2: Calculate the Kalman gain:

[0112]

[0113]

[0114]

[0115] Among them, z k is the actual measurement value at the k-th step, is the predicted measurement value calculated using is the residual, is the state transition matrix at the k-th step, R k is the measurement noise covariance matrix at the k-th step, S k kis an intermediate variable matrix, and the superscript -1 represents the inverse of the matrix, K k is the gain matrix at the k-th step;

[0116] Step 3: Update the state and error covariance matrices:

[0117]

[0118] P k|k =(I - K k H k )P k|k-1 (16)

[0119] where, is the state estimation result at the k-th step, and P k|k is the error covariance matrix estimation result at the k-th step, and I is the identity matrix.

[0120] Other steps and parameters are the same as those in any one of the specific embodiments 1 to 5.

[0121] Before identifying the time-varying parameter a2(α) of the present invention, it is necessary to analyze the observability of the nonlinear filter. The specific process is as follows:

[0122] According to the weak observability theory of the nonlinear system, the system has 6 state variables and 1 observable variable, so it is necessary to calculate the fifth-order Lie derivative:

[0123] Based on the nonlinear system formed by the above formula

[0124] The zero-order Lie derivative is:

[0125] d0 = H = x1 (17)

[0126] The first-order Lie derivative is:

[0127]

[0128] The second-order Lie derivative is:

[0129]

[0130] The third-order Lie derivative is:

[0131]

[0132] The fourth-order Lie derivative is:

[0133]

[0134] The fifth-order Lie derivative is:

[0135]

[0136] Let:

[0137]

[0138] Then the observation matrix can be written as:

[0139]

[0140] After the system observability analysis, through the filter design method of the present invention, the identification problem of the time-varying parameter a2(α) is transformed into the identification problems of two non-time-varying parameters a 20 and k 21 so that the on-line identification of the time-varying aerodynamic parameters can be carried out.

[0141] Simulation part

[0142] The Kalman filter of the present invention approximates the time-varying parameter as an nth-order polynomial through the Weierstrass approximation theorem, so as to realize the transformation of the time-varying parameter estimation problem into a constant parameter estimation problem. The specific implementation scheme of the present invention is explained below through a simulation example:

[0143] In the simulation experiment, let the process noise matrix Q = 0.001·I 6*6 , the measurement noise R = 0.01, the true aerodynamic parameter a1 is 0.075, a3 is 24.5, and the true direct force parameter l z is 0.05. According to past experimental experience, the aerodynamic parameter a2 increases with the increase of the angle of attack α, and when the angle of attack changes from 0° to 30°, the value of a2 increases to twice the initial value. The change of the angle of attack is as Figure 2 shown, so the designed a 20 true value is 10, and k 21 true value is 18. In order to better reflect the filtering effect, let the pitch angle command θ c of the aircraft projectile change as shown in Equation (25). The change of the pitch angle command is as Figure 12 shown, that is:

[0144] θ c = p1cos(p2t) + p3t (25)

[0145] where p1, p2, and p3 are adjustable parameters. To estimate the direct force parameter l z and to conform to the actual flight conditions of the aircraft, an open-loop pulsed lateral force signal is applied to the aircraft, and the ignition condition of the side thrust engine is as Figure 3 shown, its signal amplitude is 200N, the ignition duration is 4ms, the interval between two ignitions is 36ms, and the pulse engine startup time is from 1s to 5s after the start of the simulation.

[0146] The filter state variable is defined as x1 = ω z , x2 = a1, x3 = a20 , x4 = k 21 , x5 = l z , x6 = a3, then the system state vector is:

[0147] X = [x1 x2 x3 x4 x5 x6] T (26)

[0148] Define the state equation as:

[0149]

[0150] Among them, the system state transition function is:

[0151]

[0152] The measurement equation is:

[0153] h = [1 0 0 0 0 0] (29)

[0154] According to the difference method, the discrete system X k+1 = f(X k );

[0155] Among them:

[0156]

[0157] Then the EKF algorithm process is as follows:

[0158] 1) Predict the state and error covariance matrix:

[0159]

[0160]

[0161] 2) Calculate the Kalman gain:

[0162]

[0163]

[0164]

[0165] Among them is the state transition matrix, R k is the measurement noise covariance matrix.

[0166] 3) Update the state and error covariance matrix:

[0167]

[0168] P k|k=(I - K k H k )P k|k-1 (37)

[0169] By the above filter design method, the identification problem of the time-varying parameter a2(α) is transformed into the identification problems of two non-time-varying parameters a 20 and k 21 . Thus, the on-line identification of the time-varying aerodynamic parameters can be carried out. The comparison between the estimated value and the true value of a2 is as shown in Figure 5 . The estimated result of the parameter a 20 is as shown in Figure 6 . The estimated result of the parameter k 21 is as shown in Figure 7 . Similarly, for a1, a3 and l z , the estimated result of a1 is as shown in Figure 8 . The estimated result of a3 is as shown in Figure 9 . The estimated result of l z is as shown in Figure 10 .

[0170] The above numerical examples of the present invention are only to illustrate in detail the calculation model and calculation process of the present invention, rather than to limit the implementation manners of the present invention. For those of ordinary skill in the art, other different forms of changes or modifications can be made on the basis of the above description. It is impossible to list all the implementation manners here. Any obvious changes or modifications derived from the technical solutions of the present invention still fall within the protection scope of the present invention.

Claims

1. A time-varying parameter identification method for a direct gas composite flight control system, characterized in that The method specifically includes the following steps: Step 1: Design of a time-varying parameter filter system for direct gas composite flight control Step 1.1: Establish the body coordinate system ox1y1z1 of the aircraft and the velocity coordinate system ox3y3z3 of the aircraft Step 1.2: Based on the coordinate systems in Step 1.1, establish the attitude dynamics equation of the pitch channel; the attitude dynamics equation of the pitch channel is: where ω z is the pitch rate, is the first derivative of ω z , δ z is the rudder deflection angle, F z is the direct force, α is the angle of attack, and the aerodynamic parameters a1, a2, and a3 are defined as follows: The direct force parameter l z is defined as: l z =-l / J z , J z is the moment of inertia of the aircraft about the coordinate axis oz3, is the partial derivative of the pitching moment with respect to the pitch rate, is the partial derivative of the pitching moment with respect to the angle of attack, is the partial derivative of the pitching moment with respect to the rudder deflection angle, and l is the distance from the side thruster to the center of mass of the aircraft; Step 1.3: Design a time-varying parameter filter based on the attitude dynamics equation of the pitch channel in Step 1.2 For the aerodynamic parameter a2, approximate a2(α) using a first-order polynomial with respect to α, that is: a2(α)≈a 20 +k 21 α (3) where a 20 and k 21 are constants; Substitute Equation (3) into the attitude dynamics equation (1) of the pitch channel to obtain: The state variables x1, x2, x3, x4, x5, and x6 of the filter system are respectively defined as: x1 = ω z , x2 = a1, x3 = a 20 , x4 = k 21 , x5 = l z , x6 = a3; Then the state vector X of the filter system is: X = [x1 x2 x3 x4 x5 x6] T (5) Define the state equation of the filter system as: Among them, is the first derivative of X, and f(X) is the state transition function of the filter system; The measurement equation h of the filter system is: h＝[1 0 0 0 0 0] (8) Step 2: Based on the filter system designed in Step 1, perform online identification on a 20 and k 21 ; Similarly, identify the aerodynamic parameters a1 and a3 2. The time-varying parameter identification method for the direct gas composite flight control system according to claim 1, characterized in that In Step 1.1, the specific process of establishing the body coordinate system ox1y1z1 of the aircraft and the velocity coordinate system ox3y3z3 of the aircraft is: Take the centroid of the aircraft as the origin o, take the central symmetry axis of the aircraft as the ox1 axis, take the aircraft velocity direction as the ox3 axis, both the oy1 axis and the oy3 axis are in the vertical plane of the aircraft, and the relationship between the oy1 axis and the oy3 axis is determined by the angle of attack and the yaw angle. The oz1 axis satisfies the right-hand rule with the ox1 axis and the oy1 axis, and the oz3 axis satisfies the right-hand rule with the ox3 axis and the oy3 axis 3. The time-varying parameter identification method for the direct gas composite flight control system according to claim 2, characterized in that The specific process of Step 2 is: Discretize Equation (7) to obtain a discrete system: Among them, f(X k ) is the filter system state transition function at the k-th step, T is the sampling interval, and X k is the filter system state vector at the k-th step; Then, based on the discrete system, the parameters a 20 and k 21 are identified online.

4. The time-varying parameter identification method for the direct gas composite flight control system according to claim 3, characterized in that, The discretization of Equation (7) uses the difference method 5. The time-varying parameter identification method for the direct gas composite flight control system according to claim 4, wherein The online identification of parameters a 20 and k 21 is carried out based on a discrete system, and the EKF algorithm is adopted.

6. The time-varying parameter identification method for the direct gas composite flight control system according to claim 5, characterized in that, The process of the EKF algorithm is: Step 1: Predict the state and error covariance matrix Among them, is the one-step prediction result of the state, is the estimation result of the filter system state transition function at the (k - 1)-th step, P k|k-1 is the one-step prediction result of the error covariance matrix, P k-1|k-1 is the estimation result of the error covariance matrix at the (k - 1)-th step, F k is the Jacobian matrix at the k-th step, and the superscript T represents the transpose of the matrix, Q k is the process noise at the k-th step; Step 2: Calculate the Kalman gain where z k is the actual measurement value at the k-th step, is the predicted measurement value calculated using , is the residual, is the state transition matrix at the k-th step, R k is the measurement noise covariance matrix at the k-th step, S k is the intermediate variable matrix, and the superscript -1 represents the inverse of the matrix, K k is the gain matrix at the k-th step; Step 3: Update the state and error covariance matrix P k|k = (I - K k H k )P k|k-1 (16) Among them, is the state estimation result at the k-th step, and P k|k is the error covariance matrix estimation result at the k-th step, and I is the identity matrix.

Images