Method for designing multi-input system identification signals using a compact time-frequency representation

A flight test campaign for system identification is a costly and time-consuming task. Models derived from wind tunnel experiments and CFD calculations must be validated and/or updated with flight data to match the real aircraft stability and control characteristics. Classical maneuvers for system identification are mostly one-surface-at-a-time inputs and need to be performed several times at each flight condition. Various methods for defining very rich multi-axis maneuvers, for instance based on multisine/sum of sines signals, already exist. A new design method based on the wavelet transform allowing the definition of multi-axis inputs in the time-frequency domain has been developed. The compact representation chosen allows the user to define fairly complex maneuvers with very few parameters. This method is demonstrated using simulated flight test data from a high-quality Airbus A320 dynamic model. System identification is then performed with this data, and the results show that aerodynamic parameters can still be accurately estimated from these fairly simple multi-axis maneuvers.


List of symbols
Continuous signal x HT , z HT Distance between horizontal tail neutral point and wing-body neutral point

Introduction
Aerodynamic models used for certification, performance evaluations, handling qualities evaluations, flight simulators, control law design, etc. must be validated and possibly updated with flight test data to match the real aircraft stability and control characteristics. For this purpose, an extensive system identification flight test campaign is usually performed, in which dedicated maneuvers are conducted at distinct pre-defined flight conditions. Maneuver design for aircraft parameter estimation is usually done using an a-priori model of the aircraft, typically derived from wind tunnel experiments and CFD calculations. It is common practice to create inputs that excite the aircraft at its expected eigenmodes. As described in [8], Marchand [14,15] and Plaetschke et al. [22] showed that evaluating the frequency response magnitude of the terms of each equation of the aircraft's linear system is one possibility to identify the regions of identifiability of each derivative. These regions lie in the vicinity of the natural frequencies of the aircraft's eigenmodes. However, a-priori aircraft models are subject to uncertainties, and therefore the maneuver design must also consider frequencies slightly above and below the expected eigenfrequencies.
Multistep inputs are classical maneuvers to excite the aircraft at its expected eigenmodes. The length of the steps chosen herein is such that the natural frequency of the excited mode lies in the center or upper third of the input spectrum [8]. Multistep inputs like doublets or 3-2-1-1 signals, as in Fig. 1, are easy to execute manually and are therefore widely used for system identification purposes. Figure 2 3-2-1-1 input signal doublet input signal shows the normalized power spectral density PSD norm ∕ t for the impulse, doublet and 3-2-1-1 inputs. To ease the comparison, different values of t are used for each of these signals: t = 2 s for the impulse, t = 1.44 s for the doublet, and t = 0.85 s for the 3-2-1-1. Additionally, the respective PSDs were normalized such that each of the signals has unit energy (1/Hz) and is divided by the respective reference step length t , leading to a easily comparable unitless density function. The values of the reference step lengths were chosen such that the latter two signals (doublet and 3-2-1-1) are roughly centered around the same frequency. It can be seen that the 3-2-1-1 has a wider spectrum compared to the other inputs and consequently provides a more robust excitation in the presence of variations of the frequencies of the aircraft modes, e.g. due to different flight conditions. Frequency sweeps are another type of system identification maneuvers. They allow the evaluation of a complete frequency response of the system to the input signal. Usually this type of input signals is applied to one control surface at time.
Using maneuvers that excite multiple control surfaces at the same time has a great potential to reduce flight test time and costs. Such multi-axis maneuvers were already investigated in the late 70's, where Ramachandran and Wells [26] investigated the identification of aerodynamic parameters for a light aircraft by applying inputs on rudder and aileron simultaneously to minimize the correlation between the model parameters during identification. Further multi-axis maneuvers have been designed to excite the aircraft in all six-degrees of freedom.
Morelli described a method for optimal input design using dynamic programming, which yields globally optimal multi-input square wave signals [18,21]. This means that the signals maximize the information content in the data for a fixed maneuver time. The optimum criteria used is based on the Fisher information matrix, which requires an a-priori model of the aircraft.
Morelli also introduced a method to create multi-axis input signals based on orthogonal optimized multisine waves that cover a broad range of frequencies [10,19,20]. Those inputs are mutually orthogonal in time and frequency and can be designed manually or optimized, for both uniform and nonuniform power spectra [12]. In those references, accurate parameter estimates using the equation error method in the frequency domain were obtained from flight test data with this type of input applied continuously during the maneuvers. Note, however, that data obtained using multisine inputs can be processed using all kinds of system identification methods, including the output and equation error in time domain; using the multisine excitation signal does not force the use of frequency-domain methods.
Another approach for multi-axis maneuver design is to use the method proposed by Lichota in [11,12], which uses genetic algorithms to optimize the power spectra for multi-axis multisine inputs or to define the switching times for multi-axis multi-step inputs, based on cost functions computed from the information matrix associated with a prior model. In these studies the prior model was linear but can be also nonlinear, and there is no time information about the frequencies excited.
The goal of this research is to propose new ways to design maneuvers, which can, furthermore, be used to reduce flight test time or improve the accuracy of airplane simulation models. These simulation models must at least fulfill the Qualification Test Guide (QTG) criteria for level-D full flight simulators [5]. With the increase in computational power, new Virtual Flight Testing (VFT) techniques based on coupled CFD/CSM simulations are becoming available and are expected to play an increasing role in the future. The present work is part of a DLR (German Aerospace Center) project named VicToria investigating the use of virtual flight testing in aircraft early design stages [4]. This technique will eventually become cheaper and a standard aircraft design tool, but it is still fairly expensive and time-consuming. As a consequence, maneuvers used in virtual flight testing must be kept as short as possible. Contrary to real flight test, virtual flight test permits assessing any physical quantity that is simulated (aircraft motion, pressure distribution, flow velocities, etc.) without sensor noise or calibration errors, which are ideal conditions for system identification. At the same time, virtual flight testing will always be as good as the simulation models used, and the need for real flight testing is not eliminated.
A method for designing multi-axis excitation maneuvers by specifying the frequency content and the times when the frequencies are excited using very few parameters was developed for this virtual flight testing application. Naturally, the method and the maneuvers designed with it can also be used for real flight tests. A parametrization is performed using the wavelet transform, which has been applied for many applications including image processing, seismic signal denoising and analysis of diverse other physical phenomena [2]. The wavelet transform yields a Time-Frequency Representation (TFR) of a signal, and the signal can be reconstructed from its TFR. The idea behind the proposed method is to start by specifying the desired TFR and to generate the input signals by inverse wavelet transform.
This work is structured as follows: the next section gives a short overview over the applied rigid-body equations of motion and the aerodynamic model. Section 3 introduces the wavelet transform. The method used for the signal definition is explained in Sect. 4 together with its applicability to multi-axis signals. In Sect. 5, parameter estimation results are presented that were obtained from simulated flight test data using the new signal generation method. A comparative discussion to already existing multi-axis maneuver designs is given in Sect. 6. Finally, Sect. 7 sums up the work in this paper giving the next steps and future applicability for the method introduced.

