Abstract
When the observed input-output data are corrupted by the observed noises in the aircraft flutter stochastic model, we need to obtain the more exact aircraft flutter model parameters to predict the flutter boundary accuracy and assure flight safety. So, here we combine the instrumental variable method in system identification theory and variance matching in modern spectrum theory to propose a new identification strategy: instrumental variable variance method. In the aircraft flutter stochastic model, after introducing instrumental variable to develop a covariance function, a new criterion function, composed by a difference between the theory value and actual estimation value of the covariance function, is established. Now, the new criterion function based on the covariance function can be used to identify the unknown parameter vector in the transfer function form. Finally, we apply this new instrumental variable variance method to identify the transfer function in one electrical current loop of flight simulator and aircraft flutter model parameters. Several simulation experiments have been performed to demonstrate the effectiveness of the algorithm proposed in this paper.
1. Introduction
Flutter is a large-scale vibration phenomenon in which an elastic structure is coupled by aerodynamics, elastic force, and inertial force in a uniform airflow. It is a most interesting issue in aeroelastic dynamics, which can damage aircraft structures and collapse buildings and bridges. The flutter phenomenon can occur during the flight of an aircraft. At this time, due to the action of the airflow, the elastic structure of the aircraft (such as the wing, tail, or operating surface) will generate additional aerodynamic forces. As an excitation force, additional aerodynamic forces will exacerbate structural vibration; meanwhile, the damping force of the air on the aircraft structure tries to weaken the vibration. At low speeds, the flutter after the disturbance gradually disappears due to the dominant damping force. When at a certain flight speed, that is, the critical speed of the flutter (the flutter boundary) is reached, the exciting force is dominant and the equilibrium position is unstable [1], which causes a large vibration or makes aircraft to be destroyed in a short time.
In order to avoid the occurrence of flutter accidents, the new aircraft development must undergo one flutter test to determine the stable flight envelope without flight flutter. The main content of the flutter test is to apply excitation to the aircraft structure under different flight conditions (different flight altitudes and speeds) and to identify the model parameters such as the flutter frequency and damping of the aeroelastic structure based on the dynamic response data. Then, flutter boundary will be predicted by virtue of the identified model parameters. Due to the complex dynamic characteristics of the aircraft structure, the dense model, and various disturbances (such as turbulence) and measurement noise during the flight process, these above factors likely cause many results to the dynamic response signal, such as low signal-to-noise ratio, short effective samples, and nonstationary processes. Then, it is very difficult to accurately identify the model parameters of aircraft flutter. Therefore, the problem on how to effectively process the test data and accurately identify the model parameters becomes an important research direction of the current flight flutter test; it means that the aircraft flutter model parameter identification can accurately predict the flutter boundary to ensure flight safety.
The identification of current flutter model parameters has attracted wide attention from various countries. For example, when developing a new type of drum-type small rotor excitation device, the United Kingdom used wavelet to filter the test signal in the time-frequency domain and applied advanced subspace methods to identify model parameters [2]. These new methods significantly improve the estimation accuracy of flutter model parameters and the accuracy of flutter boundary prediction [3]. It is particularly worth mentioning that the Free University of Brussels proposed a large number of new model parameter identification strategies based on the theory of frequency-domain system identification [4] and analyzed the accuracy of model parameter identification by means of the asymptotic theory of parameter estimation. It is well known that when to analyze and process input-output observation data in the frequency domain, firstly, discrete Fourier transform (DFT) must be used to transform the number of finite data samples for the input-output observation data. Usually, when doing discrete Fourier transform, for the simplicity of analysis, the initial state and terminal state of the observed data are ignored, that is, the negative influence of the transient term on the frequency-domain response function is neglected. This kind of neglect will inevitably affect the identification accuracy of the flutter model parameters [5]. So, it is necessary to introduce a variety of frequency-domain windowing functions to compensate or avoid the phenomenon of aliasing spectrum and leakage spectrum, generated after discrete Fourier transform [6]. The identification of aircraft flutter model parameters is less studied in China. In China, the whole flutter test process was systematically studied in [7], where a wavelet method for flutter test data processing was proposed to improve the signal-to-noise ratio of flight test under the condition of small rocket excitation. Then, aiming at improving the effect of pulse excitation response, a method based on support vector machine for flight test response data was proposed [8]. A wavelet time-frequency domain algorithm and a fractional Fourier domain method were proposed, respectively, for the rudder surface sweep excitation of the telex aircraft [9]. In order to make up for the shortcomings of the traditional least squares frequency-domain fitting identification algorithm, the global least squares identification algorithm in the frequency domain is proposed in [10], which avoids the complex nonlinear optimization and the dependence on initial value in the iterative algorithm. Based on the aforementioned four research contents, according to the whole framework of the system identification theory [11], the four research contents are just two aspects of the parameter estimation and experimental design in system identification theory. So far, there is a lack of works in the literature that seeks to solve the problem of parameter estimation of the flutter phenomenon.
In recent years, the authors conduct in-depth research on the identification of flutter model parameter based on the parameter estimation strategy. The research results are summarized as follows. In [12], the nonlinear separable least squares algorithm is used to accurately identify the flutter model parameters of the aircraft within the noise environment. Combined with the transfer function model, the identification problem with noisy systems is transformed into one nonlinear separable least squares problem, and the variance of the two noises and the unknown parameters in the transfer function are separately identifiable [13]. Furthermore, a simplified form of the maximum likelihood cost function of the aircraft flutter model is derived by means of the principle of frequency-domain maximum likelihood estimation [14]. In order to reduce the possibility of convergence to the local minimum, the iterative convolution smoothing identification method is derived by using global optimization theory. Based on the state-space model of the aircraft flutter stochastic model, the maximum value of the likelihood function is calculated in an iterative form using an expectation-maximization method suitable for the noise environment [15]. There is no need to calculate the second-order partial derivative and its approximation corresponding to the log-likelihood function, and then this can increase the likelihood that the likelihood function converges to a stationary point. Recently, the first author combines the bias-compensated algorithm and the instrumental variable algorithm to promote a new bias-compensated instrumental variable algorithm. Combined with the transfer function model, the identification problem of the noisy system is transformed into an iterative problem, which is used to solve a complex identification problem with white input noise and colored output noise. It is the generalization of the input-output observation noise, where the input noise and output noise are all white noises. Furthermore, the instrumental variable and subspace identification algorithm are combined to obtain the optimal instrumental variable subspace identification algorithm, which is used to accurately identify each system matrix in the aircraft flutter stochastic model under the random state-space form, and then the desired aircraft flutter model parameters are obtained. The detailed derivation of various methods for identifying the aircraft flutter model parameters proposed by the author can be referred to reference [11].
Based on the aforementioned research results, this paper continues to study the problem of identifying aircraft flutter model parameters in [12]. The main contribution of this paper is to combine instrumental variable method in system identification theory with variance matching method in modern spectrum estimation theory to form a new identification strategy—instrumental variable variance method. The instrumental variable method is derived from the system identification theory, accompanied by the development of the system identification discipline. At present, the research on instrumental variable identification is mature, and it has been successfully applied to closed-loop system identification and cascade system identification. The purpose of introducing instrumental variable is to ensure the independence between the two vectors, and then, the consistency and unbiased parameter estimates are guaranteed. Variance is a second-order statistical property used to describe random vectors in probability theory. It has been proved that second-order statistical property cannot fully reflect the statistical properties of random vectors, but higher-order statistics such as third-order and fractional orders are needed. Statistical characteristics are not limited to single time or frequency domain, but a combination of time domain and frequency domain, i.e.,time-frequency domain. In this paper, according to the stochastic model of aircraft flutter, by introducing instrumental variable to form the variance function, a certain norm of the difference between the theoretical value of the variance function and its actual estimated value is taken as a criterion function. Then, the criterion function based on the variance function is used to identify the unknown parameter vector. A stochastic model of aircraft flutter is represented as a model with input-output observation noise, which is used to describe the random factors generated in the flutter test. The unknown parameter vector in the transfer function form is obtained by minimizing this criterion function. As accurate transfer function estimation is the premise of model parameter identification, the detailed process of minimizing the criterion function is deduced, and the corresponding partial derivative expression is also given.
2. Problem Description
An important algorithm for flutter model parameter identification is a frequency-domain algorithm based on transfer function model. It usually needs to establish a parameterized frequency-domain transfer function (frequency response function) model and obtain the model parameters of each order by fitting the estimated frequency response function from the experiment. The rational fractional method proposed by Tanaskovic [16] is used in flutter model parameter identification. The rational fractional method is to represent the frequency response function as a rational fractional form, and then apply the linear least squares method to fit the estimated flutter frequency and damping. Although it is simple and easily implemented, as the traditional least squares algorithm fails to fully consider the influence of noise disturbance on the fitting result, it is difficult to accurately identify the model parameters when dealing with the noised test data, especially the damping parameters.
For this reason, the transfer function model is still used in this paper. The stochastic model of the flutter test experiment is shown in Figure 1.
In Figure 1, and are input and output signals, respectively, represents the artificial excitation applied to aircraft, and is atmospheric turbulence excitation. is one transfer function of the considered aircraft, and it is unknown and needed to be identified; is one time shift operator, i.e., . Since the flight test is inevitably affected by atmospheric turbulence excitation, then is regarded as an unmeasured excitation, and the random response generated by will be included as process noise in the measured response signal. is the flutter acceleration signal, and observed noises and are generated by the sensor. The processing method of the flutter test data and the choice of the excitation method are closely related with each other. At present, the commonly used excitation methods mainly include control-surface pulse excitation, small rocket excitation, frequency sweep excitation, and atmospheric turbulence excitation. The excitation of the flight test uses rocket excitation. The principle of excitation point and sensor arrangement is to effectively stimulate the first-order symmetric bending mode of the wing of interest, the first-order antisymmetric bending mode, and the first-order symmetric mode. And it is convenient to measure the response signal corresponding to these third-order modes. The flutter characteristics of an aircraft can be determined through estimating the frequency and damping of various flight states as a function of altitude and speed.
In time domain, the following relationships hold:
Substituting the expression into , we obtain that
After introducing one new observed noise , the expression of can be rewritten as
It means that the effect coming from unmeasured excitation is included as process noise in the measured response signal, so is neglected in the above equation.
Using the rational transfer function model for our considered aircraft structure, we havewhere and are coefficients of polynomial and and are orders of polynomial.
Define as the poles of the transfer function , where is the number of poles; then, it means thatand by using the correspondence relation between continuous time model and discrete time model, poles are obtained in continuous time domain, i.e., , where is the sampled period. So, the poles for continuous time model are given as
Then, the model frequency and damping coefficient can be solved as follows:
The main goal of aircraft flutter experiment is to identify the model frequency and damping coefficient . Obviously, accurate transfer function estimation is the premise of model parameter identification, so we give some assumptions as follows.
Assumption 1. All zeros of polynomial are outside the unit circle, and and have no common factor.
Assumption 2. Two observed noises and are two stochastic processes with uncorrelated zero-mean stationary Gaussian distribution, and their corresponding power spectrums are and . These two observed noises are independent of two noiseless variables and .
3. Analysis Process
The purpose of the identification problem in this paper is to determine the unknown parameter vector in the rational transfer function model. When given input-output observed datawith observed noises, then the considered aircraft flutter model parameters can be identified from equation (4).
Let the unknown parameter vector to be identified be
As the second-order statistic of the stochastic process can be used to describe the extent to which the dynamic system deviates from the equilibrium or working point, we introduce one parameter to describe the statistical characteristics of two observed noises. This introduced parameter corresponds to variance value. When observed output noise is colored noise, we set the parameter vector corresponding to observed noise aswhere denotes the variance value corresponding to the observed output noise at different time instant , and similarly, is the variance for observed input noise .
Furthermore, if observed output noise is white noise, then equation (10) reduces to
Combing equations (9) and (10), the unknown parameter vector to be identified in this paper can be obtained as
To formulate equations (1) and (4) as one more condensed form, we construct the following three regressor vectors:where is consisted by input-output data and is one regressor vector, consisted by input-output data with no noise. The elements of regressor vector are the observed noises. In the latter derivation process, is the true parameter vector and is its parameter estimation. and are two true polynomials for two polynomials and , respectively.
According to the above three regressor vectors, the following relation holds:
After a straightforward calculation,
Then, formulate equations (1) and (4) as the following linear regression form:
Based on this linear regression form (17), many classical identification methods are applied to identify the unknown parameter vector .
For convenience, set
For one stationary stochastic process , define its variance function as
Then, the covariance function between two stationary stochastic processes and are also defined as
Similarly, the variance function and covariance function for stationary stochastic process can be extended to variance matrix and covariance matrix for stochastic vector. E denotes the expected value. And in practice, variance function (matrix) and covariance function (matrix) can be estimated by observed data, i.e.,where is the number of observed sequence.
The above defined variance estimation is one of the main key technologies in the latter proposed identification method.
4. Instrumental Variable Variance Method
Set one new instrumental variable and consider the following parameterized system:
As the unknown parameter vector must satisfy the above equation, the classical instrumental variable estimation is denoted as , i.e.,whereand the above inverse matrix exists in case of Assumptions 1 and 2 hold.
From the result of equation (23), we see that the classical instrumental variable estimation is an unbiased estimation. The constructed instrumental variable must be independent of the random disturbance term , and it is also dependent on the regressor vector in order to make the matrix (23) well posed. Generally, the conditions that need to be met can be summarized as
Based on the related results from [6], the solution process of equation (22) depends on the following deterministic equation:
The choice of elements in the instrumental variable will determine the structure of and , and these two variance matrices are functions of the parameter vector .
By choosing , the obtained equation reduces to the bias-compensated least squares equation. Also, if the elements of are independent of the noise , the model parameter is named as the extended instrumental variable estimation. In order to obtain estimations of the parameter vector and , some elements of the instrumental variable must be selected at least in relation to noise .
Equation (26) can be regarded as one nonlinear least squares problem, where unknown parameter vectors and in equation (12) can be identified through minimizing the following criterion function:where the instrumental variable is chosen as
Regardless of whether the observed output noise is white noise or colored noise, the choice of the regressor vector must be satisfied:
The above equation shows that when observed output noise is one white noise, then regressor vector is consisted by the delayed input-output data. Furthermore, if observed output noise is one colored noise, then regressor vector is consisted only by the delayed input data.
To give our new instrumental variable variance method, we rewrite equation (26) as the following criterion function:where is one nonnegative weighting matrix with unknown parameter vector , i.e., . Our defined instrumental variable is chosen aswhere and are two selected variables by the researcher, and they satisfy
Using the defined three regressor vectors and random noises, we continue to abbreviate equation (26) as follows:where the new instrumental variable is chosen as
From equation (26), the new instrumental variable variance method can be changed as the following minimization problem:where the nonnegative weighting matrix depends on the unknown parameter vector .
For the criterion function in the minimization problem (35), when the parameter estimation is consistent, then it holds that
Linearizing the criterion function around the true parameter estimation , we havewhere
Neglecting the difference term , we haveand when parameter estimation approaches its corresponding true parameter value , it holds that
From equation (40), the gradient matrix of the criterion function with respect to the unknown parameter vector is obtained.
The solution to the optimization problem (35) can be solved iteratively with each other using the nonlinear separable least squares method proposed by the authors in [11]. Here, the simplest Newton method can be used to get a general understanding of the parameter estimation. This approximation is sufficient in practical applications. As the nonlinear separable least squares method needs to calculate an optimization solution for each iteration, this requirement may be solved several times, so it is not easy to be used in engineering practice, but only used in theoretical analysis.
Generally, the simplest Newton method is formulated here for solving that minimization problem (35). Step 1. Given the initial estimation Step 2. Compute and : Step 3. Let Step 4. Compute and set where is one forgetting factor, chosen by the researcher Step 5. Estimate and choose Step 6. Update Step 7. Repeat the above iterative process until the following inequity is satisfied: where is one arbitrary small scalar
5. Simulation Examples
To verify the feasibility and effectiveness of the proposed instrumental variable variance method, the transfer function identification of the current loop in flight simulation turntable and aircraft flutter model parameter identification are, respectively, used.(1)Firstly, the proposed instrumental variable variance method is applied to identify the transfer function of the current loop in flight simulation turntable, which is seen in Figure 2. Generally speaking, in identifying the transfer function of the current loop, we only consider white noise in output signal. However, in the actual test process, disturbances caused by natural wind, human factors, or signal acquisition instruments are inevitably mixed in the input-output signals, so the stochastic model in this paper needs to be used in the whole system identification process. In flight simulation turntable, the function of the current loop is to control the current of the motor not exceed the maximum stall current of the motor. At the same time, it is necessary to make the armature current strictly follow the change of the control voltage command, as the torque of the motor can be accurately controlled to eliminate the effect from the back electromotive force to the torque. The closed-loop structure of the DC motor current of the flight simulation turntable is shown in Figure 3, where the current loop includes pulse width modulation (PWM), electric motor, current regulator, low-pass filter, and current detector. And the combination of low pass filter and current detector is regarded as the feedback part. In Figure 3, and are the input and output in the whole current loop, respectively. In order to be able to utilize the identification method described in this paper, it is necessary to convert the above closed-loop system into an equivalent open-loop system through certain simplifications and calculations. The equivalence means that the outputs are same in case of the same input According to the actual test and calculation of the flight simulation turntable, the open-loop linear part can be obtained in Figure 4. In Figure 4, the first term is the filter, the second is PWM, and the third part corresponds to the electric motor. The electromechanical time constant in the motor is , and the electrical time constant is . Multiplying these three terms to obtain one rational transfer function form, i.e., In Figure 4, is regarded as the true input, is regarded as the true output, and then it is similar to Figure 1. Moreover, after two observed noises are considered, we obtain one linear system identification with two input-output noises in Figure 5. In Figure 5, as the rational transfer function form is unknown, the next step is to estimate it. Within the system identification theory, we only deem it as a true value to measure the identification accuracy. Assume the number of input-output data is 10000, and both observed noises are randomly independent and identically distributed random white noise with their variances ; furthermore, set the initial value as . Now, our proposed instrumental variable variance method is applied to identify each parameter in numerator and denominator polynomials from that rational transfer function. The transfer function consisted by identified parameters is named as the identified model. To verify the accuracy of the identification model, step input signal is used to excite the true model and identified model, respectively; then, the obtained step responses are in Figure 6. Also, the step response is given based on the classical least squares method in Figure 6, where the black curve denotes the true step response, the blue curve is step response based on our proposed method, and the green curve is step response in case of the classical least squares method. Through comparing these three curves in Figure 6, as our proposed instrumental variable variance method is one extended form on the basis of the classical least squares method, all the properties of the least squares identification method hold for our proposed instrumental variable variance method. The estimation accuracy based on our proposed method is greatly improved. As the influence of the two observed noises on the identification method and accuracy can be considered in our proposed method, the variance of the observed noise is also taken as part of the estimation parameters. So, our proposed instrumental variable variance method is more superior than the classical least squares method.(2)Secondly, the proposed instrumental variable variance method is applied to identify the aircraft flutter model parameter The flutter test data of a certain aircraft is used to verify the effectiveness of the proposed method. Aircraft flutter model is very complex, as it is full of flexible structure and aerodynamics. So, in this simulation part, we only use wind tunnel test to construct the flutter mathematical model of the two-dimensional wing. In the whole wind tunnel test for two-dimensional wing, the input is chosen as an artificially applied excitation signal, and the output is an accelerometer measurement, collected from many sampling data points. The simulation is based on 100 independent experiments and 500 data points. When comparing our proposed method with the classical instrumental variable method, the standards for measuring the performances are chosen as follows:(1)Computation time (S)(2)The estimated error of the model parameters on the test set is
True system is that
Define the input without noise aswhere is a white noise source, and the colored noise model is defined as follows:
Define the variances of white noises as and , and these two variances correspond to SNR of the input-output signal.where is the average power of the signal.
Then, we apply our proposed instrumental variable variance method and classical instrumental variable method on the same system until the optimal parameters are obtained, and then the identification results are shown in Tables 1 and 2, where, to simplify notation, our proposed instrumental variable variance method is simplified as IVC and classical instrumental variable method is simplified as IV.
From the two tables, we see as IVC is obtained on the basis of IV, its estimated performance result is better than IV, and the generalization performance is better too. But IVC needs to solve an optimization problem, so its computational complexity is obviously much more complicated than IV. In our recent advanced times, this more computational complexity is tolerable for us.
The transfer function's poles can be used to estimate the model frequency and damping coefficient. To improve the efficiency of our proposed method, we only calculate the pole of the transfer function during online flutter analysis, i.e., the poles of the transfer function is used to be a comparison criterion to verify the validity of our proposed method. From Figure 7, the true poles corresponding to the transfer function are very consistent with the identified poles from our proposed method. Then, the accurate model parameters can be obtained by using equation (7). Furthermore, to test the identification result for the molecular polynomial, comparison of Bode plots between the true system and identified model is given in Figure 8, where the approximation performance is acceptable.
(a)
(b)
The BACT wind tunnel model is a rigid rectangular wing with NACA 0012 airfoil section. The wing is mounted to a device called the Pitch and Plunge Apparatus, which is designed to permit motion in principally two modes-rotation (or pitch) and vertical (or plunge). The BACT system has dynamic behavior very similar to the classical two-degree-of-freedom problem in aeroelasticity. The preliminary analysis, control surface sizing, and flutter suppression control law design were based on the analytical state-space equations of motion of the BACT wing model. These equations were developed analytically, using structural dynamic analysis and unsteady doublet-lattice aerodynamics with rational polynomial approximations. These linear state-space equations consisted of 14 states, 2 inputs, and 7 outputs. This state-space equation is used for classical control law design and for performance simulation and verification purposes.
The analytical open-loop flutter dynamic pressure in air was 128 pounds per square feet at a flutter frequency of 4.5 Hz. The sinusoidal sweep signal is commonly used as input for system identification, and this sinusoidal signal is seen in Figure 9. Figure 10 shows the response of the wing trailing-edge and leading edge accelerometers due to 1 degree step input of the trailing-edge control surface in air at 225 per square feet dynamic pressure. Also from Figure 10, we can see that the primary plunge motion with small pitch diverges rapidly.
6. Conclusion
In this paper, aircraft flutter model parameters are identified based on one stochastic model of aircraft flutter. Through introducing instrumental variables to form the variance function, a certain norm of the difference between the theoretical value of the variance function and the actual estimated value is taken as a criterion function. The unknown parameter vector in the transfer function form is obtained by minimizing this criterion function, where this transfer function form corresponds to the aircraft flutter model. But we do not consider the convergence and accuracy of our proposed method, so these two aspects can be considered as our future work.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Natural Science Foundation of China (grant no. 60874037).