Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 4, pp. 1193-1203, Warsaw 2017
DOI: 10.15632/jtam-pl.55.4.1193

POWER SPECTRUM OPTIMIZATION IN THE DESIGN OF MULTISINE

MANOEUVRE FOR IDENTIFICATION PURPOSES

Piotr Lichota, Joanna Szulczyk

Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Poland

e-mail: plichota@meil.pw.edu.pl; jszulczyk@meil.pw.edu.pl

Daniel Agudelo Noreña, Felipe A. Vallejo Monsalve

Universidad de San Buenaventura, Bogotá, Colombia

e-mail: dagudelon@gmail.com; fvallejo@usbbog.edu.co

In this paper, two sets of multisine signals are designed for system identification purposes.
The first one is obtainedwithout any information about systemdynamics. In the second ca-
se, the a priori information is given in terms of dimensional stability and control derivatives.
Magnitude Bode plots are obtained to design the multisine power spectrum that is optimi-
zed afterwards. A genetic algorithm with linear ranking, uniform crossover and mutation
operator has been employed for that purpose. Both designedmanoeuvres are used to excite
the aircraft model, and then system identification is performed. The estimated parameters
are obtained by applying twomethods: Equation Error andOutput Error. The comparison
of both investigated cases in terms of accuracy andmanoeuvre time is presented afterwards.

Keywords:flightdynamics, optimal inputdesign,multisine excitations, systemidentification,
least squares, maximum likelihood

1. Introduction

System identification (Sys-ID) techniques play an important role in the development of high-
-quality aircraft-simulation models from flight-test data (Grauer, 2016; Albisser et al., 2017;
Viana, 2016). Themain purpose of those approaches is to develop amathematical model of the
aircraft. To achieve this aim, an adequate input has to be designed in order to obtain dynamic
characteristics of the object from flight tests.

The research on this topic started in the 1960s when Levin (1960) made the first systematic
attempt at obtaining optimal inputs. Later, in the late 60s, Nahi and Wallis (1969) made a
significant step by using an optimality criterion to design inputs in the time domain. Since then
the use of optimality criteria was a great deal of interest and was investigated by, e.g., Fedorov
(1972), Kalaba and Springarn (1982) or Goodwin and Payne (1997). At first, designed inputs
were used for manoeuvres in which only one flight control was deflected. Further investigation
made by, e.g., Wells and Ramachadran (1977) or Morelli and Klein (1990) showed that desi-
gning sets of excitations in which flight-control surfaces were deflected simultaneously was also
possible. Even though, none of the theoretical developments seemed to be fully adapted for so-
lving complex issues. Some recent progress has beenmade in applying numerical algorithms and
incorporating themwith other tools to design optimal inputs (Seren et al., 2006, 2013; Lichota,
2016; Lichota et al., 2017).

To shorten the time of the flight test campaign, multisine excitations that allow simultane-
ous flight control deflections can be used. Planning manoeuvres with those inputs requires the
designing of their power spectrum. In our study, a genetic algorithm has been used to optimize



1194 P. Lichota et al.

that energy distribution. Themain purpose is to improve the input effectiveness bymaximizing
the power stored in certain harmonics.
The organization of this article is as follows. First, we describe the process of designing

multisine signals for simultaneous aileron and elevator deflections. That section is divided into
two parts. In the first one, a method for designing the inputs without a priori knowledge of
the object is presented. In the second part, the optimization process (based on the available
information of the aircraft dynamics) that leads to obtaining energy optimized multisine exci-
tations is described. Parameter estimation methods used for estimating stability and control
derivatives from the model response are presented in the next Section. Comparison of inputs
in terms of manoeuvre time and Sys-ID accuracy is shown thereafter. The paper closes with a
short summary of conclusions.

2. Multisine input design

One of the approaches that allows one to shorten the time spent during flight tests in order to
reduce the total cost of the Sys-ID process is to design a manoeuvre with simultaneous flight
control deflections. Multisine input signals that are mutually orthogonal in the time and frequ-
ency domain can be easily designed for those purposes. To achieve this aim different harmonics
need to be assigned to each flight control (Morelli, 2009)

u=
M∑

k=1

Ak sin(2πfkt+φk) (2.1)

