Abstract
This article focuses on the problem of parameter estimation of the uncoupled, linear, short-period aerodynamic derivatives of a “Twin Squirrel” helicopter in level flight and constant speed. A flight test campaign is described with respect to maneuver specification, flight test instrumentation, and experimental data collection used to estimate the aerodynamic derivatives. The identification problem is solved in the time domain using the output-error approach, with a combination of Genetic Algorithm (GA) and Levenberg-Marquardt optimization algorithms. The advantages of this hybrid GA and gradient-search methodology in helicopter system identification are discussed.
1. Introduction
System identification techniques applied to aerodynamic parameter estimation of fixed and rotary wing aircrafts are well developed and commonly used in many research centers and aeronautical industries around the world [1, 2].
More specifically for helicopter parameter estimation, Hamel and Kaletka [3] presented a general vision of the progress in this field up to 1997 and Padfield [4] described a comprehensive flight dynamic theoretical model, flight qualities criteria development, flight test techniques, and several results of this research in the United Kingdom.
More recently, Tischler and Remple [5] described the state of the art in rotorcraft dynamic system identification in their book Aircraft and Rotorcraft System Identification. In this reference work, emphasis is given on engineering methods and interpretations of rotorcraft flight test results using a system identification software package in the frequency domain, known as CIFER (Comprehensive Identification from Frequency Responses). The frequency response methodology is in reality a hybrid method, since accurate mathematical modeling of the aircraft begins with flight test data collected in the time domain and the reliability of the identified model is also validated by a time-domain verification method using different flight test data.
It is important to notice that all above related works use local optimization algorithms based on gradient search methods such as Gauss-Newton and Levenberg-Marquadt search methods for finding a local minimum of the prediction error function at the system output. Concerning the use of global optimization algorithms and, more specifically, Genetic Algorithms they have been used by Hajela and Lee [6] in rotor blade design, by Wells et al. [7] in the acoustic level reduction rotor design, and by Zaal et al. [8] in the parameters estimation of multichannel pilot models, among others.
Regarding helicopter system identification techniques, very few articles have used Genetic Algorithms for global optimization of a cost function based on the prediction error. In this framework one can cite Cruz et al. [9, 10] in the longitudinal mode system identification of the Twin Squirrel helicopter and del Cerro et al. [11] in the identification of a small unmanned helicopter model. It should be emphasized that previous works [9, 10] address the short-period dynamics, considering as input data the longitudinal body axis velocity component, while in this work one deals specifically with the longitudinal short-period system identification of an helicopter, using the time-domain output-error approach combined with genetic and Levenberg-Marquardt optimization algorithms.
The rotorcraft identification methodology used in this work is the well-known (Quad-MV) methodology, proposed by Jategaonkar [12] and shown schematically in Figure 1. This methodology takes into account the main elements of rotorcraft system identification, including the rotorcraft excitation maneuvers, the aerodynamic data measurements, the mathematical model of the helicopter equations of motion, and the parameter estimation methods used to minimize the predicted output-error between the model and the real data. Finally, the identified rotorcraft model validation is done using new data from complimentary flight test maneuvers. Each one of these elements will be discussed in the following sections in order to clarify the proposed system identification methodology.
2. Maneuvers and Validation
The choice of the proper flight test maneuvers, by shaping the excitation signals, is of great importance to minimize the uncertainties in the parameter estimation procedures and to maximize the flight test data content. The optimization of the excitation signal is realized from the a priori knowledge of the dynamic content of the system that described the parameters of interest [13]. These optimization procedures, however, are of difficult application in the AS355-F2 (Twin Squirrel) helicopter, since no prior studies were available for its aerodynamic modeling. In this way, the maneuvers had to be specified by applying conventional flight test procedures and taking into account specific flight safety constraints for this class of rotorcraft.
A special sequence of sharp-edged pulses known as the 3-2-1-1 [12] was applied to the longitudinal cyclic inputs to excite the forward short-period flight dynamics with fixed collective control at an indicated airspeed of 80 kts and 5,000 ft of altitude pressure.
The validation procedure was done with a second 3-2-1-1 test sequence and also of a sinusoidal frequency sweep of the longitudinal cyclic inputs at the same airspeed and altitude, as described in more detail in Section 6.
3. Flight Test Measurements
The tested helicopter was equipped with the Aydin Vector Data Acquisition System (AVDAS), ATD-800 digital recorder, and a flight-test air data system, mounted on a nose boom, as depicted in Figure 2. This system measures a total of thirty-five different parameters with sampling rates varying from 18 to 72 Hz.
Some of the measured data channels include fuel quantity in each tank, nose boom static and dynamic pressures, external stagnation temperature, aerodynamic angle of attack () and sideslip (), roll, pitch, and yaw rates (, , and , resp.), load factors, longitudinal () and lateral () body attitudes, heading, collective, longitudinal and lateral cyclic, and pedal command deflections (, , , and , resp.).
The Earth axis speeds (, , and ) are obtained with the aid of a Z12 Differential Global Positioning System (DGPS), from Astech, whose antenna is fixed in the top of the vertical fin. The DGPS and AVDAS data synchronization is made by inserting a simultaneous event in both systems. The DGPS data is represented with the same AVDAS data sampling rate by means of linear interpolation procedure.
The wind direction and intensity are obtained comparing the body axis speeds with the aerodynamic speed from the flight-test air data system, mounted on a nose boom, at trim conditions. Consequently, the body axis speeds (, , and ) are easily calculated adding wind vector to the Earth axis speeds.
Thus, the observation vector () is given by (3.1), where is the state vector and is a constant vector included in order to compute offsets in steady state due to sensor misalignments errors and other interfering influences:
4. Helicopter Longitudinal Model
The helicopter equations of motion are deduced from direct application of Newton’s law for translational and rotational movements, as given by Prouty [14] and Cooke and Fitzpatrick [15]. This work deals only with the estimation of the aerodynamic parameters of the longitudinal motion. In this particular case, the longitudinal dynamic equations are decoupled and given as [14] where and represent the longitudinal and vertical forces, respectively; is the pitch moment and corresponds to the moment of inertia of a rotating body. The kinematic relation for the pitch rate about the -axis is written as
These equations of motion are clearly nonlinear, but a meaningful analysis can be realized by converting them into linear differential equations, by considering only small perturbations about a trimmed equilibrium point (represented by subscript 0) in the rotorcraft flight envelope. In matrix notation, the linearized dynamic equations are given by
Considering a pure longitudinal cyclic control input, with collective control fixed, and assuming that the short-period mode takes place at an approximately constant speed, the above equations are simplified to: which can be written as a state space equation: where is the constant and unknown bias vector that was included in the state-space equation in order to provide a first-order correction for the random errors due to nonmeasured inputs such as turbulences and noise, as proposed by Tischler and Remple [5].
In the state equation above, the time-delay can be associated with unmodeled dynamics, such as actuator dynamics, control linkages, and transient rotor dynamics, for instance.
Therefore, taking into account , there are a total of thirteen parameters to be estimated, as collected in the following vector of unknown parameters:
5. Parameter Estimation Method
5.1. Output Error Approach
The output error method was applied to estimate the aerodynamic parameters and the time delay between the cyclic inputs and the helicopter response. As shown in Figure 1, the basic principle of this methodology is to minimize a quadratic cost function related to the error between the in-flight measured responses and the simulated responses obtained with the identified mathematical model submitted to the same measured inputs. The cost function to be minimized is a function of the parameters of the dynamic model, such as the helicopter aerodynamic stability and control derivatives, sensor bias, and sensitivities.
Therefore, a cost function is defined which measures the matching of the real and simulated data, when a certain group of aerodynamic parameters (e.g., the coefficients of the proposed dynamic model) is varied. An iterative search method is employed to vary this group of parameters such as to minimize the adopted cost function, usually taken as a quadratic function of the predicted matching error.
Let be the cost function, and consider as being the vector of parameters to be estimated with M elements, the input vector of the system, y the vector of observations (sensor measurements), and let , where , be the simulated outputs of the proposed model. Then, the problem becomes how to determine the vector of adjustable parameters such that it minimizes the cost function given by .
There are many methods to calculate the cost function, including the maximum likelihood (MLE) and the least-square method. The least-square optimization procedure is basically a minimization of sum of the squared error at each measured point, given by
On the other hand, the MLE method takes into account the measurement noise variance to weight the system outputs during the optimization procedure through the measurement noise covariance matrix.
Suppose that each output value, , has an associated random measurement error with a normal distribution around the “true” value . Then, the conditional probability density function is given by where is the covariance matrix given as
From the equation of conditional probability density function given by (5.2), there can be obtained a modified cost function, applying the negative of the natural log on both sides:
Equation (5.4) shows that the parameter identification problem becomes a quadratic function minimization problem, since the maximization of the likelihood function is equivalent to the minimization of a weighted least-square equation.
Taking into account that the covariance matrix is constant, the cost function can be simplified to the following form:
The determination of a parameter vector that minimizes the cost function given by (5.5) is equivalent to finding a parameter vector that maximizes the probability of measuring, , .
Comparing (5.1) and (5.5), it can be concluded that the least squares cost function, except for the multiplication factor, is given by (5.5), replacing by the identity matrix. Therefore, in order to implement the weighted least-square cost function, it is enough to replace by a diagonal matrix to consider different weights to each state. In this work, a Matlab program was developed to minimize the cost function defined by both weighted least square and MLE errors.
5.2. Genetic Algorithm Optimization
Classic methods to solve the maximum likelihood minimization problem are Levenberg-Marquardt, Gauss-Newton, and Newton-Raphson, among others. This article proposes a hybrid search-gradient approach by using the Matlab Genetic Algorithm (GA) and the Direct Search Toolbox in order to estimate the helicopter longitudinal short-period derivatives combined with the Levenberg-Marquardt optimization algorithm, as shown schematically in Figure 3.
In case there is no a priori information available, this hybrid search-gradient methodology states that there should be applied initially the genetic algorithm and after that a classic method to solve (5.5). As large initial errors may lead to parameter estimation divergence, this method is especially important since the parameters estimation procedure becomes less sensitive to these initial values. Besides, another interesting characteristic is that the solution obtained by the genetic algorithm is globally optimal when the nonlinear model is used.
Based on biological evolution principles, the genetic algorithms are structured random search techniques for optimization. Therefore, they utilize concepts such as population size, individual characterization and processes related to selection, reproduction, crossover, and mutation of the individuals.
The GA optimization process starts by creating a matrix, representing an initial population IP, with each line formed by randomly chosen individual with fixed-length character strings (), as shown in the following: where represents a thirteen-component vector as shown in (4.6). This work uses two hundred individuals at each generation due to adequate computation time and convergence of the fitness function.
The following step is to assign a fitness value to each individual in the population using a fitness measure, represented by the cost function of (5.5). A new population is created by applying three genetic operations, depicted in Figure 4, to each individual string of the current generation. The individuals generated by the genetic processes are as follows
The GA iteratively performs the above steps until a termination criterion is satisfied. In this work, the maximum number of generations was used for this purpose.
Concerning the GA implementation in the Matlab environment, each current generation is composed of 2 elite individuals, with the rest of the population formed by crossover and mutation processes in the proportion of 80% and 20%, respectively.
The functions that create the initial population and produce mutation children were developed from a uniform distribution of random numbers, generated on a specified interval: ±100% of the BO-105 derivatives presented by Heffley et al. [16].
The output function is evaluated over a subdomain of individuals that satisfies the limits associated with bounds and stability constraints. While the bounds constraints refer to the above interval of the BO-105 derivatives, the real parts of the eigenvalues of (stability matrix presented on (4.5)), must be negatives (stability constraint) since the short-period dynamics is known to be stable on this helicopter.
5.3. Levenberg-Marquardt Optimization Algorithm
After genetic search verification (with a certain number of generations taken so that the fitness function tends to converge), the Levenberg-Marquardt optimization algorithm as presented by Mulder et al. [17] and Press et al. [18] is initialized. Therefore, the cost function given by (5.5) continues to be minimized iteratively according to where represents the estimate vector at th iteration, is the scalar chosen in accordance with the algorithm proposed by Press et al. [18] for a reduced value of , and is based on information about the cost function acquired at previous iterations. Numerically, is often obtained using the values of the first- and the second-order gradients of f() with respect to parameters vector , according to
A difficulty with (5.8) is that Hessian matrix may not be positive definite and thus not point in a “downhill” direction. Mulder et al. [17] propose to replace the by its expectation , known as the Fisher information matrix, that can be shown to be nonnegative definite. Therefore, (5.7) becomes
The first-order gradient of the likelihood function in (5.9) can be derived from (5.5) as
Applying the following equivalency presented by Goodwin and Payne [19] the Fisher information matrix may be obtained as
The Levenberg-Marquardt algorithm is an iterative technique that can be thought of as a combination of steepest descent and the Gauss-Newton method. When the current solution is far from the correct one, the algorithm behaves like a steepest descent method: slow, but guaranteed to converge. When the current solution is close to the correct solution, it becomes a Gauss-Newton method. Therefore, the iteration starts with initial estimates of all unknown parameters and computes parameter updates according to (5.9) based on (5.10) and (5.12). Assuming initially a modest value for , as 0.001, on (5.9), the cost function is evaluated, and if its value is greater than or equal to the last value, then is increased by a factor of 10 and the program returns to evaluate the cost function until the attainment of a value lower than the last value. Also necessary is a condition for stopping; so all these procedures repeat until the cost function values reach at least five significant digits in accordance with the last value.
6. Parameter Estimation Results
The 3-2-1-1 excitation maneuver using longitudinal cyclic inputs with fixed collective input was implemented in the flight test, under a nominal flight operation point of 80 kt and 5,000 ft.
Figure 5 shows the measured deflection inputs and Figure 6 presents the measured and simulated responses of vertical body axis speed w, pitch rate , and horizontal attitude .
(a)
(b)
(c)
As shown in Figure 5, five 3-2-1-1 maneuvers were performed in sequence in order to ensure a persistent signal and consequently to increase the signal power density spectrum at low frequencies.
The correlation coefficient between measured () and simulated data (), defined as the normalized cross-covariance function , is given by (Bendat and Piersol [20]) can be used to estimate how well the estimated signals can reproduce the measured data. If the correlation coefficient is close to unity, one may conclude that the estimation algorithm can provide a good fit to the experimental data, but on the other hand, if the coefficient is close to 0, the estimation was poor.
The fitting index, represented by the normalized correlation coefficient, is depicted in Figure 7. In this figure, each group of parameters represents the fitting index of the estimated outputs (, , and ), obtained with three different types of maneuver at same airspeed and altitude: (a) 3-2-1-1 sequence with longitudinal cyclic inputs used to estimate the flight derivatives (b) second 3-2-1-1 sequence and (c) sinusoidal sequence at increasing frequencies. The (b) and (c) groups of fitting index correspond to complementary flight data used only to validate the estimated parameters.
Considering the assumptions and limitations of the theoretical helicopter dynamic model with 2 degrees of freedom (DOF), it can be concluded that the obtained results are satisfactory, reaching fitting coefficients for the observed response in excess of 80% for all of cases, despite the lack of a priori knowledge about the stability and control derivatives of the vehicle.
This conclusion is confirmed analyzing the time response of all parameters and the Cramér-Rao Lower Bound (CRLB). It expresses a lower bound on the variance of each estimated parameter and, as presented by Tischler and Remple [5], normalized CRLB lower than 20% suggests a satisfactory estimation.
As all flight derivatives identified are dimensional, Table 1 presents the mass and inertia parameters used.
The estimated aerodynamic derivatives of the AS355-F2 short-period longitudinal mode at 80 kt and 5,000 ft and the respective CRLB are given in Table 2.
The CRLB presented in Table 2 was calculated considering that the time delay , the bias vector , and the vector , also identified, are constants inserted in the model. Therefore their values are written in Table 3 without CRLB. Finally, it should be noticed that the time delay, as shown by (4.5), was identified as the shift of the outputs measurements that minimizes the cost function given by (5.5).
In order to show the advantages obtained by the hybrid methodology described in this work, the optimization results with and without the parameter estimation using genetic algorithm (GA) are compared in Table 4, considering, as initial estimate, the derivatives of stability and control of BO-105.
The hybrid method, as shown in Table 4, improves the results in terms of the average correlation function coefficient of the three states identified and substantially reduces the cost function.
Finally, Table 5 presents the same estimated aerodynamic derivatives showen at Table 2 and the respective CRLB using only Levenberg-Marquardt approach.
It should be observed that considering, as initial estimate, the derivatives of stability and control of BO-105, the Levenberg-Marquardt algorithm leads to inadequate CRLB values as shown in Table 5.
7. Concluding Remarks
From the methodology , it was possible to obtain the longitudinal short-period stability and control derivatives in forward flight of the AS355-F2 helicopter at 80 kt and 5,000 ft.
The results showed that the hybrid search-gradient methodology is useful for parameter estimation problems by the output-error method, especially when the level of information regarding the studied system is reduced. As stated at previous sections, the adopted methodology initially estimates the flight derivatives using the genetic optimization algorithm and then, when the cost function tends to converge, the parameters are improved using Levenberg-Marquardt algorithm. This approach leads to better results in terms of correlation function coefficients, cost function, and CRLB values and at the same time increases the possibility of parameter convergence especially for large errors in the initial parameter values.
Concerning the GA, the more relevant aspects are as follows: (i)easiness to implement parameter restrictions: this point is of fundamental importance in the rotorcraft parameter identification algorithms, since it avoids divergence problems and estimated parameters without physical meaning;(ii)capacity to deal with the parameter estimation problem with low a priori knowledge level: the genetic algorithm for the stability and control derivative estimation is a powerful tool as a first parameter estimation approach.
On the other hand, the Levenberg-Marquardt optimization algorithm guides the parameters to a precise local minimum and also requires less computational time.
However, the helicopter dynamics has nonlinear terms, with parameters that may vary abruptly with the alteration of the flight conditions, or with helicopter configurations (large flap hinge-offsets) which require a larger number of degrees of freedom in order to introduce the rotor dynamics. Thus, in spite of satisfactory results obtained with linear models for the longitudinal short-period, more complex models will be required to improve curve fitting and to add other dynamic modes.
This is especially relevant to obtain a full flight simulator model and to develop control system strategies that depend on a comprehensive knowledge of the rotor and fuselage dynamic interaction in the entire flight envelope (weight, CG, altitude and airspeed) of the rotorcraft.
Concerning flight test maneuvers, the need for an appropriate specification of the inputs is also evident, taking in to account a priori model knowledge as necessary to maximize the level of information contained in the flight test data, minimizing the uncertainties associated with the identification process. So, the observed results show that the 3-2-1-1 sequence was suitable to identify the longitudinal short-period mode.
For the next studies, the following points will be taken in to consideration: (i)a more complex model will be proposed to include both longitudinal and lateral-directional dynamic modes; (ii)application of other optimization techniques such as Levenberq-Marquadt in conjunction with the maximum likelihood estimation problem and the Extended Kalman Filter (EKF);(iii)application of optimization techniques to generate maneuvers in order to maximize the level of information contained in the flight test data and, consequently, to make the stability and control derivatives estimation more robust and precise.
Nomenclature
| : | State-space representation |
| : | Center of gravity |
| , : | External force (longitudinal and vertical), N |
| : | Gross weight, N |
| : | Moment of inertia about the Y-axis, kgm2 |
| : | Mass, kg |
| : | Pitching moment, Nm |
| : | Dimensional moment stability derivatives, rad/ms, 1/s, or rad/s2cm for , , and , respectively |
| : | Fuselage angular rates, rad/s |
| : | Body axis velocity components, m/s |
| : | Input vector |
| : | State vector |
| : | Bias vector |
| : | Output vector |
| : | Output bias vector |
| : | Dimensional vertical force stability derivatives, 1/s, m/srad, or m/s2cm for , , and , respectively |
| : | Angle of attack, rad |
| : | Angle of sideslip, rad |
| : | Longitudinal cyclic position, cm or % |
| : | Collective position, cm |
| : | Variation from initial trim condition |
| , : | Fuselage attitudes (pitch and roll), rad |
| : | Identification vector |
| : | Correlation function coefficient |
| : | Time delay, s |
| : | Initial trim condition |
| : | Simulated response |