Dynamical Model Identification for a Small-Scale Unmanned Helicopter Using an Integrated Approach

This article presents an integrated approach for the parameter identification of a small-scale unmanned helicopter. With the flight experiment data collection and preprocessing, a hybrid identified algorithm combining the improved artificial bee colony algorithm and prediction error method is proposed to obtain the unknown dynamical parameters of the linear model. The proposed algorithm is valid to use thanks to an adaptive search equation, a novel probability-scaling method, and a chaotic operator and has a good performance in search speed and quality. Afterwards, we design a wind tunnel test to modify the main rotor time constant of the identified model. The identified accuracy and feasibility of the proposed approach are verified by making a time-domain comparison with three other algorithms. Results show that the dynamical characteristics of the helicopter can be determined accurately by the identified model. And the proposed approach is propitious to enhance the reliability and availability of the identified dynamical model.


Introduction
Small-scale unmanned helicopters play an essential role in civilian applications and military applications, such as environmental monitoring, emergency response, crime prevention, and urban mapping [1][2][3].The agility and high maneuverability of the helicopters require a high-performance control scheme for full-authority autonomous flight.And the design of a high-performance control scheme is based on a high-fidelity mathematical model that captures the key dynamics of the helicopters.However, obtaining the exact dynamical model of the helicopter system is not an easy mission due to its nonlinearities and unstable and high degree of coupling.
In general, there are two methods to the dynamical model establishment of the helicopter.Stating precisely, the first method is the establishment of a model using first principles.This method needs accurate instruments and a priori knowledge to obtain the parameters of the helicopter.The result of the first principles modelling is a set of high-order nonlinear differential equations that cover a wide portion of the flight envelope.Although several literatures [4,5] on the first principles modelling of the helicopters are reported, the technique is regarded as very time consuming and costly.Hence, these constraints render the first principles modelling method impractical for many applications.
The other method is the parameter identification technique.This is a procedure of deriving a mathematical model based on flight data including the control inputs and measured outputs.This method is relatively simpler and does not require lots of instruments and a priori knowledge of the system.Generally, parameter identification is classified into two categories: frequency-domain identification and time-domain identification [6].The frequency-domain identification uses frequency sweep experiments to compute the input-output frequency responses.With these responses, an estimated model can be obtained.Many works exist in the frequency-domain identification for helicopters.Metter et al. [7] and Raptis and Valavanis [8] use a software called Comprehensive Identification from FrEquency Responses (CIFER) to identify the linear model for the Yamaha R50 and X-Cell 60 helicopters.However, such an approach may be unsatisfactory because it needs the accreditation of CIFER and lots of frequency sweep experiments.Recently, there has been a surge of interest in the time-domain identification for helicopters.Kim et al. [9] uses the prediction error minimization (PEM) algorithm to estimate the unknown parameters of a Raptor E620 helicopter.However, PEM suffers from a high sensitivity of initialization.Combining a different evolution (DE) algorithm, a hybrid DE-PEM algorithm is developed to overcome the insufficiency of PEM and enhance its robustness [10].To a certain extent, the hybrid algorithm may be a viable option that can combine the advantages of subalgorithms to acquire a better performance.Moreover, some swarm intelligence algorithms have been employed in the time-domain identification for helicopters.Bian et al. [11] proposes a particle swarm optimization (PSO) to identify the linear model of a helicopter.According to the literature [12], the entire model is divided into two subsystems and the genetic algorithm (GA) is introduced to determine the parameters.While conducting the identification with large-scale parameters, the traditional warm intelligence algorithms may get stuck in a local optimum.To avoid the premature convergence of the GA, Wang et al. [13] apply a chaotic operator to improve its ability of information exploitation and identified the dynamical model which is successfully used for controller design.Therefore, the aforementioned studies make a major contribution to research on the parameter identification for the helicopter by displaying the possible methods.
In our study, a hybrid algorithm combining the PEM with an improved artificial bee colony (IABC) algorithm is proposed to identify the unknown dynamical parameters of a Trex-600 small-scale unmanned helicopter.In other words, the incipient scopes of the unknown parameters are firstly determined by the PEM.Furthermore, the more accurate parameters are found by the IABC.To overcome the premature convergence of the IABC, an adaptive search equation, a novel probability-scaling method, and a chaotic operator are adopted.The application of traditional ABC in parameter identification for an industrial robot has been explored and interpreted in our previous works [14,15].Plainly, the purpose of this study on parameter identification of a helicopter is to develop generic methods for achieving a high-accuracy dynamical model by the PEM-IABC algorithm.
A good aerodynamic performance can guarantee the satisfactory flight qualities of the helicopters.For modelling of the helicopters, previous published studies are limited to an ideal dynamical model without aerodynamic characteristics.To obtain these aerodynamic characteristics, computational fluid dynamic (CFD) or the wind tunnel test are the useful methods.While the CFD has had a big development over the past decade, it is treated as a supplement and may never be a substitute for the wind tunnel test.See [16] for deeper insight of the wind tunnel test.Consequently, we will design a wind tunnel test to modify the main rotor time constant of the identified dynamical model.
The remainder of the paper is organized as follows.In Section 2, the nonlinear dynamical model of a small-scale unmanned helicopter is established, including its rigid body kinematics and dynamics, main rotor dynamics, and yaw dynamics.Subsequently, this nonlinear dynamical model is linearized under a trimmed flight condition.In Section 3, we will highlight the process of the PEM-IABC and give the methods of experimental data preprocess.Section 4 and Section 5 present the experimental result of parameter identification and model modification with the wind tunnel test.The conclusions are discussed in Section 6.

