A method for identifying time-varying global aerodynamic parameters of an aircraft

Abstract

The application discloses a kind of aircraft global time-varying aerodynamic parameter identification method, comprising: S1. establishing aircraft rigid body dynamics model, and obtaining time-varying parameter to be identified;S2. the time-varying parameter to be identified is written as polynomial form about time;S3. state sensitivity, covariance sensitivity, gain sensitivity, filter state sensitivity and covariance sensitivity are estimated;S4. according to aerodynamic parameter characteristics or time period, carry out local subarea;S5. in current interval, obtain state estimated value, estimated covariance matrix, Kalman gain, state filter value, covariance matrix by parameter value in combination with Kalman filter;S6. according to the above obtained criterion function, when criterion function is not minimum value, parameter is updated using Newton-Raphson algorithm;S7. using moving average filter to smooth all parameters, obtain a smooth global aerodynamic data curve.The application avoids the problem that fitting is not accurate when parameter variation amplitude is too large.

Description

Technical Field

[0001] This invention relates to the field of system identification technology, and more specifically to a method for identifying global time-varying aerodynamic parameters of an aircraft. Background Technology

[0002] During actual flight, due to its high speed and the high dynamic pressure environment, the external disturbances are very significant, and the interference between different parts of the aircraft is obvious. Therefore, the mathematical model of the aircraft is complex and variable, and the values ​​of aerodynamic parameters are not fixed constants, but variable parameters that change with factors such as angle of attack, sideslip angle, Mach number, and Reynolds number.

[0003] Currently, there are two methods for identifying globally time-varying aerodynamic parameters: one uses different aerodynamic models for different aerodynamic characteristics, and the other is local linearization, which assumes that the changes in aerodynamic characteristics are linear over a short period of time. Traditional methods lack the ability to identify variable parameters. Previous work combined the existing output error method with a "moving window" approach. This involves refreshing the data in the current sliding window at each sampling time, assuming that the aerodynamic parameters are constant within each window, and estimating the aerodynamic parameters using the output error method based on the maximum likelihood criterion. This process is repeated to estimate time-varying parameters. While this method can identify time-varying parameters to some extent, it relies heavily on the length of the "moving window" and has poor anti-interference capabilities, especially when the observed data is close to zero. Summary of the Invention

[0004] The purpose of this invention is to propose a method for identifying global time-varying aerodynamic parameters of aircraft that combines the filtering error method based on the maximum likelihood criterion with the approximation theorem. This method avoids the problem of inaccurate fitting when the parameter variation range is too large.

[0005] To achieve the above objectives, this application provides a method for identifying global time-varying aerodynamic parameters of an aircraft, comprising:

[0006] S1. Establish a rigid body dynamics model of the aircraft and obtain the time-varying parameters to be identified;

[0007] S2. Write the time-varying parameter to be identified in polynomial form with respect to time, and substitute the polynomial into the rigid body dynamics model of the aircraft to obtain the state equation;

[0008] S3. Estimate the state sensitivity, covariance sensitivity, gain sensitivity, filter state sensitivity, and covariance sensitivity using the state equation;

[0009] S4. Divide the area into local zones based on aerodynamic parameter characteristics or time period.

[0010] S5. In the current interval S i(i = 1 2 ... n) through parameter value ξ m (m refers to the number of iterations in the Newton-Raphson algorithm) combined with Kalman filtering to obtain the state prediction value. Predicted covariance matrix P(t|t) j-1 Kalman gain K(j), state filter value Covariance matrix P(j|j), (j=1 2 ... k), where j represents the current time and k represents the length of the time series of the current time period;

[0011] S6. Based on the state estimate Predicted covariance matrix P(t|t) j-1 Kalman gain K(j), state filter value The covariance matrix P(j|j) yields the criterion function J. m When the criterion function J m When the value is not a local minimum, use the Newton-Raphson algorithm to update the parameters (ξ). s ) m+1 Adjust (ξ) s ) m For (ξ) s ) m+1 Return to S6; when the convergence condition |J is satisfied. m+1 -J m When |<ε1, J m+1 The corresponding parameter (ξ) s ) m+1 This is the optimal estimate for the current interval; ε1 is the convergence threshold.