Basic model formulation
The basic aircraft equations of motion are described by a six degree-of-freedom dynamic model. The translational motion is given by and the rotational motion is given by The aerodynamic equations are described in the following by Eqs. 3 to 5. In these equations, colors are used to distinguish different groups of parameters. The parameters in black are included in both, the simulation model used to generate the data and the model identified from this data. The parameters in violet may either be useful for system identification of airplanes exhibiting some asymmetrical behavior ( C Y 0 , C n 0 , C l 0 ) or have been found useful in past work from the authors, but are neither in the simulation model nor required here to achieve a good match with the data (e.g. C Y r or C L e 0,HT ). Finally, the parameters in green are included in the simulation model used to generate the data, but are not required to achieve a good match with the data, and are therefore not included in the identified model (e.g. k, C L e q,WB and C Y p ). The determination of a suited model structure is an important step in the system identification process and has been investigated in many prior works, see for instance Reference [10] (Sect. 5.4-Model Structure Determination) and references therein. This paper focuses on the input signal generation and the selected model structure is not further discussed.
The coefficients of the aerodynamic forces and moments for the lateral-directional dynamics ( C Y e , C l , C n ) are derived by simple Taylor series expansion. (1) where p * and r * are the nondimensional roll and yaw rates. The reference lengths l ref,p and l ref,r are both set to the semispan b/2. The equations for the longitudinal motion are based on the two-point model described in [8,17]. This model separates the wing and horizontal tail influences and allows to account for the downwash lag effect from the wing to the horizontal tail.
The lift coefficient is separated into a wing-body component and a horizontal tail component. The drag coefficient is calculated via the simple polar equation. The pitching moment coefficient is also separated into a wingbody component and a horizontal tail component calculated via the body-axis components C X HT and C Z HT and the respective lever-arms.
The separated influences of wing and tail for the longitudinal motion are calculated by where q * is the nondimensional pitch rate. The reference length l ref,q is defined as the mean aerodynamic chord c 1 . The downwash angle HT describes the influence of the wing on the tail, where 0,HT is the downwash angle at the horizontal tail for = 0 • , calculated by 0,HT = , and is the time delay = x HT V to consider the transport lag of a change in the flow from the wing to the tail.
The body-axis force coefficients C X and C Z from Eq. 1 are finally derived from the lift and drag coefficient by Note that the intermediate system 'e', as defined in the ISO 1151-1 standard [1] and in which the lift and drag coefficients C L e and C D e are here expressed, is such that there is no rotation with the sideslip angle in the transformation between the intermediate system and the body system.