Dynamical Model
In general, the first assumption toward dynamic modeling of a small-scale unmanned helicopter is to treat it as a rigid body with six degrees of freedom (DOF) [2]. Figure 1 displays the linear velocities u, v, w , the angular rates p, q, r , the attitude angles ϕ, θ, ψ and the main rotor flapping angles a s , b s with respect to the body-fixed frame O B X B Y B Z B .The nonlinear rigid body equations of motion for the helicopter that describe the vehicle's translational and rotational motion about the three referenced axes are deduced by the Newton-Euler equations [17]: where F x , F y , F z represents the components of the force acted on the helicopter.Similarly, the sum of the external moments is denoted by τ x , τ y , τ z .And the variable I σσ is the inertial moment along the σσ axis.
Obviously, the above equations are described in the body-fixed frame.As we know, the real motion of the helicopter is defined by the orientation and the position of the body-fixed frame relative to the inertial frame [18].A rotation matrix using the Euler angles can express the transformation between the two frames.Denote by v I x , v I y , v I z T , the linear velocity vector with respect to the inertial frame, one obtains where s ⋅ and c ⋅ are the abbreviations for sin ⋅ and cos ⋅ , respectively.Furthermore, the relation between the Euler angles and the angular rates is given by The control input vector of the helicopter is defined as u = u lat , u lon , u col , u ped T .Here, u lat , u lon , u col , and u ped are lateral cyclic, longitudinal cyclic, collective control, and tail rotor control inputs, respectively.The simplified main rotor dynamics can be derived from a couple of first-order flapping dynamic equations [19]: The term τ s represents the main rotor time constant, and it is given by where Ω mr denotes the ultimate rotate speed of the main rotor.I β denotes the blade moment of inertia about the flapping hinge.C lα,mr0 is the lift curve slope, ρ is the air density, c mr is the blade chord length, and R mr is the rotor radius.
A lon and B lat are constants associating the stick commands with the blade's pitch angles.Finally, the main rotor cross-coupling terms A b and B a are where λ β is the frequency.
In the real control of a helicopter, a feedback yaw rate gyro installed on the helicopter is used to reduce the effect of the antitorque fluctuation on the yaw response.The yaw dynamics can be found by the following form [20]: where r fb is the additional rate feedback term and K r and K fb are the parameters to be identified.Now that all the dynamical models have been obtained, the equations of motion can be transformed into a state-space form: with the state vector where M contains the rotor constants, the system matrix F contains the stability derivatives, and the control matrix G contains the input derivatives.
As already mentioned, equations describing the helicopter motion are nonlinear differential equations.Parameter identification using this model is difficult and meaningless.A linearized state-space model of the perturbed states and control input is obtained by small disturbance theory under a trimmed flight condition (e.g., hovering/cruise) [9].Note that the angular rate of yaw can be considered as a derivative of the yaw angle in the trim operating condition.In order to reduce the computational complexity of identification, the helicopter dynamics are separated into two decoupled subsystems, namely, the lateral-longitudinal subsystem and yaw-heave subsystem.These subsystems are expressed by where In the matrices above, there are more than 30 unknown parameters to be identified.Moreover, some parameters are unable to be directly identified from the actual flight data.According to the report [8], some assumptions can be obtained: T mr H mr /I yy , N fb = −N ped , and k fb = −2N r .Here, g is the gravitational acceleration, K β is the hub stiffness, T mr is the main rotor thrust, and H mr is the distance from the center of the helicopter to the center of the hub.