[0012] S7. Apply a moving average filter to all parameters (ξ) s ) i (i = 1 2 ... n, where n refers to the total length of the entire time series) is smoothed to obtain a smooth global aerodynamic data curve.

[0013] Furthermore, the rigid body dynamics model of the aircraft is specifically as follows:

[0014]

[0015] y(t)=g[x(t),u(t),ξ(t),t]+v(t)

[0016] Where x represents the state variables of the dynamic system, u represents the control variables, ξ represents the parameters to be identified, w represents the process noise of the dynamic system, y represents the output equation of the dynamic system, and v represents the observation noise vector. Let x be the first derivative.

[0017] Furthermore, the state equation is:

[0018]

[0019] y(t) s =g[x(t),u(t),ξ s ,t]+v(t)

[0020] Where the fixed parameter ξ s =(a2,a1,a0,b2,b1,b0,).

[0021] Furthermore, the criterion function expression is as follows:

[0022]

[0023] Wherein, the new information v is represented as:

[0024]

[0025] The innovation matrix B is represented as follows:

[0026] B(j)=h(j)P(t j |t j-1 )h T (j)+R v (j)

[0027] Where y(j) represents the observation, u(j) represents the input, h(j) represents the linearized observation matrix, and Rij v The variance matrix represents the observation noise.

[0028] Furthermore, the Newton-Raphson algorithm is used to update the parameters (ξ). s ) m+1 The implementation method is as follows:

[0029] (ξ s ) m+1 =(ξ s ) m +Δ(ξ s ) m

[0030]

[0031] Furthermore, the criterion function with respect to parameter ξ s The first and second derivatives are:

[0032]

[0033]

[0034] Where l = 1, 2, ..., p; p represents the parameter vector ξ s The length of the first and second derivatives; the sensitivity to innovation contained in the first and second derivatives. With the sensitivity of the innovation matrix These two sensitivities are obtained through the predicted state sensitivity. Predicting covariance sensitivity Gain Sensitivity Filter state sensitivity Covariance Sensitivity Obtain.

[0035] Furthermore, the predicted state sensitivity is:

[0036]

[0037] The sensitivity for predicting covariance is:

[0038]

[0039] in, The acquisition method is as follows:

[0040]

[0041]

[0042]

[0043] in, The acquisition method is as follows:

[0044]

[0045] The sensitivity of the innovation matrix is:

[0046]

[0047]

[0048]

[0049] The gain sensitivity is:

[0050]

[0051] The sensitivity of the new information is:

[0052]

[0053] The filter state sensitivity is:

[0054]

[0055] The above formula yields the filter sensitivity at the current moment. This value will be used as the initial value for the integral of the filter state sensitivity equation at the next moment;

[0056] The covariance sensitivity is:

[0057]

[0058] The above formula yields the filter covariance sensitivity at the current time. This value will be used as the initial value for predicting the integral of the covariance sensitivity equation at the next time step.

[0059] Furthermore, the moving average filter moves along a window of length windowSize and obtains the average value of the data contained in each window; the output of the moving average filter is:

[0060]

[0061] Wherein, input is the input vector.

[0062] Compared with existing technologies, the technical solutions adopted in this invention have the following advantages: This application proposes a global time-varying aerodynamic parameter identification method for launch vehicles that combines the output error method with the approximation theorem, based on prior knowledge of aerodynamic parameter changes. By referencing prior knowledge of aerodynamic parameter variations and dividing flight data into several intervals according to the strength of parameter time-varying characteristics, identification is performed within each interval. This method offers higher fitting accuracy and avoids the problem of inaccurate fitting when parameter variations are too large. Furthermore, by performing nonlinear fitting on different regions based on prior knowledge, the time-varying parameter identification is transformed into polynomial coefficients with respect to time, avoiding the problem of insufficient effective information caused by using the "moving window" output error method where the window cannot be too large. This method has a very wide range of applications and value. Attached Figure Description

[0063] Figure 1 This is a flight status data diagram for the first stage of a launch vehicle.

[0064] Figure 2 A diagram showing the aerodynamic coefficients for the first stage of a launch vehicle.

[0065] Figure 3 The data includes observation data and an observation noise graph.

[0066] Figure 4 Flowchart of the filtering error algorithm;