where fk are consecutive harmonic frequencies,Ak is the amplitude of k-th harmonic, φk is the
phase shift angle and t denotes time.

2.1. Uniform power spectrum

When there is no knowledge about the aircraft dynamics, it should be assumed that the
same amount of information can be obtained from each power spectrum component. Therefore,
the amplitude of each harmonic is designed for the energy content that is uniformly distributed
along the frequency range of interest

Aj,k =Aj
√
pj,k (2.2)

where pj,k is the power of the k-th harmonic component assigned to the j-th flight control.
To increase the efficiency of the input signal, the phase shift angles φk have to be adjusted

byminimizing the Relative Peak Factor

RPF(uj)=
max(uj)−min(uj)
2
√
2rms(uj)

(2.3)

Initial values of the phase angles are obtained by using the Schroeder formula (Schroeder, 1970)

φj,k =φj,k−1+2π(fj,k−1−fj,k)
k−1
Mj
T (2.4)

The Nelder-Mead simplex algorithm (Lagarias et al., 1998) has been used tominimize the non-
-linear cost function.
In order to obtain non-corellated inputs, even harmonics have been assigned to the elevator

and the rest of the set to the ailerons. The first harmonic is omitted due to signal optimization
which is performed later, therefore the lowest frequency is 0.1Hz. The total set is sampled with
a resolution of 0.05Hz and the upper bound of the frequency range of interest is set to 2Hz.
Moreover, to start and finish the control vector with zero amplitude, each input signal is shifted
along the time axis. As a result, we obtain a control vector u shown in Fig. 1.



Power spectrum optimization in the design of multisine manoeuvre... 1195

Fig. 1. Multisine input signals

2.2. Non-uniform power spectrum

In the second case, themultisine inputs are designed by including the a priori knowledge of
the system dynamics. The initial information is given in terms of dimensional stability and con-
trol derivatives. Improvement of the signals efficiency is achieved byusing the a priori knowledge
in order to design apower spectrumwith the energy that ismaximized at significant frequencies.
To select the frequencies at which more power should be stored, the magnitude Bode plots for
all aerodynamic force andmoment components are created. One of the Bode plots used for the
reference power spectrum design is presented in Fig. 2.

Fig. 2. Rolling moment Bode plot for aileron deflection

The reference power spectrumdesigned on theBode plots base andMarchandmethod (Mar-
chand, 1974) is showed in Fig. 3.
In order to optimize the design, the power stored at different harmonics is adjusted by

introducing weighting factors so that the energy distribution minimizing the RPF maximizes
the input efficiency.
A genetic algorithm with linear ranking, uniform crossing and 5% mutation probability has

been used to obtain the final values of the weighting factors.We chose the genetic algorithm for
optimization because the cost function has multiple local minima and the search space is large.
Due tomutual orthogonality of the multisine inputs, this procedure is performed separately for
each control surface. The obtained optimized power spectrum shown in Fig. 4 is used later to
design the control vector u that is presented in Fig. 1.



1196 P. Lichota et al.

Fig. 3. Reference Power Spectrum

Fig. 4. Optimized Power Spectrum

3. System identification

The general aim of the parametric Sys-ID is to determine an adequate mathematical model
which contains unknownvalues of certain system coefficients that have to be obtained indirectly
frommeasured data. To achieve this goal, specific experiments have to be done in order to excite
the system sufficiently and register inputs and object responses. From the modelling aspect, it
means that determined equationsmust provide themodel responseywhichmatches adequately
the measured system response z.
In this paper, to perform the Sys-ID experiment, simultaneous multsine inputs for both

investigated cases are used as excitations for a transport aircraft model with cross-coupling de-
rivatives. To perform parameter estimation, the input signals (elevator and ailerons deflections)
and aircraft response (longitudinal velocity, angle of attack, sideslip angle, roll, pitch, yaw ra-
tes, roll and pitch angles and linear acceleration components) are recorded. As the signals are
corrupted with measurement noise, a 15-th point low-pass digital filter is used (Kendall et al.,
1983)

yn =
1

320
(−3un−7−6un−6−5un−5+3un−4+21un−3+46un−2+67un−1+74un

+67un+1+46un+2+21un+3+3un+4−5un+5−6un+6−3un+7)
(3.1)