Wavelet transform
The wavelet transform is used for the analysis of diverse physical phenomena, e.g. in the denoising of seismic signals, climate analysis, heart monitoring, amongst others [2]. The complete wavelet transform theory goes far beyond the scope of this paper and will not be described in detail. A few essentials are provided to give the reader some insight into the transform applied in the method. More details can be found in [2,9,13].
Roughly speaking, the wavelet transform is a convolution of a signal to be analyzed with a wavelet function, also known as wavelet. The wavelet is a small wave-like function that begins and ends at zero amplitude. Some commonly used wavelets are depicted in Fig. 3. As will be shown in more detail later, for the present study several wavelets from the the well-known bior wavelet family [13] were used to generate the input signals for the primary control surfaces. This is to avoid sharp edges, which would add undesired high frequency content to the signal. The parameter A is the dilation parameter, which is used to scale (stretch or squeeze) the wavelet. The parameter B is used to translate (shift) the wavelet to various locations, see Fig. 4a.
The continuous wavelet transform of x(t) at any scale A and position B generates a two-dimensional transform plane as indicated in Fig. 4b. The x axis represents the location (e.g. time shift) of the wavelet function while the y-axis indicates the current scale of the wavelet. If the signal matches well with a wavelet at a certain position, a large value in the transform plane is expected.
Mallat [13] has shown that a multiresolution representation of the signal is achieved, applying the Discrete Wavelet Transform (DWT) using a certain set of parameters A and B. Therefore, a wavelet is characterized by wavelet filters g and h. The decomposition of a discrete signal is then obtained by its convolution with these wavelet filters, in a filter bank, using the following equations: where the index m (respectively m + 1 ) denotes the decomposition level.
The multiresolution representation transforms the signal in a combination of approximation a m [n] and detail d m [n] coefficients, see Fig. 5. The approximations contain the lowfrequency information and can be interpreted as a general trend of the signal. They are obtained using the filter coefficients g d , whereas the details contain the higher frequency information, obtained using the filter coefficients h d . Indices d and r stand for decomposition and reconstruction. Another method to analyze discrete signals is the Discrete Wavelet Packet Transform (DWPT). It is a generalization of the discrete wavelet transform, however in this case the signal details and the approximations are further decomposed at each level. The direct representation derived from the DWPT of a signal is in the so-called filter bank ordering. A better insight into the nature of the signal is given by the natural frequency ordering, which can easily be obtained from the filter bank ordering by swapping some of the subdivisions of the high frequency/detail part. Both ordering methods and the swapping operations required are extensively described in reference [9]. The DWPT allows for finer decomposition at higher frequencies creating a complete decomposition tree structure as represented in Fig Fig. 6b, is decomposed and represented by a TFP in natural frequency ordering. As expected, the frequency of the signal increases with increasing time (x axis).
The decomposition coefficients in a m [n] and d m [n] are represented in the TFPs using so-called Heisenberg boxes. Each of these boxes defines the center frequency and location of the scaled and shifted wavelet function used to decompose the signal. In Fig. 6b these Heisenberg boxes are represented by colored rectangles, where the color indicates the value of the coefficient (cf. colorbar). This representation is similar to a spectrogram calculated from the time signal using the Fourier transform, where the frequency spectrum of a signal is depicted as it varies with time.
Applying the Inverse Discrete Wavelet Packet Transform (IDWPT), a signal a m [n] can be reconstructed using an approximation a m+1 [n] and a detail d m+1 [n] portion and the respective recomposition filter coefficients g r and h r . The reconstructed sampled signal a 0 [n] can then be transformed back to a reconstructed continuous signal x(t) by using, for instance, a zero-order hold. For more details on this see [13]. In this work, two wavelet bases were used: a so-called bior3.1 and a bior3.3. The filter coefficients of these wavelets are provided in Eq. (10) for the bior3.1 wavelet and in Eq. (11) for the bior3.3 wavelet. In the method described in this paper, the signal is directly defined through the approximation and detail coefficient at some chosen level and the inverse transform is used to create the continuous time signal that is provided as input to the simulation. In the examples from this paper, the  In an application to CFDbased virtual flight testing (DLR project VicToria [4,24]) a cubic spline interpolation was used instead as it permitted to prevent numerical integration issues with the CFD solver.