[0067] Figure 5 Here is a flowchart of the Newton-Raphson algorithm;

[0068] Figure 6 This is a schematic diagram of local time-varying parameter polynomial fitting. Specific implementation methods

[0069] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of this application and are not intended to limit the application; that is, the described embodiments are only a part of the embodiments of this application, and not all of them.

[0070] Example 1

[0071] Taking the aerodynamic coefficient identification of a certain launch vehicle in the first stage of flight as an example, a motion model of the rocket's center of mass was built based on project data. The rocket's aerodynamic parameters are fixed values ​​for the first 25 seconds, and then become time-varying parameters. The initial flight state is as follows: The specific steps for simulating a rocket's flight from 10 seconds to 114 seconds are as follows:

[0072] S1. Establish a rigid body dynamics model of the rocket and obtain the time-varying parameters to be identified;

[0073] Its rigid body dynamics model of the rocket is as follows:

[0074]

[0075] Mid-state variable x = [V θ σ xyz] T These represent the three positional parameters: velocity, trajectory inclination, and flight path angle, respectively. Input ξ(t) represents the parameter to be identified. In the dynamic equation, g represents gravitational acceleration, m represents mass, and P represents mass. e Indicates effective thrust. ψ represents the pitch and yaw program angles, q represents the dynamic pressure, and S represents the pitch and yaw program angles. m R represents the maximum cross-sectional area of ​​the rocket, R represents the Earth's radius, and r represents the distance from the rocket to the Earth's center.

[0076]

[0077] in, Let represent the transition matrix from the arrow system to the velocity system, and let α be the angle of attack and β be the sideslip angle.

[0078] S2. Write the time-varying parameter to be identified in polynomial form with respect to time. Here, we take a quadratic polynomial as an example:

[0079] C x (t)=a2t0 2 +a1t0+a0

[0080]

[0081] Substituting the above polynomials into the rocket rigid body dynamics model, we obtain the state equations;

[0082]

[0083] S3. Estimate the state sensitivity, covariance sensitivity, gain sensitivity, filter state sensitivity, and covariance sensitivity using the aforementioned state equation; the parameter to be identified is changed from the original C. x (t), C y (t) becomes a2, a1, a0, b2, b1, b0;

[0084] Specifically, a segmented identification method is adopted, assuming that the aerodynamic parameters are fixed in each time period. In this way, the fixed parameter identification algorithm in the previous section is combined with the Weierstrass approximation theorem to perform global time-varying parameter identification.

[0085] Taking a2 as an example, the state sensitivity equation is derived as follows:

[0086]

[0087]

[0088]

[0089]

[0090]

[0091]

[0092] S4. Divide the area into local zones based on aerodynamic parameter characteristics or time period. Set the window length L = 501 (first segment 0:0.01:5).

[0093] S5. In the current interval S i (i = 1 2 ... n) through parameter value ξ m Combine Kalman filtering to obtain state prediction values Predicted covariance matrix P(t|t) j-1 Kalman gain K(j), state filter value Covariance matrix P(j|j), (j=1 2 ... k), where j represents the current time and k represents the length of the time series of the current time period;

[0094] S6. Based on the state estimate Predicted covariance matrix P(t|t) j-1 Kalman gain K(j), state filter value The covariance matrix P(j|j) yields the criterion function J. m When the criterion function J mWhen the value is not a local minimum, use the Newton-Raphson algorithm to update the parameters (ξ). s ) m+1 Adjust (ξ) s ) m For (ξ) s ) m+1 Return to S5; when the convergence condition |J is met. m+1 -J m When |<ε1, J m+1 The corresponding parameter (ξ) s ) m+1 This is the optimal estimate for the current interval; ε1 is the convergence threshold.

[0095] Specifically, the correction value for each step is:

[0096]

[0097]

[0098]

[0099]

[0100]

[0101] The specific iterative process is as follows: Based on the wind tunnel test and theoretical calculation results of the aircraft (which can be based on previous calculations based on least squares aerodynamic parameters), initial values ​​of the aerodynamic parameters (ξ) are given. s 0, the rocket's state values, observation values, and sensitivity matrix are obtained by integrating the state equation, observation equation, and sensitivity equation. Then solve the system of linear algebraic equations to obtain Δ(ξ). s ), and then with (ξ s )1=(ξ s )0+Δ(ξ s )0 replaces the original (ξ) s )0; Repeat the above process until convergence, at which point the parameter (ξ) s ) k , which are the aerodynamic parameters we are looking for.