Parameter Identification
3.1.Identification Algorithm.PEM is the usual algorithm for parameter estimation by minimizing the prediction error.However, it is extremely sensitive to parameter initialization which affects it convergence to the global minimum.This disadvantage can be overcome by integrating the PEM into ABC optimization.The brief mechanism of the proposed identified algorithm is given below: Step 1.The incipient dynamical parameters of the helicopter are estimated by the PEM.And the detailed process of the PEM is reported in the literature [21]; it is unnecessary to repeat the work here.
Step 2. In the ABC algorithm, food sources x ij are D dimensional vectors and their values are randomly distributed between the lower parameter bound x L ij and the upper parameter bound x U ij .x ij are randomly generated by where i = 1, 2, … , N p and j = 1, 2, … , D. N p is the size of the food sources.rand ⋅ is a random function.
Step 3. New food sources are generated by employed bees.The size of the food sources is equal to the amount of employed bees.A new food source v ij is close to the current food source x T ij , which is generated by an adaptive search equation: where r ∈ 1, 2⋯,N p .Φ is a random scale factor between -1 and 1. ω min and ω max are the minimum and maximum weight factors, respectively.T and T max are the current and maximum iterations, respectively.To ensure the ABC convergence to the global optimum, a selection of the new and current food source is based on the fitness value: where fit i is the fitness value.The objective function f i is described as where N is the data length, y i is the processed actual flight data, y i is the mean of y i , and y mi is the estimated data.
Step 4. The food source with optimal fitness value selected by the onlooker bees is according to a probability.In traditional ABC, the probability is proportional to the fitness value.After several iterations, the weak food sources with low quality will be eliminated.However, these food sources also contain useful information.Accordingly, a novel probability-scaling method is proposed to improve the population diversity.
Step 5.If the fitness value of a food source cannot be improved during a predetermined number of cycles called limit, then this food source will be abandoned.And the relevant employed bees become the scout bees.In the scout bee phase, 4 International Journal of Aerospace Engineering a chaotic operator is introduced to prevent the ABC falling into local optimum.The chaotic sequence is generated by where y nj is the chaotic sequence between 0 and 1 except 0.25, 0.5, or 0.75.R is an amplification factor.x New i+1 j is the new food source.
Step 6.Output the best solution parameters achieved at the present time and go back to Step 3 until termination criteria T max is met.2, a modified Trex-600 unmanned helicopter is adopted as the experimental platform.It is equipped with a flight controller, an electronic aileron called KBAR, a couple of data links, a global positioning system (GPS), and so on.To determine the dynamical parameters of the model presented in equations (10) and (11), we propose an integrated method of parameter identification.The method involves three steps: direct measurement, parameter estimation using flight data, and model modification by the wind tunnel test.The third step will be described in Section 5.
Direct measurement mainly focuses on these physical parameters such as mass of the helicopter, main rotor radius, rotor speed, air density, and inertial moment.These parameters are generally related to the platform geometry, experimental environment, and loading.For instance, it is assumed that the principal axes coincide with the axes of the body-fixed frame; the inertial moment of the helicopter is measured via the three wire pendulum method [22].As depicted in Figure 3, the helicopter is suspended by three steel wires with equal length.All the directly measured physical parameters of the Trex-600 are listed in Table 1.
In the flight experiment, the helicopter is controlled to fly in hover via the ground station.Then, the pilot excites the roll and pitch channels by the Futaba to obtain the flight data in the lateral-longitudinal subsystem, and the altitude and yaw channels are kept in a small range to weaken the coupling effect from the yaw-heave subsystem.Finally, the flight data of the yaw-heave subsystem is collected only by exciting the yaw and altitude channels.The raw flight data consisting of the input and output are all sampled at the 50 Hz frequency and transmitted to the ground station.To ensure that the sample size is enough, the flight experiment is repeated several times.
Since the raw flight data contains some measured noised and outliers, it is immediately unfit for the parameter identification.In order to remove the outliers and attenuate the effect of the noises, a five-spot triple smoothing method is introduced to process the raw flight data [23].The method is given as where k = 3, … , m − 2, n i is the raw flight data, and n i is the processed flight data.The more equation ( 19) is used, the smoother the processed curves get.Obviously, excessive use of the smoothing method will lead to the rise of estimated error.A processed curve after using the method 5 times is shown in Figure 4.