Methodology
The proposed maneuver design method is based on the idea of creating input signals which have distinct frequency band excitations at predefined maneuver times. An overview of the method is given in Fig. 7.
Aircraft a-priori information can be used to select the desired frequency band excitations for the signal. Heisenberg boxes in the TFP are used to represent this information in the time-frequency domain. An inverse discrete wavelet packet transform using a selected wavelet yields the desired signal in the time domain, which can then be assigned to any of the control surfaces.
As will be shown in Sect. 5, this method also allows the defintion of a single maneuver that uses several control surfaces at the same time. Through proper signal definition (10) in the TFP, a time and frequency decorrelation between the different control inputs can be assured.

Single input case
As the DWPT uses a dyadic scale the time-frequency plane must be defined as a square matrix having 2 N t time segments and 2 N f frequency segments. Each combination of a time segment and a frequency segment is a Heisenberg box. A value can be assigned to any of these Heisenberg boxes and the absolute value represents the local energy for the reconstructed signal, schematically depicted in Fig. 8. Furthermore, due to the Inverse Discrete Wavelet Packet Transform (IDWPT), this value defines the amplitude of the input. An example of signal generation starting from values specified in a TFP definition is shown in Fig. 8. The upper diagram of Fig. 8 shows the TFP with 16 time segments and 16 frequency segments. The signal length was set to 30 seconds, resulting in Heisenberg boxes with a width of 2 s and a height of approximately 0.28 Hz. In this example, an input with frequency excitations in the 0.2844 ± 0.1422 Hz band at 4 seconds and around 0.8533 Hz between 18 and 20 seconds was specified by assigning values for the corresponding Heisenberg boxes. The generated signal is the direct result of the IDWPT and is shown in the center diagram of Fig. 8. The lower diagram shows the frequency content of the signal, evaluated by a Fast Fourier Transform (FFT). This information can be used to evaluate if the combination of any wavelet chosen for the IDWPT introduces undesired frequencies.

Multi-input case
The multi-input design method used in this work follows the same steps as for a single input described above. For each signal, a TFP is specified. To assess the correlation in time and frequency between different input signals, an overlap diagram was created. It is a superimposition between the TFPs of the selected signals and gives a visual representation of the information used to define the inputs. The overlap diagram can be applied to get input signals for any type of system, where the aim is to avoid concurrent frequency band excitations, and this evaluation does not require any system response for the design. Figure 9 shows the overlap diagram and the resulting time histories for the part of the selected example in this paper. It is clearly visible that there is no concurrent frequency band excitation between these two control surfaces, even though both control surfaces are deflected simultaneously.