3.1. Equation Error Method

Multiple different approaches canbeused toperformtheaircraft Sys-ID, e.g., EquationError
Methods, Output ErrorMethods, Filter ErrorMethods or Artificial Neural Networks. Equation



Power spectrum optimization in the design of multisine manoeuvre... 1197

ErrorMethodsminimize the cost function defineddirectly by the input-output equation.Among
this class, the Ordinary Least Squares (OLS) is the most popular due to its mathematical
simplicity

Y=XΘ+ǫ (3.2)

where Y are dependent variables, Θ denotes unknown parameters, ǫ represent equation errors
and X are independent variables. To obtain the estimates, the cost function which is the sum
of the squares of the residuals is minimized

J(Θ) =
1

2

N∑

k=1

ǫ
2(k) =

1

2
[Y−XΘ]T[Y−XΘ] (3.3)

where k are discrete time points indices.
As a result, the estimates of the unknown parameters are given by

Θ̂=(XTX)−1XTY (3.4)

The block schematic of the OLS is presented in Fig. 5.

Fig. 5. Block schematic of the Least SquaresMethod

Before using theOLS to estimate stability and control derivatives, a data preprocessing step
is performed. It is required to differentiate selected measured object responses and eliminate
biases from the signals. To increase the accuracy of the estimates by eliminating parameters
that have no physical meaning, Eq. (3.2) is solved independently for each output signal. In
the parameter estimation process, the backward elimination is used to determine the model
structure.

3.2. Output Error Method

After performing the Sys-ID using the OLS it turned out that the obtained results are not
satisfying. Regardless of the excitations set, the absolute relative error is considered too high
for Nr. When the identification experiment is designed without a priori knowledge Mu, Mq
and Nβ are estimated with large uncertainty. When the a priori information is included, it is
observed for Np.
As the estimates are biased due to measurement noise in the registered signals, to increase

accuracy of the estimates, the Output ErrorMethod (OEM) is used. This approach is themost
widely applied Sys-ID technique. The main reason for its popularity is its good representation
of a natural formulation for a dynamic system. The block schematic of the OEM is shown in
Fig. 6.
The OEM used in this study is based on the Maximum Likelihood Estimation (MLE),

therefore it seeks for maximizing the probability of observing themeasured responses for model
parameters (Jategaonkar, 2015)

Θ̂=argmax
Θ
p(z|Θ) (3.5)



1198 P. Lichota et al.

Fig. 6. Block schematic of the Output ErrorMethod

where z denotes themeasured response,Θ aremodel parameters and the hat symbol stands for
the estimates. To obtain the estimates, the cost function which is negative log-likelihood, has
been minimized

L(Θ|z) = 1
2

N∑

k=1

[z(tk)−y(tk)]TR−1[z(tk)−y(tk)]+
nN

2
ln(2π)+

N

2
ln(det(R)) (3.6)

where tk is the discrete time at the k-th point, k = 1, . . . ,N, R is the measurement noise
covariance matrix and n denotes the number of model outputs. The unknown measurement
covariance matrixR is estimated from

R̂=
1

N

N∑

k=1

[z(tk)−y(tk)][z(tk)−y(tk)]T (3.7)

Substituting themeasurement covariancematrix intoEq. (3.6) andneglecting the constant terms
allows one to simplify the cost function to the form given by

J(Θ) =
m∏

l=1

1

N

N∑

k=1

(
zi(tk)−yi(tk)

)2
(3.8)

where l=1, . . . ,m is the number of system outputs.
Similarly to the OLS, when theMLE is applied, backward elimination is used to determine

the model structure.

4. Results

To determine if including the a priori knowledge in the multisine inputs design improves the
identification experiment, several criteria are used. For that purpose, deviations from the trim
condition, magnitude coherence function and estimates accuracy are investigated.

4.1. Manoeuvre analysis

To increase themodel efficiencydesignedmanoeuvres shouldallow theaircraft to remainclose
to the equilibriumpoint inwhich it is excited, and the time spent to performaflight experiment
should beminimized. Small deviations from the trim conditions decrease the stabilization time,
so the next manoeuvre can be performed faster.
Deviations of various flight parameters fromthe equilibriumpoint for both investigated cases