Experiment Results
With the processed flight data, the unknown dynamical parameters of the Trex-600 helicopter is estimated by the proposed PEM-IABC algorithm.The parameter identification is conducted in MATLAB 2012b programming environment on an Intel Core i7-3770 PC running Windows 7. To avoid the influence of the inherent search randomness on the identified result, the two subsystems of equations ( 10) 6 International Journal of Aerospace Engineering and ( 11) are repeatedly identified ten times under the same condition.Moreover, three other algorithms are taken for comparison: the genetic algorithm (GA), traditional artificial bee colony algorithm (ABC), and chaotic artificial bee colony algorithm (CABC) proposed in the literature [24].It should be noted that all the proposed algorithms will identify the dynamical model based on the incipient parameters determined by the PEM.The parameters of the algorithms are listed in Table 2.In GA, the parameters p m and p c are the cross factor and mutation factor, respectively.
Figures 5(a) and 5(b) show the evolution curves of GA, ABC, CABC, and PEM-IABC regarding equation (13).Taking the yaw-heave subsystem as an example, PEM-IABC begins to converge in the third iteration, faster than the other algorithms.And PEM-IABC achieves a better result with a higher objective value of 0.8957.In this regard, the proposed algorithm has a stronger ability for information exploitation and a faster rate of evolution than the three other algorithms.
The dynamical parameters of the Trex-600 helicopter estimated by PEM-IABC are listed in Table 3.
With the estimated parameters and raw input data, the relevant predicted output can be calculated.The responses of attitudes, the three-axis linear velocities, and angular velocities from the four different identified models are shown in Figures 6, 7, and 8.It is obvious that the four algorithms can all get a model that matches the trends of the actual data.Nonetheless, the predicted output generated by PEM-IABC approximates the actual data best.From the results, the conspicuous error appears in the crest and the trough of a curve, since the algorithm may be trapped into the local optimum.For p response, the GA, ABC, and CABC trend to the local optimum from 1 second to 3 seconds.With the adaptive search equation, novel probability-scaling method, and chaotic operator, the PEM-IABC has the ability to jump out of the local optimum.
To make the comparison clearer, take correlation coefficient [13] in the identified result for accuracy analysis.As shown in Table 4, the closer the correlation coefficient is to unity, the better the identified model is, while if the coefficient is close to zero, the identified model is poor.Taking the pitch angle as an example, the correlation coefficient of PEM-IABC is 14.43%, 8.79%, and 4.02% higher than those of GA, ABC, and CABC, respectively.Seeing the results, PEM-IABC can guarantee a high-accuracy model, where the accuracy issue is of utmost priority.
Remark 1.In our previous work [25], CABC had already been proposed for the parameter identification of the small-scale unmanned helicopter.We have discovered that the population diversity along with data mining ability of the CABC would decline during the iteration.In particular, while dealing with the large scale of flight experimental data, the CABC will become sluggish and may be stuck in the local optimum.Hence, an adaptive search equation and a novel probability-scaling method are introduced to enhance the performance of the CABC, namely, the PEM-IABC.From the results in Figures 6-8, it can be observed that the identification precision of the PEM-IABC is higher than the CABC.Similarly, the correlation coefficient in Table 4 has also confirmed the conclusion.