Application to airplane identification maneuver design
The proposed formulation of input signals can, in principle, be applied to all kinds of systems. For the design of input signals aiming at performing airplane identification, various guidelines can be found in the extensive literature, see for instance Reference [8] (Chapter 2-Data Gathering) and Reference [10] (Chapter 8-Experiment Design). Note that whilst a quite different formulation of the input signals is proposed in this paper, the physical relationships between model parameters and the cost function (here maximum likelihood in output-error in the time domain) has not changed. Therefore all the preexisting guidelines are still valid. For instance, the regions of identifiability of the aerodynamic parameters described in [8,14,15,22] provide a good basis for choosing the wavelet scales in the framework Multi-input signal generation scheme of the proposed method. Overall, it also remains crucial to excite the aircraft at and around the frequencies of its natural modes.
Regardless of the frequencies that users of the method choose for their identification program, the aim of the method is to provide a compact and easy way to generate and describe such input signals. It is important to note here that the time-frequency-based formulation used in the proposed method relies on excitations of frequency bands rather than combinations of pure frequencies, as in the case of multisines. This forces the user to separate the excitation of overlapping frequency bands of different control surfaces over time. For instance, in the example shown in Fig. 9 the frequency band around 0.3 Hz is excited by the rudder at the beginning of the signal and later excited again by the ailerons. When considering several separate maneuvers simultaneously during the parameter estimation, the respective excitations of a particular frequency band could also be present in only a subset of those maneuvers.
In the multisine case it is possible and common practice to define orthogonal excitations in the form of interleaved frequencies for the different control surfaces. This means that the same "frequency bands" can be excited at the same time without impacting the quality of the parameter estimates. As a consequence, for the same quality of the parameter estimates, it should be expected that the length of the required multisine excitation signals will be shorter than for the proposed method.

Simulation and results
To illustrate the properties and practical application of the proposed method, simulations with a high-quality Airbus A320 dynamic model were performed in replacement of real flight test data. This model has been derived from an extensive flight test campaign during the DLR-internal project OPIAM (Online Parameter Identification for Integrated Aerodynamic Modeling). It complies with criteria for simulator validation and qualification, and includes compressibility effects, ground effects, and high-alpha characteristics [23].
The reference aircraft data used for the calculation of aerodynamic forces are given in Table 1.
Maneuvers using the new design method described above were first developed for the longitudinal motion using the elevator and horizontal tail as control inputs. Therefore only the aerodynamic parameters of the longitudinal motion could be estimated from a subsequent estimation.
After this proved to be successful, a multi-axis maneuver using elevator, horizontal tail, aileron, and rudder was designed to allow estimation of the full set of aerodynamic parameters. Compressibility, nonlinear lift behavior, and ground effect were neglected for all identification runs and the parameters were always assumed to be constant for the respective maneuver.