[0102] S7. Apply a moving average filter to all parameters (ξ) s ) i (i = 12...n, where n refers to the total length of the entire time series) is smoothed to obtain a smooth global aerodynamic data curve.

[0103] The foregoing description of specific exemplary embodiments of the invention is for illustrative and explanatory purposes. These descriptions are not intended to limit the invention to the precise forms disclosed, and it will be apparent that many changes and variations can be made in accordance with the foregoing teachings. The exemplary embodiments were chosen and described in order to explain the specific principles of the invention and its practical application, thereby enabling those skilled in the art to implement and utilize various different exemplary embodiments of the invention, as well as various different choices and variations. The scope of the invention is intended to be defined by the claims and their equivalents.

Claims

1. A method for identifying global time-varying aerodynamic parameters of an aircraft, characterized in that, include: S1. Establish a rigid body dynamics model of the aircraft and obtain the time-varying parameters to be identified; S2. Write the time-varying parameter to be identified in polynomial form with respect to time, and substitute the polynomial into the rigid body dynamics model of the aircraft to obtain the state equation; S3. Estimate the state sensitivity, covariance sensitivity, gain sensitivity, filter state sensitivity, and covariance sensitivity using the state equation; S4. Divide the area into local zones based on aerodynamic parameter characteristics or time period. ; S5. In the current interval Through parameter values Combine Kalman filtering to obtain state prediction values Predicting covariance matrix Kalman gain State filter value Covariance matrix , ,in Indicates the current moment. Indicates the length of the time series for the current time period; S6. Based on the state estimate Predicting covariance matrix Kalman gain State filter value Covariance matrix Obtain the criterion function When the criterion function When the value is not a minimum, use the Newton-Raphson algorithm to update the parameters. ,Adjustment for Return to S6; when the convergence condition is met. hour, Corresponding parameters This is the optimal estimate for the current interval; This is the convergence threshold; S7. Apply a moving average filter to all parameters. Smoothing is performed to obtain a smooth global aerodynamic data curve.

2. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 1, characterized in that, The rigid body dynamics model of the aircraft is as follows: in, For the state variables of the dynamic system, To control variables, For the parameters to be identified, For process noise in dynamic systems, The output equation of the dynamic system is... To observe the noise vector, for The first derivative.

3. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 1, characterized in that, The state equation is: Among them, fixed parameters .

4. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 1, characterized in that, The criterion function expression is as follows: Among them, new information Represented as: New Information Matrix Represented as: in, Representing the observed quantity, Indicates the input quantity. This represents the linearized observation matrix. The variance matrix represents the observation noise.

5. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 1, characterized in that, Update parameters using Newton-Raphson algorithm The implementation method is as follows: 。 6. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 1, characterized in that, The criterion function has parameters The first and second derivatives are: in, ; Represents the parameter vector The length of the first and second derivatives; the sensitivity to innovation contained in the first and second derivatives. With the sensitivity of the innovation matrix These two sensitivities are obtained by predicting state sensitivity. Predicting covariance sensitivity Gain and sensitivity Filter state sensitivity Covariance sensitivity Obtain.

7. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 6, characterized in that, The estimated state sensitivity is: The sensitivity for predicting covariance is: in, , The acquisition method is as follows: in, The acquisition method is as follows: The sensitivity of the innovation matrix is: The gain sensitivity is: The sensitivity of the new information is: The filter state sensitivity is: The above formula yields the filter sensitivity at the current moment. This value will be used as the initial value for the integral of the filter state sensitivity equation at the next moment; The covariance sensitivity is: The above formula yields the filter covariance sensitivity at the current time. This value will be used as the initial value for the integral of the covariance sensitivity equation at the next time step.

8. The method for identifying global time-varying aerodynamic parameters of an aircraft according to claim 1, characterized in that, The moving average filter is along the data shift length of... The window is defined, and the average value of the data contained in each window is obtained; the output of its moving average filter is: in, The input vector.

Images