Model Modification with the Wind Tunnel Test
In this section, the identified dynamical model of the helicopter will be modified through the wind tunnel test.The main rotor time constant τ s is related to the curve slope C lα,mr0 regarding equation (5).To enhance the model accuracy, the   To obtain the curve slope of the main rotor, the wind tunnel tests are divided into three steps, namely, the blow air test of the bracket, the blow air test of the full helicopter, and the blow air test of the main rotor.In the schedule, the test would be conducted when the wind speed is 0 or 10 m/s and the pitch angle and roll angle of the helicopter which is installed in the bracket are changed with the step of 2 degrees.The specific test steps are given as follows: Step 1 (blow air test of the bracket).As shown in Figure 10, the bracket without the helicopter is blown in the wind tunnel when the speed is 0 or 10 m/s.The motion of the bracket is given in Table 5.All the test data are collected by the control center.
Step 2 (blow air test of the full helicopter).As shown in Figure 11, the full helicopter is blown in the wind tunnel when the speed is 0 or 10 m/s.Meanwhile, the rotate speed of the main rotor is set at 1000 rpm, 1400 rpm, and 1800 rpm.The motion of the bracket is given in Table 6.Hence, 12 groups of experimental data are collected.It should be noted that the balance must be calibrated before each test.
Step 3 (blow air test of the main rotor).Since the tail rotor has a little influence on the aerodynamic performance, the curve slope of the main rotor can be measured only when blowing the helicopter without the tail rotor.As shown in Figure 12, the main rotor is blown in the wind tunnel when the speed is 0 or 10 m/s.The motion in this step is the same as that in Step 2.
The aerodynamic data of the main rotor are collected and processed at a frequency of 1000 Hz via the TQD-1 processor and DJS-131 PC machine in real time.The result of the curve slope C lα,mr0 is shown in Figure 11.With three different speeds of the main rotor, the curve slope has a linear relationship with the pitch angle.When the helicopter flies in hover in the flight experiment, the speed of the main rotor measured by the tachometer is 1397 rpm and the wind speed measured by the anemometer is 9.82 m/s.Also, the pitch angle and roll angle are approximately equal to zero at this time.According to Figure 13, the value of the curve slope C lα,mr0 is about -0.7691 rad -1 under the circumstances.The τ s is calculated as -0.3085 according to equation ( 5), and the estimation of τ s is -0.3361.As a result, the value of the main rotor time constant in equation ( 10) is rewritten as -0.3085.Note that the error between them is about 8.22%.

Conclusion
Dynamical model identification of a small-scale unmanned helicopter based on an integrated approach is given in this paper.The integrated approach roughly has two parts: the parameter identification process and the wind tunnel test.The parameter identification uses a hybrid algorithm PEM-IABC to optimize the unknown dynamical parameters of the Trex-600 small-scale unmanned helicopter.The wind tunnel test is conducted to modify the main rotor time constant of the identified model.Some conclusions are summarized in this paper.One is that the proposed PEM-IABC has higher identified accuracy when compared with GA, ABC, and CABC.With the three improved strategies, the robustness and search efficiency of PEM-IABC have been enhanced.Another is that the process of the wind tunnel test for the helicopter is presented.Based on the experimental result, it proves the credibility of the proposed parameter identification.In the future, we will continue to study the design of the control law based on the identified dynamical model.

Data Availability
The figures and tables used to support the findings of this study are available from the corresponding author upon request.
τ x = I xx p − I yy − I zz qr, τ y = I yy q − I zz − I yy pr, τ z = I zz r − I xx − I yy pq, 1

Figure 1 :=
Figure 1: Illustration for the helicopter and coordinates.

Figure 5 :
Figure 5: Evolving curves of the four algorithms.

Figure 8 :
Figure 8: Responses of the angular velocities.

Figure 9 :
Figure 9: Flow chart of the wind tunnel test.

Figure 11 :
Figure 11: Blow air test of the full helicopter.

Table 6 :
The experimental scheme of the bracket experiment.Wind speed (m/s) ϕ (degree) θ (

Table 3 :
The values of the estimated parameters.

Table 4 :
The comparison of the correlation coefficient.