are presented in Fig. 7. From the cross-plot, it can be seen that the dispersion of the motion



Power spectrum optimization in the design of multisine manoeuvre... 1199

parameters from their initial values are slightly smaller for the uniformpower spectrum. It leads
to a conclusion that designing the inputs by only optimizing their power spectrum does not
result in shortening the manoeuvre time. Nevertheless, the obtained deviations from the trim
point are small for both cases. This means that both designs are short in time in comparison
to the conventional approach, and even the global aerodynamic parameter estimation can be
performed.

Fig. 7. Cross-plots of flight parameters

To asses the amount of information that can be obtained from a specific manoeuvre, the
magnitude squared coherence function (known as coherence) can be used (Young and Patton,
1988) as it is a measure of the linear dependence between the input and output signals

Γ
2
xy(fk)=

|Gxy(fk)|2
Gxx(fk)Gyy(fk)

(4.1)

whereGxy,Gxx,Gyy are cross, input and output spectral densities.
When there is an ideal linear dependencybetween the input andoutput signal, the coherence

equals 1, whilst 0 means that there is no dependence between those signals. Typically, a value
above 0.6 provides a good object response mapping for multi-output systems (Tischler and
Remple, 2012).
Coherence for the roll rate due to ailerons and pitch rate due to elevator deflection for both

investigated designs is shown in Fig. 8. In our study, we observed that the coherence obtained
for bothmanoeuvres is comparable. As it can be seen from the plots, the coherence is almost 1
in the whole frequency range when the elevator is deflected, and over 0.6 when the ailerons
are used. For the roll rate due to ailerons deflection, the design with the a priori knowledge is
slightly better in the low frequency range, however both manoeuvres provide sufficient amount
of information as the coherence is above the assumed threshold.

4.2. Estimates accuracy

After the whole data set containing excitations (aileron δA and elevator deflections δE) and
model response (longitudinal velocity u, angle of attack α, sideslip angle β, angular rates p,
q, r, attitude angles φ, θ and linear acceleration components ax, ay, az) has been collected and
filtered, Sys-ID is performed. The OLS are used for that purpose. To obtain a mathematical
model of the object, the recorded signals are preprocessed and the backward elimination is
used to determine the model structure. To compare the results obtained for both manoeuvres
(uniform and non-uniform power spectrum), an absolute relative error for each aerodynamic
coefficient is evaluated. The obtained outcomes are presented in Table 1.



1200 P. Lichota et al.

Fig. 8. Coherence function

Table 1.Absolute relative error for Equation ErrorMethod, [%]

Name Uniform Optimized Name Uniform Optimized

Xα 2.74 1.35 Lβ 2.03 2.89

Zα 5.07 5.42 Lp 7.33 4.55

Mu 16.29 8.83 Lr 3.64 1.42

Mα 6.13 3.68 LδA 1.14 0.83

Mq 13.31 8.58 Nα 0.06 1.23

MδE 2.26 1.72 Nβ 10.51 6.23

Yβ 0.15 0.36 Np 6.29 14.47

Yp 0.23 0.03 Nr 31.35 34.18

Yr 1.09 0.02 NδA 2.88 4.13

YδA 0.03 0.05

From Table 1 it can be observed that the absolute relative error values are high for both
investigated cases. It could be caused by the fact that there is a measurement noise in the
registered signals. As it was mentioned earlier, the Sys-ID is performed again using a more
sophisticated method, which is the MLE, because it allows one to include the measurement
noise. The outcomes are presented in Table. 2.

Table 2 shows that the parameters are estimated with high accuracy in both cases. It can
be seen that the non-uniform power spectrum allows one to evaluate slightly better the out-
comes. On the other hand, one has to have in mind that for a non-uniform case, additional
computational time has to be spent to performpower spectrumoptimization. The results of the
Sys-ID performed using theMLE in form of lateral and longitudinal acceleration time histories
are presented in Fig. 9.

In the plots, the recorded object response obtained for both sets of the input signals is
represented by blue lines, whilst the red lines with black markers show the results of the Sys-ID
process. Almost an ideal match between the measurements for all flight parameters is observed