Longitudinal motion
The method described in Sect. 4 was used to define a maneuver to identify the aircraft's aerodynamic parameters for the longitudinal motion. Two control surfaces were used for the maneuver: the elevator and the horizontal tail.
Inputs for the elevator contain excitations with frequencies close to the expected short period motion as well as low-frequency excitations to obtain a phugoid response. The horizontal tail is being deflected simultaneously to estimate its characteristics and the downwash parameters. For the elevator, a bior3p3 wavelet, from the bior wavelet family as described in Fig. 3 and Eq. (11), was chosen to generate a signal with no sharp edges.
To counteract the deviation from the trim point, the horizontal tail is being deflected in the opposite way as the elevator. A variation of less than 1.5 degrees of angle of attack and of no more than 4% of true airspeed allows the estimation of parameters without any need to account for additional compressibility effects. Therefore, the aerodynamic parameters can be assumed constant throughout the complete maneuver at a given flight test point. Figure 10 shows the input signal generated for the elevator. The upper diagram shows the definition of the timefrequency plane. A 32 × 16 TFP was chosen resulting in Heisenberg boxes with a width of 2 s and height of approximately 0.2753 Hz. The center diagram shows the resulting signal from the inverse wavelet packet transform. The bottom diagram in Fig. 10 shows the frequency content of the final signal. It can be seen that no frequency above 2 Hz is excited. This is desired when creating signals for rigidbody identification because the elastic modes are generally Usually, some a-priori information for a new aircraft is available from wind tunnel results, CFD calculations, or after the preliminary design. With the proposed design method, the expected natural frequencies can easily be excited by using corresponding frequency bands. In the current evaluation, the expected short period natural frequency is 0.35 Hz. To account for this, using the 32 × 16 TFP tiling, three frequency bands, centered at approximately 0.2753, 0.5505 and 0.8258 Hz are included for the signal generation. For the phugoid mode, the expected natural frequency is much lower, on the order of 0.02 Hz. Therefore, a classical pulse is emulated using the lowest frequency band over a few time tilings. For display purposes, only relevant portions of the chosen time-frequency plane tilings are shown, meaning that all other coefficients are zero and do not influence the results of the IDWPT.
Parameter estimation was performed using the output error method in the time domain and the aerodynamic model for a rigid body aircraft for the longitudinal motion as described by Equations (4)- (5). For the simulated flight test data, Gaussian noise with a standard deviation according to Table 2 was added. The noise levels in Table 2 are representative of the levels that can typically be observed on flight test data of large transport airplanes. Note that some of these channels (e.g. ṗ , ̇q , ̇r ) are usually not directly measured in practice, but rather derived from the other quantities. In those cases, the noise levels used remain representative, but no correlation between their noise and the noise added to the quantities that they were derived from was modeled here; their noise was generated independently from each other here.
The comparison between the simulated flight data and the identified longitudinal motion model is shown in Fig. 11 and the results of the comparison of the time histories are given in Table 3. The maximum absolute difference between the identified model and the simulated flight test data is shown in the second column. The third column gives the relative difference to the overall output amplitude during the maneuver, and the last column shows the standard deviation of the output errors. The maximum absolute difference max , the relative difference rel,max , and the standard deviation are defined as follows: where z(t i ) are simulated flight test data outputs and y(t i ) are identified model responses.
Identified parameters for the longitudinal motion are shown in Table 4. The relative difference rel = estimate − true true between the estimated value for the parameter for the reduced model and its true value in the simulation model is given in the last column. The initial values used are also indicated in Table 4. To illustrate the good convergence of the estimation, these values were chosen reasonably far from the true values, even if it is usually easy to define better start values than the ones used here.
For simplicity, only one maneuver at one flight condition was used in this example. The results show that the lift parameters were obtained with the highest accuracy, however, the indetermination of the drag parameters arises. This yields less accurate estimates for C D e 0 and the Oswald factor e. The indetermination of C D e 0 and the Oswald factor e is a well-known phenomenon while  identifying the aerodynamic model with only data from one flight condition, and is similar to the one described in [6,17]. In the current example, only a small portion of the drag polar is covered during the maneuver, what can be seen depicted by the black markers in Fig. 12. In Fig. 12 the drag polar derived from the identification (red curve) is compared to the drag polar from the simulation model (blue curve). Note also that the drag polar of the simulation model contains the linear term k which is not included in the identified model. The small portion of the drag polar covered yields an additional indetermination if the linear term k is added to model structure. In practice, several flight conditions should be considered simultaneously, which removes this issue. These results show that it was possible to design a maneuver using the method in Sect. 4 that allows to successfully identify a nonlinear model with nine parameters that accurately describes the longitudinal motion for the rigid body aircraft, only using the a-priori information of desired frequencies and a proper wavelet function.