Power spectrum optimization in the design of multisine manoeuvre... 1201

Table 2.Absolute relative error for Output ErrorMethod, [%]

Name Uniform Optimized Name Uniform Optimized

Xα 0.3346 0.2409 Lα 1.8495 1.8454

Zα 0.3134 0.2233 Lβ 0.1540 0.0573

ZδE 2.1896 0.8917 Lp 0.1689 0.1440

Mu 0.5759 0.1646 Lr 0.3220 0.0698

Mα 0.0598 0.1160 LδA 0.1409 0.1221

Mq 0.3982 0.0529 Nα 0.0205 0.0304

MδE 0.2597 0.0846 Nβ 0.0369 0.0471

Yβ 0.0568 0.0228 Np 0.2287 0.0540

Yp 0.2089 0.3757 Nr 0.0329 0.0748

Yr 0.1619 0.0261 NδR 0.2930 0.0426

YδR 0.3917 0.3380

Fig. 9. Time histories forMLE

for both cases. The flight parameters time history plots confirm also our previous conclusion
regarding the stabilization time, i.e.whentheuniformpower spectrumhasbeenused, theaircraft
returned to the trim point slightly faster. Moreover, if the non-uniform energy distribution has
been selected, a new equilibrium point was reached.

5. Conclusions

In this paper, it is shown how to designmultisine input signals with a priori information of the
system dynamics in order to improve the Sys-ID accuracy. Two types of multisine excitations
are used in our investigation.

In the first case, the input signals are obtainedwithout any information of the aircraft dyna-
mics. However, for the majority of Sys-ID tests, there is an a priori knowledge of aerodynamic
coefficients provided that can be used to improve the excitations effectiveness. In our study, this



1202 P. Lichota et al.

knowledge is given in terms of dimensional stability and control derivatives. Those parameters
are used to obtain magnitude Bode plots that allow one to design an optimized input ener-
gy distribution. In order to achieve this aim, a genetic algorithm has bee used. The obtained
non-uniform power spectrum has been used later to design multisine signals.
The model is excited with both sets of inputs, and the Sys-ID is performed for both cases

separately. For estimating the stability and control derivatives, twomethods are used: Equation
Error andOutputError. Themodel parameters obtained from theOLS are biased as the appro-
ach is not resistant to the measurement noise in the independent variables. The secondmethod
allows one to obtain stability and control derivatives with a very high accuracy. On the basis of
those outcomes, it is found that both manoeuvres allow one to obtain accurate estimates. Al-
most an idealmatch between the results andmeasurements is observed. Introducing the a priori
information in the manoeuvre design enables evaluation of slightly better estimates.
However,whenthe initial knowledgeof the systemdynamics is not included, a slightly shorter

time required for stabilization is observed. Nevertheless, in both cases the time that is required
for performing the experiment is greatly reduced in comparison to the conventional approach, so
both designs lower the overall cost of the Sys-ID process. The only noticeable drawback of the
novel approach is that an additional time has to be spent on performing the energy distribution
optimization.

Acknowledgement

This work has been supported by the Universidad de San Buenaventura in Colombia, Bogota and
Deutsches Zentrum für Luft- und Raumfahrt e.V., Flugsystemtechnik Institut, Braunschweig, Germany.

References

1. Albisser M., Dobre S., Berner C., Thomassin M., Garnier H., 2017, Aerodynamic co-
efficient identification of a space vehicle from multiple free-flight tests, Journal of Spacecraft and
Rockets, in press, DOI: 10.2514/1.A33587

2. Fedorov V.V., 1972, Theory of Optimal Experiments, Probability and Mathematical Statistics,
Acade, NewYork, USA

3. GoodwinG.C., PayneR.L., 1977,Dynamic System Identification: Experiment Design andData
Analysis, Academic Press Inc., NewYork, USA

4. Grauer J.A., 2016, Aircraft fault detection using real-time frequency response estimation,
AIAA Guidance, Navigation, and Control Conference, San Diego, USA, AIAA 20160372, DOI:
10.2514/6.2016-0372

5. Jategaonkar R.V., 2015, Flight Vehicle System Identification: A Time Domain Methodology,
2nd ed., Progess in Astronautics and Aeronautics, AIAA, Reston, USA, DOI: 10.2514/4.102790

6. Kalaba R.E., Springarn K., 1982, Control, Identification and Input Optimization, Mathema-
tical Concepts and Methods in Science and Engineering, Plenum Press, New York, USA, DOI:
10.1007/978-1-4684-7662-0

7. KendallM.G., StuartA., Ord J.K., 1983,The Advanced Theory of Statistics, 4th ed., Vol. 3,
Charles Griffin (Edit.), London, UK, DOI: 10.1002/for.3980040310

8. Lagarias J.C., Reeds J.A., Wright M.H., Wright P.E., 1998, Convergence properties of
the Nelder-Mead simplexmethod in low dimensions, SIAM Journal of Optimization, 9, 1, 112147,
DOI: 10.1137/S1052623496303470

9. Levin M.J., 1960, Optimum estimation of impulse response in the presence of noise, IRE Trans-
actions on Circuit Theory, 7, 1, 5056, DOI: 10.1109/TCT.1960.1086622

10. Lichota P., 2016, Inclusion of the D-optimality in multisine manoeuvre design for aircraft para-
meter estimation, Journal of Theoretical and AppliedMechanics, 54, 1, 8798,DOI: 10.15632/jtam-
-pl.54.1.87



Power spectrum optimization in the design of multisine manoeuvre... 1203

11. Lichota P., Sibilski K., Ohme P., 2017, D-optimal simultaneous multistep excitations for air-
craft parameter estimation, Journal of Aircraft, in press, DOI: 10.2514/1.C033794

12. Marchand M., 1974, Untersuchung der Bestimmbarkeit von Fleugzeugderivativen aus Flugver-
suchen, DLR, Report DFVLR-IB 154-74/32, Braunschweig, Germany

13. Morelli E.A., 2009, Flight test experiment design for characterizing stability and control of hy-
personicvehicles,Journal ofGuidance,Control andDynamics,32, 3, 949959,DOI: 10.2514/1.37092

14. Morelli E.A., Klein V., 1990, Optimal input design for aircraft parameter estimation using
dynamicprogrammingprinciples,AIAAAtmospheric FlightMechanicsConference, Portland,USA,
DOI: 10.2514/6.1990-2801

15. Nahi N.E., Wallis D.E.J., 1969, Optimal inputs for parameter estimation in dynamic systems
with white observation noise,Proceedings of Joint Automatic Control Conference, 506512

16. SchroederM., 1970, Synthesis of low-peak-factor signals and binary sequenceswith lowautocor-
relation, IEEE Transactions on Information Theory, 16, 1, 8589, DOI: 10.1109/tit.1970.1054411

17. Seren C., Bommier F., Verdier L., Bucharles A., Alazard D., 2006, Optimal experi-
ment and input design for flight test protocol optimization,AIAAAthmospheric Flight Mechanics
Conference, Keystone, USA, DOI: 10.2514/6.2006-6280

18. Seren C., Hardier G., Roos C., 2013, Swarm intelligence and system identification: a hybrid
discrete jumping frogs algorithm for optimal input design, 3rd IFAC International Conference on
Intelligent Control, Chengdu, China

19. Tischler M.B., Remple R.K., 2012, Aircraft and Rotorcraft System Identification, 2nd ed.,
AIAA Education Series, AIAA,Washington D.C., USA, DOI: 10.2514/4.868207

20. Viana M.V.P., 2016, Time-domain system identification of rigid-body multipoint loads model,
AIAA Atmospheric Flight Mechanics Conference, Washington D.C., USA, AIAA 20163706, DOI:
10.2514/6.2016-3706

21. Wells W.R., Ramachandran S., 1977,Multiple control input design for identification of light
aircraft, IEEE Transactions on Automatic Control, 22, 985987, DOI: 10.1109/TAC.1977.1101653

22. Young P., Patton R.J., 1988, Frequency domain identification of remotely-piloted helicopter
dynamics using frequency-sweep and Schroeder-phased test signals,AIAAAtmospheric Flight Me-
chanics Conference, Minneapolis, USA, AIAA 884349CP,DOI: 10.2514/6.1988-4349

Manuscript received February 27, 2017; accepted for print May 4, 2017