Multi-axis maneuver
To be able to estimate the aerodynamic parameters for a complete aircraft model, a multi-axis maneuver was   The input signal generation for the other control surfaces is shown for the rudder in Fig. 8 and for the aileron in Fig. 13. The signal designs were performed using a 16 × 16 time-frequency tiling and the bior3p1 wavelet from the bior wavelet family, see coefficients in Eq. (10). This results in Heisenberg boxes with a width of 2 s and a height of 0.2844 Hz. The Dutch roll of the aircraft is excited in the beginning of the maneuver by a rudder input with low frequency at 0.2844 Hz. Furthermore, the aileron is deflected excite the roll motion. This is the same approach as for the classical maneuver design criteria, described in [8]. Additionally, other frequency band inputs are applied to the aileron at the beginning of the maneuver and to the rudder in the second half of the maneuver, shown in Fig. 9. The signals generated for the lateral directional control surfaces, similar to the elevator, do not contain any frequencies above 2 Hz, thus avoiding significant excitation of elastic modes.
Parameter estimation was again performed using the output error method in the time domain and the aerodynamic model for a rigid body aircraft, as described in the complete set of Equations (3)-(5).
Overall, a multi-axis maneuver was successfully designed and allows the estimation of 21 aerodynamic parameters that describe the rigid body aircraft dynamics in six degrees of freedom. Plots for the longitudinal and lateral-directional motion of the aircraft can be seen in Fig. 14. Gaussian noise according to Table 2 was added to the measurement data to account for measurement noise as it would be expected in the real flight experiment. The differences between the identified model outputs and simulated flight test data are given in Table 5. Table 6 shows that all main parameters are estimated with good accuracy. It should be emphasized that, analogous to the longitudinal maneuver, the model used for identification is a simplified model with some selected parameters to describe the longitudinal and lateral-directional motion. The comparison plots in Fig. 14 show that a good match between simulated flight test data and identified model is obtained.
To evaluate the quality of the obtained model, its prediction capability is typically evaluated with a subset of the flight test data which was not used for the system  Table 6. The respective deviations between the real and identified values are not always easy to interpret in terms of prediction capability of the identified model. Additional simulations were performed with both the original and identified models and with a different multi-axis maneuver to compare the reactions of both models on this maneuver. Figures 15, 16 show the comparison of the response of both models. As the simulations are not used for identification but only for comparison, no measurement noise was included in either simulation. Both figures correspond to the same simulations but were split for readability reasons. Figure 15 shows only the first 25 s of the time responses, which permits the observation of the good match of the identified system on the short period, dutch roll, and roll modes. Figure 16 shows these responses until 150 s and shows that the phugoid is not perfectly identified. This results from the constraints from the virtual flight testing with coupled CFD/CSM environment, which restrict the maximum signal length for the identification, and from the fact that only one maneuver at one particular flight point was used. Even if the difference between both responses is observable in Fig. 16, in practice any pilot or control law would easily cope with such modeling errors. Note also that the phase and amplitude difference seen in Fig. 16 cause the error on the phugoid mode to appear larger than it is in reality. The difference between both simulations is mainly due to a difference in the strength of the initial phugoid excitation during the pitch maneuver (i.e. in the coupling from the short period motion to the phugoid mode) and is only marginally due to the phugoid dynamics itself. Overall, considering that only one maneuver was used to identify the model, the obtained match between the predicted response from the identified model and the response of the baseline simulation on this verification maneuvers is very satisfying.

Comparison with other maneuver design methods
In this section, a qualitative comparison between the proposed method and the main approaches from the literature is attempted. Note that, in the general case, such comparisons can hardly be made fairly. The quality of the estimated parameters during the identification process might indicate the superiority of one method on another one when considering a specific application and the opposite result could be obtained for another application. In addition to that, there is no single and commonly accepted metric to judge of the quality of input signals, not even for a single application. It is crucial to note that the "quality of input signals" usually involves a combination of desirable properties whose relative importance might be viewed very differently by different practitioners. The presented methodology was used for the definition of a multi-axis maneuver for a so-called CFD-based virtual flight test for system identification of flexible aircraft (DLR project VicToria [4,24]). The virtual flight testing using URANS with Fluid-Structure Interactions (FSI) is very time-consuming and therefore the defined maneuvers need to be efficient. In the case of the VicToria project, only a 10-second maneuver could be defined to serve as a testcase for the identification of the flexible aircraft stability and control characteristics in all axes. A further aspect, which will not be addressed in the following, concerns the selection of the operating/flight points at which the maneuvers are performed. Indeed, to obtain a model of the system/airplane that is valid across all relevant operating conditions, data covering the entire range of operating conditions (e.g. flight envelope for an aircraft) must be gathered and analyzed during the system identification process. For conciseness, this paper focuses on the design of maneuvers for gathering relevant data for a single operating point. This same constraint was also assumed for the multiaxis maneuver designed for the work published in [24].
Maneuver designs for parameter estimation must provide flight test data with rich information content so that an accurate model can be extracted from the system responses. A  maneuver with high information should include all frequencies as expected by analyzing the a-priori model to estimate the aerodynamic parameters accurately. Additionally, for input signals to enable an effective and precise identification of the system's internal dynamics and of the effectiveness of the control effectors, the excitation input signals must (among other things) be uncorrelated. Classically, separate excitations/maneuvers using only one control surface are performed, and the response of the aircraft to each of these input signals is recorded. This approach necessarily leads to uncorrelated inputs as only one input is used, but has two main drawbacks. First, the total duration of all these maneuvers is unnecessarily 3 large and therefore the overall system identification is more costly. Second, interaction effects between the control surfaces cannot be captured, unless even more maneuvers with simultaneous combined excitations are performed as well.
In order to reduce the total experiment time required for the system identification and/or to enhance the parameter estimate results, many researchers have proposed input signals and input signal generation methods aiming at exciting the system on several of/all its inputs at the same time, see for instance [11,16,19,25]. One crucial aspect when exciting the system through several inputs at the same time is that the input signals must be uncorrelated as this ensures that the effects of the respective inputs can be separated afterward.
One possibility is to use so-called multisine inputs, which are built by summing several pure sine signals for each input, however with disjointed sets of frequencies for each input. This strategy has been already successfully used for airplane system identification, e.g. by Morelli [19] and Grauer et al. [7]. When properly designed, orthogonal optimized multisines have been shown to be highly effective and efficient in flight tests.
A simple way to generate multi-axis maneuvers is to combine several classical single-axis maneuvers. Such maneuvers designs can be found in the literature, see for instance [11,25,26]. It remains, however, challenging to properly design "good" multi-axis maneuvers this way, as the classical single-axis maneuvers often involve fairly large frequency ranges. Consequently, they will often not be as well uncorrelated as desired. Besides, when the number of controls grows, it becomes rapidly harder to define such signals with a reasonably low level of correlation.
Using wavelets (i.e. time-bounded "wave-like" signals with zero mean), an alternative basis of the signal space is used instead of pure sine waves. Each wavelet corresponds to an excitation around a given time and within a given frequency band. This time-frequency representation of waveletbased multi-inputs signals can easily be visualized with the overlapping diagram described in 4.1.2, on which the signal designer can directly see whether the inputs are uncorrelated or not. In this regard, the proposed methodology reuses an idea which already underpins the multisine input signal generation as used in [19,20].

Conclusion and outlook
In this work, a new method to design single-or multi-input signals that can be used as maneuvers for aircraft parameter estimation was presented and discussed. Based on classical design criteria, the methodology developed allows the user to design complex signals both with and without a-priori information on the system to be identified (here the aircraft), even though only a quite restricted number of parameters is used to describe these signals.
The new method for maneuver design has shown promising results for both, single-axis and multi-axis excitations. The proposed method was demonstrated based on a 30-second maneuver which permitted to successfully estimate the aircraft's aerodynamic parameters for the longitudinal motion. It was also successfully demonstrated that it was possible to design a 30-second maneuver for which 21 parameters describing the complete aircraft rigid body dynamics for a single flight test point were accurately estimated.
A small number of design variables for the input signal definition, especially the choice of the type of wavelets and scale, seem very well suited for further automatic optimization of the maneuvers aiming at maximizing the quality of the identified parameters. Whilst such an optimization remains a medium-term goal at this stage of the development, it should be noticed that the parametrization of the experiment design used herein will ease the interpretation of the aerodynamic parameter values that would be obtained through the optimizer. On the other hand, the a-priori knowledge available about the system can easily be included in the design of the signal in order to be certain to excite the relevant dynamic modes of the system.
The proposed methodology also permits the design of signals, rich in information, very quickly, which can be fully or almost fully uncorrelated. Additionally, this design method can easily account for specific constraints, such as control surface deflection rates, amplitudes and excitation frequencies.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long 3 By comparison with the other approaches mentioned later. as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.