Nonparametric System Identification of an Autonomous Helicopter

of an Autonomous Helicopter

FREDRIK SVENSON

The unmanned helicopter Skeldar relies on model based control to per-form its tasks. System identification is an integral component in the process of deriving models of helicopter dynamics. This thesis aims to investigate how nonparametric methods of system identification can support the current modelling and system identification practices at Saab.

Nonparametric system identification does not require a pre-defined model structure. Models estimated with this methodology may be used to validate parametric models, which are necessary for the implemen-tation of the model based control system. This thesis examines several nonparametric methods of system identification in both the frequency and time domains. The theory of these methods is presented and their performance is analyzed on data from flight tests as well as from simu-lated systems.

Analysis of the results shows that models are highly dependent on the choice of input signal spectrum. To best take advantage of nonpara-metric system identification in this application, experiments should be designed with special regard to the system properties sought to be mod-elled. Nonparametric system identification can then be used to provide a good understanding of the system properties in the excited frequency region.

Den obemannade helikoptern Skeldar litar till modellbaserade styrla-gar för att utföra sina uppgifter. Systemidentifiering är en viktig del i processen för att härleda modeller av helikopterdynamik. Detta exa-mensarbete syftar till att undersöka hur icke-parametriska metoder för systemidentifiering kan stödja de nuvarande metoder för modellering och systemidentifiering på Saab.

Icke-parametrisk systemidentifiering kräver inte en fördefinierad mo-dellstruktur. Modeller skattade med denna metod kan användas till att validera parametriska modeller, som är nödvändiga för implementering-en av de modellbaserade styrlagarna. Detta examimplementering-ensarbete undersöker flera icke-parametriska metoder för systemidentifiering i både frekvens-och tidsdomänen. Teorin för dessa metoder presenteras frekvens-och dess pre-standa analyseras på data från flygprov samt simulerade system.

Analys av resultaten visar att modellerna är mycket beroende av va-let utav insignalspektrum. För att bäst utnyttja icke-parametrisk syste-midentifiering i denna tillämpning bör experimenten utformas med sär-skild hänsyn till de systemegenskaper man avser modellera. Då kan icke-parametrisk systemidentifiering sedan användas för att ge en god förståelse av systemegenskaper i det exciterade frekvensområdet.

I wish to send my gratitude to Sam Nicander, for accepting me as an intern at Saab and introducing me to the Skeldar project. My special thanks go out to my supervisors at Saab, Dr. Ola Härkegård and Dr. Patrick Borgqvist, for always giv-ing thoughtful feedback and shargiv-ing their extensive knowledge of automatic control matters and its application to helicopters. I would also like to thank Fredrik Eskils-son, Erik Backlund and the rest of the software team at Skeldar for always being ready to answer questions and providing an interesting and fun work environment. I would also like to express my gratitude towards my supervisor and my ex-aminer at KTH, Per Hägg and Prof. Håkan Hjalmarsson. Though most of my time was spent in Linköping, the weeks I did spend at KTH allowed for good dis-cussions and valuable feedback. Valuable disdis-cussions were also had at LiTH with Prof. Lennart Ljung, whom I would like to thank for the special insight into system identification he brought.

Finally, my thanks go out to my family for their help with proof-reading and for always supporting me in the choices I make in life.

Fredrik Svenson

Nomenclature

1 Introduction 1

1.1 Background . . . 1

1.2 Problem Description . . . 2

1.3 Report Outline . . . 3

2 Nonparametric system identification 4 2.1 Preliminaries . . . 4

2.1.1 Time-domain methods . . . 4

2.1.2 Frequency-domain methods . . . 5

2.2 Spectral analysis and frequency windowing . . . 7

2.3 Local Polynomial Method . . . 8

2.4 Transient Impulse Response Modelling Method . . . 9

2.5 Regularized Finite Impulse Response . . . 10

3 Helicopter dynamics 12 3.1 Rigid-body dynamics . . . 12

3.2 Main rotor dynamics . . . 14

3.3 Tail rotor yaw dynamics . . . 17

4 Results 18 4.1 Experimental setup . . . 18 4.1.1 Pre-processing of data . . . 19 4.2 Yaw dynamics . . . 20 4.2.1 Frequency-domain methods . . . 20 4.2.2 Time-domain methods . . . 21

4.3 Simulated first order time-delayed system . . . 24

5.1.1 Frequency-domain approach . . . 37

5.1.2 Time-domain approach . . . 38

5.2 Conclusion . . . 38

5.3 Future work . . . 39

Ω Rotor blade angular speed Ψ Rotor blade angular position

β Rotor blade flapping angle δlat Cyclic roll control input δlon Cyclic pitch control input δped Tail rotor control input

φ Roll angle ψ Yaw angle θ Pitch angle p Roll rate q Pitch rate r Yaw rate

FIR Finite Impulse Response FRF Frequency Response Function LPM Local Polynomial Method

Introduction

1.1

Background

Recent years have seen a great increase in the application of unmanned aerial vehi-cles (UAVs) and the development of such vehivehi-cles. Most applications for UAVs are still seen within the defense industry. However, as flying and airspace regulations are re-evaluated it is likely that an entrance of UAVs into the civilian market will be seen in the coming years.

While fixed-wing UAVs in many shapes and forms are already common, rotor-craft UAVs such as helicopters have seen comparatively little use. Helicopters, with their ability to take off and land vertically, as well as to hover, already have a unique role amongst aircraft. Hence, translating these qualities of the helicopter into a smaller-scale autonomous platform can be considered an attractive venture.

Saab AB, the company at which this thesis has been conducted, has identified the need for a helicopter UAV and has developed the Skeldar platform. The main focus of Saab has previously been on fixed-wing aircraft whose flight characteristics are radically different to that of helicopters. An initial study and development of a helicopter simulation model was done in [1].

Skeldar is capable of carrying out missions fully autonomously, thus putting high demands on the control systems. Model based control systems have reached the point where they can be considered to be standard in industrial applications. This is true also for Skeldar, where a model based control system is implemented. The use of a model based control system puts strict requirements on the models used for control design to be as accurate as possible. The feasibility of using an adaptive method of control has been studied in a previous master thesis [7]. Well performing control systems also put great reliance on high quality sensor data for use in the feedback loop. The sensors used in Skeldar and methods of data fusion were studied in [8].

to the system. For helicopters, such systems of equations tend to end up being high-order, include many parameters and be nonlinear. For control design applications it is often sufficient to use linear models. The theory for control of linear systems is very well developed compared to that of nonlinear control. First-principles models of helicopter dynamics can often be well approximated by linearizing around a working point, such as hover flight, and be used to great success in control design [12].

System identification makes use of experimental input-output data to produce a mathematical model of the system. This approach is thus more direct, as the actual behaviour of the system is used to derive the model. In the event that extensive experimental data is available, system identification can be a great help to recognize certain characteristics. These characteristics can then be taken into account at an early stage of more comprehensive first-principles modelling.

Neither of the two approaches to system modelling excludes the other. In gen-eral they should be seen as a complement to one another and each be used to their strengths. This thesis will focus on a subset of system identification, namely nonparametric system identification.

1.2

Problem Description

The most common approach to system identification is the use of parametric model structures. In this type of system identification a model structure with a number of free parameters is pre-defined. These parameters are then chosen, often through solving a least squares problem, such that the user is left with the model producing the best fit to experimental data. One of the main difficulties in parametric system identification is to choose the best model structure.

In nonparametric system identification, as implied by the name, the identifi-cation is not restricted to a small amount of parameters. These methods make use of the information contained in the data directly to estimate frequency and impulse responses. An advantage of using nonparametric methodologies is thus to gain knowledge of the true system without dependency on model structure choice. The aim of this thesis is to investigate the system identification process which is used today at Saab to produce models of Skeldar. More specifically, this thesis aims to find how nonparametric system identification methods can be applied to current methodology. Some specific goals of this thesis are summarized below:

• study available nonparametric system identification methods, • evaluate the performance of these methods on available data,

• consider how nonparametric system identification can be used to validate the current system identification methodology used at Saab.

1.3

Report Outline

The second chapter of this thesis describes the investigated methods. A more in depth analysis of the differences between parametric and nonparametric methods is conducted. Methods to be analyzed include the frequency response function meth-ods Local Polynomial Method (LPM) and Transient Impulse Response Modelling Method (TRIMM). A method using spectral analysis and frequency windowing is also investigated. Finally, a method of identifying a regularized finite impulse re-sponse is investigated.

The third chapter describes the relevant helicopter dynamics to be identified. There are two specific subsystems from which data for system identification has been gathered, the cyclic dynamics of the main rotor as well as the dynamics of the tail rotor.

In the fourth chapter the previously presented system identification methods are applied to flight test data. The performance of each of the methods are analyzed and compared.

Nonparametric system identification

This chapter describes the principles of nonparametric system identification. Some differences and advantages compared to parametric methods are described and a number of specific methods used in nonparametric system identification are pre-sented.

2.1

Preliminaries

In the literature it is established that a linear time-invariant (LTI) system can be described by its transfer function or corresponding impulse response [10]. Nonpara-metric system identification is a technique where transfer functions and impulse responses can be estimated without selecting a model structure. The transfer func-tion from input to output of a given system is denoted as G0(q) where q is the

time-shift operator. This is the function to be identified.

2.1.1 Time-domain methods

Consider a system described by

y(t) = G0(q)u(t) + v(t) (2.1)

where v(t) is additive noise. The transfer function for this system is given by

G0(q) = ∞

X

k=1

g_k0q−k. (2.2)

Here, the coefficients g0

k, k= 1, . . . , ∞, describe the impulse response of the system.

the model. The model to be estimated is then given as G(q) = n X k=1 g0_kq−k (2.3)

where n is the order of the FIR model. The impulse response coefficients, which are to be identified, can be gathered in the parameter vector θ = [g1, g2, . . . , gn]T.

This vector θ can now be estimated using a Least Squares (LS) approach. Least squares is a common methodology used in system identification and its application to FIR models is described in [2]. With given N samples of input-output data

y(t), u(t), t = 1, . . . , N the FIR model can be stated as

y(t) = φ(t)Tθ+ v(t), φ(t) = [u(t − 1) . . . u(t − n)]T. (2.4)

This can be rewritten on vector form as YN = ΦTNθ+ ΛN, where YN = [y(n + 1) y(n + 2) . . . y(N)]T

ΦN = [φ(n + 1) φ(n + 2) . . . φ(N)]

ΛN = [v(n + 1) v(n + 2) . . . v(N)]T.

The LS solution is now given by: ˆθLS

N = [ˆg1LSˆg2LS . . . ˆgnLS]T = arg min νN(θ) (2.5a) νN(θ) = kYN −ΦTNθk2= N X t=n+1 (y(t) − φT_(t)θ)2 _(2.5b) ˆθLS N = (ΦNΦTN)−1ΦNYN = R−1N FN (2.5c) FN = ΦNYN = N X t=n+1 φ(t)y(t) (2.5d) RN = ΦNΦTN = N X t=n+1 φ(t)φ(t)T. (2.5e)

Note that summation in equation (2.5b) starts at n + 1. This is due to the inputs

u(−n + 1), . . . , u(0) not being known.

The FIR model structure can in fact be seen as a special case of either of two common model structures, the Output Error (OE) and the Autoregressive Exoge-nous (ARX) structures [10]. Similarly to those parametric model structures, the FIR model can be evaluated in the time-domain on validation data.

2.1.2 Frequency-domain methods

It is well known that the transfer function has useful properties in the complex plane. The complex number G(eiω_{) carries information about the amplitude and}

these properties are especially interesting for system identification purposes as they provide information about the system’s behaviour in certain frequency ranges. As shown in [10], an estimate of G(eiω_{) can be gained by applying the discrete Fourier}

transform (DFT) to time-domain data.

When only a finite record of data is available, a system as described by (2.1) can be stated as

y(t) = G0(q)u(t) + v(t) + transient (2.6)

where transient terms are also taken into account. Transient terms exist as a result of the system leaving a state of equilibrium due to the excitation by the input. The frequency-domain representation YN(ω) of a signal y(t) at frequency ω is given by

applying the DFT as follows:

YN(ω) = √1 N N X t=1 y(t)e−iωt. (2.7)

By applying the DFT as in (2.7) to (2.6), the following frequency domain represen-tation of the system excited by a sinusoid is acquired:

YN(ω) = G0(eiω)UN(ω) + VN + T (2.8)

Here VN and T correspond to the transformed noise and transient terms,

respec-tively. Again, the transient term is present due to a limited data set. A pure periodic signal causes transient terms to tend to zero.

A common first approach to estimating the transfer function is given by ˆ

GN(eiω) =

YN(ω) UN(ω)

. (2.9)

Here UN(ω) 6= 0 is assumed. For frequencies where this is untrue, the estimate

is regarded as undefined. The estimate ˆGN(eiω) is referred to as the empirical transfer-function estimate (ETFE). It is called empirical since no assumptions on

the system other than linearity has been imposed [10, pp. 171-173].

Due to the fact that no assumptions apart from linearity are made on the ETFE, this estimate can be very crude. Unless the system is excited by a periodic input, the variance of the estimate does not decrease with N. Moreover, the transient terms of (2.6) are not taken into account. Due to the nature of identification data being limited, such terms will have an impact on the identification results.

2.2

Spectral analysis and frequency windowing

To alleviate the problem of great variance typically present in the ETFE, a bias is introduced which is based on the assumption that the true transfer function G0(eiω)

is a smooth function of ω. If G0(eiω) is assumed to be constant over the interval

2πk1

N = ω0−∆ω < ω0 < ω0+ ∆ω =

2πk2

N (2.10)

and ∆ω is assumed small, then a weighted average of the estimates of G at the two frequencies 2πk1/N and 2πk2/N can be formed to estimate the assumed constant

behaviour between them. The measurements are weighted according to their inverse variance which is given by

αk=

|U_N(2πk/N)|2

Φv(2πk/N) (2.11)

where the numerator is the input periodogram and Φv(2πk/N) is the noise

spec-trum. However, if the noise spectrum can be assumed to be constant across the interval of frequencies affected by the smoothing action, the impact of the noise spectrum cancels out in the estimate.

With this, and by introducing the windowing function Wγ(ξ) which is used to

pay more attention to frequencies near ω0 when smoothing, the smoothed estimate

can be expressed as ˆ GN(eiω0) = Rπ −πWγ(ξ − ω0)|UN(ξ)|2GˆET F EN (eiξ)dξ Rπ −πWγ(ξ − ω0)|UN(ξ)|2dξ . (2.12)

This estimate is connected to spectral analysis by the following reasoning as described in [10, pp. 180-181]. Due to the denominator of (2.12) being a weighted average of the periodogram |UN(ξ)|2, it is found that, as N → ∞,

Z π −π Wγ(ξ − ω0)|UN(ξ)|2dξ → Z π −π Wγ(ξ − ω0)Φu(ξ)dξ (2.13)

where Φu(ω) is the spectrum of u(t). If Φu(ω) does not change much over the

frequency interval defined by Wγ(ξ), the left side of (2.13) may be interpreted as a

estimate of the input spectrum as follows ˆΦN

u(ω0) =

Z π

−πWγ(ξ − ω0)|UN(ξ)|

2_{dξ (2.14)}

By the same reasoning, and since |U_N(ξ)|2_GˆET F E

N (eiξ) = |UN(ξ)|2 YN(ξ)

UN(ξ) = YN(ξ) ¯

UN(ξ), (2.15)

the numerator of (2.12) can be stated as the cross spectrum of input and output ˆΦN

yu(ω0) =

Z π

Therefore, the frequency response function (FRF) estimate in (2.12) is the ratio of two spectral estimates:

ˆ GN(eiω0) = ˆΦN yu(ω0) ˆΦN u(ω0) (2.17) The result of the above is implemented in Matlab as the function spafdr. A user may specify the frequency vector ωk, k = k1, . . . , kn as the set of frequency

points around which the estimate will be locally smoothed. The windowing function

Wγ(ξ) is controlled by specifying a frequency resolution. This resolution may either

be a scalar defining the resolution of the whole frequency interval, or a vector of the same length as ωk where each entry is a local, frequency-dependent resolution. The

resolution controls the trade-off of bias and variance, a finer resolution may capture finer details but also leads to greater variance.

2.3

Local Polynomial Method

Two closely related methods of nonparametric system identification in the frequency domain have recently been developed [4, 6]. In this section the Local Polynomial Method (LPM) will be presented.

Similar to what is presented in Section 2.2, the LPM also takes advantage of the smoothness of the transfer function G0. However, unlike the previously introduced

methods, the LPM takes the transient term T in (2.8) into account. This transient term, present due to initial action of the transfer function G0, can also be shown to

be smooth [4].

The idea of the LPM is that, based on this smoothness, the functions G0 and T

can be approximated by a complex polynomial in a narrow frequency band centered around a frequency Ωk= ei2πk/N. This leads to N local least squares (LS) problems

with a set of parameters to be solved for in each frequency window. The polynomial parameters are estimated from experimental data in that frequency band and Gk is

retrieved as the FRF estimate at frequency Ωkfrom this local polynomial model [4].

The polynomial representation of Gk and Tk for frequency lines k + r, with r= −L, . . . , 0, . . . , L is given as: Gk+r = Gk+ R X s=1 gs(k)rs (2.18) Tk+r = Tk+ R X s=1 ts(k)rs (2.19)

Here the polynomials are limited to order R. The terms to be identified are then the frequency response coefficients Gk and Tk as well as the complex polynomial

coefficients gs and ts.

polynomials R. Similarly to the spafdr these parameters control the bias-variance trade-off. Here, bias increases with the width 2L + 1 of the window and decreases with the degree R of the polynomials. The limitation on how these parameters can be chosen is that a full rank set of equations for the LS problem must be available. Thus the parameters must be chosen according to the constraint L ≥ R + 1 [5]. In the estimates done in Chapter 4 of this thesis, these parameters are chosen as L = 5 and R = 2. These parameters are chosen due to displaying consistent results on the test data.

2.4

Transient Impulse Response Modelling Method

Similar to the Local Polynomial Method presented in Section 2.3, the Transient Impulse Response Modelling Method (TRIMM) attempts to achieve a more accu-rate FRF estimate by taking into account the frequency leakage errors caused by transient terms.

TRIMM follows a similar approach to the LPM. The frequency response Gkand

transient response Tk are estimated by interpolating over a band of neighbouring

frequencies. The main difference is that TRIMM, instead of assuming a local poly-nomial behaviour, uses the impulse response characteristics of Gk+r and Tk+r as

follows:

Gk+r= Gk+

∞

X

t=1

g(t)[e−iωk+rt_{− e}−iωkt] (2.20)

Tk+r= √1 N N −1 X t=0 τ(t)e−iωk+rt (2.21)

Here, the frequency response Gk and impulse response coefficients g(t) and τ(t) are

estimated. The width of the frequency window is defined by r = −L, . . . , 0, . . . , L. If this frequency window is narrow, then the behaviour around Ωk may be well

approximated by low order FIR models [5]. Thus, the equations (2.20) and (2.21) can be truncated to

Gk+r= Gk+ n1 X

t=1

g(t)[e−iωk+rt_{− e}−iωkt] (2.22)

Tk+r= √1 N n2−1 X t=0 τ(t)e−iωk+rt (2.23)

where n1 and n2 are design parameters deciding how many impulse response

coef-ficients are to be estimated. Again, the design parameters n1, n2 and L govern the

trade-off between bias and variance. A larger frequency window 2L + 1 decreases variance but increases bias. Increasing the number of impulse response coefficients

n1 and n2 increases variance error. Thus, if a large amount of noise variance is

A second characteristic differentiating TRIMM from the LPM is that instead of solving N small LS problems locally, TRIMM assembles all equations for the DFT frequencies ω1, . . . , ωN and solves a large LS problem for a vector of 2N + n1 +

n2 unknowns. This means that while the LPM has different parameters for each

individual frequency, TRIMM considers the global problem and thus only involves a one set of unknown parameters for the full transfer function. One implication of this is that TRIMM may suffer from larger computation times than the LPM, but in presence of high noise levels, it may provide an additional means of smoothing by solving across the full length of the experiment. In the estimates done in Chapter 4 of this thesis, the tuning parameters are chosen as L = 3 and n1 = n2 = 15.

2.5

Regularized Finite Impulse Response

Recent development has shown that the time-domain nonparametric approach of estimating a FIR as presented in Section 2.1.1 may be refined to show good iden-tification results [2, 11]. This refinement consists of adding regularization to the FIR model estimate. The reasoning for this stems from the fact that the variance of a FIR model as in (2.3) increases linearly with model order n [2]. The regular-ization acts to counter the variance by introducing a bias, based on the reasonable assumption that impulse response coefficients bear some relation to eachother.

Regularization is added by introducing the positive semi-definite matrix D in the criterion νN(θ) in (2.5) to get

ν_NR(θ, D) = N

X

t=n+1

(y(t) − φT_(t)θ)2_{+ θ}T_{Dθ. (2.24)}

This changes the estimate and relates it to the unregularized estimate ˆθLS N as

ˆθR

N = [ˆg1Rˆg2R. . .ˆgnR]T = (RN + D)−1FN = (RN + D)−1RNˆθLSN . (2.25)

The regularization matrix D can be chosen in several ways. This is essentially a tuning parameter and the current Matlab implementation has seven different so called “regularization kernels” implemented. The intuition behind the choice of kernel is that some correlation exists amongst the regressors in matrix RN. For

instance, if the system to be identified is stable, the parameters on the diagonal of D should increase exponentially. From that assumption a kernel can be formed from a set of hyperparameters, which may be determined through solving an optimization problem. Previous results have shown that for systems such as those investigated in this thesis, i.e. slow and has a high signal to noise ratio, choice of kernel has a rather small impact on the general performance of the method [2].

added as

G(q, η, θ) = Gb(q, η) + Gr(q, θ) (2.26)

where Gb(q, η) is the base-line model and Gr(q, θ) is a FIR of order n as defined

previously. The base-line model is a low-order auto regressive on OE or ARX form, with an associated parameter vector η, which is estimated separately. Using this model a preliminary output

yb(t) = Gb(q, ˆηN)u(t) (2.27)

is first estimated which then allows for the residual output yr(t) to be defined as

yr(t) = y(t) − yb(t). (2.28)

Using this, the FIR model Gr(q, θ) is estimated from the input u(t) and residual

output yr(t) with the regularization method as in (2.24).

Helicopter dynamics

Unlike fixed-wing aircraft such as airplanes, helicopters do not have an inherent “will” to fly. While an airplane relies on its wings to produce lift, a helicopter is very likely to drop out of the sky should its systems fail. In this chapter, the dynamics allowing a helicopter to fly will be described. However, a complete description of all physical interactions and technical properties of the helicopter is outside the scope of this thesis. This chapter will thus be limited to describe the dynamics of movement in a state of hovering or low speed.

3.1

Rigid-body dynamics

Most of the states relevant to helicopter flight in its body-fxed reference frame are illustrated in Figure 3.1. These states, with the origin at the helicopter center of gravity, include:

u, v, w the velocities in the fuselage coordinates, p, q, r the roll, pitch and yaw angular rates,

φ, θ, ψ the roll, pitch, and yaw attitude angles about the principal fuselage axis, βs, βc the lateral and longitudinal rotor flapping angles,

Ψ, Ω the rotor blade angular position and the rotor angular velocity.

A parameterized linear model can be developed to provide an understanding for the physical interactions involved. This model can then be used as a base for system identification. The relationships presented here closely follow the work presented in [12] and [13].

x, u, φ, p y, v, θ, q z, w, ψ, r Ψ, Ω -βs βc hub plane c.g.

Figure 3.1. Illustration of the body-fixed frame of reference.

the Newton-Euler equations of motion. These equations , for constant mass m and moment of inertia I, are expressed in the body-fixed reference frame as follows:

m ˙v+ m(ω × v) = F (3.1)

I ˙ω+ (ω × Iω) = M. (3.2)

where v = [u v w]T _{and ω = [p q r]}T _{are the fuselage velocities and angular rates.}

The vector F = [X Y Z]T _{consists of the external forces acting on the center of}

gravity and M = [L M N]T _{is the vector of external moments.}

From equations (3.1) and (3.2) follow six first-order nonlinear differential equa-tions describing the rigid-body motion. These equaequa-tions can be represented as a nonlinear vector differential equation

˙

x= f(x, u) (3.3)

where x is the vehicle state vector and u is the control input vector. The state vector for the rigid-body dynamics is defined as

x= [u v w φ θ ψ p q r]T. (3.4)

The helicopter control variables are lateral and longitudinal cyclic rotor controls

ulat and ulon, the collective pitch ucol and the tail rotor control input uped. The

control input vector is then

As described earlier, this thesis will consider the helicopter in a hover flight condi-tion. This means that the nonlinear vector differential equation of motion in (3.3) can be linearized about an equilibrium state defined by zero velocity and angular velocity v0= ω0 = [0 0 0]T. The linearized form of the equations of motion can be

expressed in state-space form as

δ ˙x= Aδx + Bδu. (3.6)

where δ denotes a small deviation from the equilibrium state. About the equilibrium point the equations are given as

δ ˙v= ∆F /m (3.7)

δ ˙ω= ∆M/Ixx,yy,zz. (3.8)

where ∆F and ∆M are defined as in (3.1) and (3.2), with ∆ denoting a small deviation from the equilibrium. It is desired to express the external forces and moments as functions of vehicle states and control variables. This can be done with Taylor series expansion, and due to the wish for linearity, only the first order terms are used. For instance, the longitudinal force component can then be expressed as

∆X = ∂X ∂uδu+ ∂X ∂vδv+ ∂X ∂pδp+ · · · + ∂X ∂δlat δlat+ ∂X ∂δlon δlon+ · · · . (3.9)

The partial derivatives of forces or moments with respect to the helicopter states are called stability derivatives, and are abbreviated as Xu, Xv and so on. The

deriva-tives with respect to inputs are called control derivaderiva-tives, and are abbreviated as

Xδlat, Xδlon and so on. These stability and control derivatives can then be collected

in the system and control matrices, A and B, as parameters to be identified. The representation given above has the implication that the rotor control forces and moments are instantaneous and proportional to the control inputs. However, this representation is limited due the fact that the rotor is itself a dynamic system responding to control inputs and the helicopter state. This leads to a dynamic coupling between the main rotor system and the helicopter fuselage dynamics. To take the higher order effects of coupling into account, the model of the rigid-body dynamics can be extended with a simplified model of the rotor dynamics.

3.2

Main rotor dynamics

The model of the rigid-body dynamics is extended by first introducing the concept of a tip-path-plane (TPP). The TPP is the plane defined by the tips of the rotor blade. A main force responsible of generating helicopter movement is the thrust vector T . The main component of this force is perpendicular to the TPP which in turn is subject to control input.

instead produced by varying the blade pitch angle Θ, which can be expressed as a function of the rotor blade angle Ψ as

Θ(Ψ) = Θ0− A1cos Ψ − B1sin Ψ. (3.10)

Here, Θ0 is the collective (average) pitch angle while A1 and B1 is the amount of

pitch the blade undergoes when it is located above the tail and on the right-hand side, respectively. The lateral and longitudinal cyclic controls, δlat and δlon, are

related to A1 and B1 as

A1 = Blatδlat B1= Alonδlon. (3.11)

A swashplate mechanism forms the basis of the blade pitch control. The collective and cyclic pitch control as handled by the swashplate is illustrated in Figure 3.2 and Figure 3.3.

Non-rotating part

Rotating part

Increased pitch at blade Increased pitch at blade

Ω

Rod pushed up Rod

pushed up

Figure 3.2. Illustration of the principles of collective pitch control [1].

Increased pitch at blade

Non-rotating part Rotating part Ω Rod pushed up Rod pulled down Decreased pitch at blade

Figure 3.3. Illustration of the principles of cyclic pitch control [1].

flapping motion of the rotor blade. The flapping angle of the blade is denoted β and this motion causes the TPP to tilt, and with it tilts the rotor thrust vector T . This means that through the cyclic control input, the pilot has indirect control of the direction of rotor thrust. This is also the primary source of manoeuvrability of the helicopter [13]. The dynamics of the rotor blade flapping angle are thus an important aspect to take into account in the model. The TPP representation is illustrated in Figure 3.4.

β

β

β

β

β

β

Figure 3.4. Illustration of the tip-path-plane representation and flapping motion.

Due to the periodicism of the rotor blade pitch angle Θ, the flapping angle β can also be described as a periodic function in Ψ as

β(Ψ) ≈ β0(t) − βc(t) cos Ψ − βs(t) sin Ψ. (3.12)

This is a truncated Fourier series, higher order harmonics can be omitted since the first-harmonic representation provides sufficient accuracy for the purpose of deriving an identification model [3, p. 153]. The non-periodic term β0 represents the coning

angle present during hovering, while βcand βsdescribe the tilting of the TPP in the

longitudinal and lateral directions, respectively. Through some manipulation, which is omitted here, the flapping equations of motion can be described by a first-order system in terms of βcand βs [13, p. 20]:

τf ˙βc= −βc− τfq+ p Ω+ Aβsβs− B1 (3.13) τf ˙βs= −βs− τfp − q Ω− Bβcβc+ A1 (3.14)

Here τf is the rotor time constant and Aβs and Bβc are derivatives describing the

The TPP dynamics described above are coupled to the rigid-body dynamics by the forces and moments produced by the tilted thrust vector. These forces and moments can be expressed as stability derivatives by similar reasoning as in (3.9). The rigid-body fuselage equations and the rotor TPP equations are then coupled by the flapping derivatives. By extending the model with rotor dynamics, the input derivatives Xδlon, Yδlat, Mδlon and Lδlat are replaced by the flapping derivatives Xβc, Yβs, Mβc and Lβs. This means that cyclic input commands go into the rotor

dynamics.

A model suitable for identifying the roll and pitch dynamics can be described by a 4th order system of equations. The analysis and evaluation performed in Chapter 4 is based on this assumption. This simplified system can be stated as:

˙p = Lβsβs (3.15a)

˙q = Mβcβc (3.15b)

τf ˙βc= −βc− τfq+ Aβsβs+ Aδlonδlon+ Aδlatδlat (3.15c) τf ˙βs= −βs− τfp − Bβcβc+ Bδlonδlon+ Bδlatδlat (3.15d)

Here the effects of translational movement are not taken into account and some off-axis flapping response terms have been dropped. However, the cross-coupling derivatives Alatδlat and Blonδlon are added to model off-axis effects. This system

of equations describe the coupled rotor-fuselage dynamics. Coupled systems are known to exhibit second order characteristics in the frequency response which is to be kept in mind when analyzing identification results.

3.3

Tail rotor yaw dynamics

Compared to the dynamics of the main rotor and rigid-body fuselage coupling de-scribed in the previous section, the yaw dynamics is a much simpler system to model. The equation of motion for the yaw rate r is drawn from the Euler equa-tion in (3.8). When disregarding the contribuequa-tion from the main rotor rotaequa-tion, a stability derivatives model is stated as

˙r = Nrr+ Npedδped. (3.16)

The damping moment produced by the tail rotor and tail boom is modeled by the stability derivative Nr. In this case, the tail rotor dynamics do not need to

be modeled explicitly. The small size of the rotor as well as the high rotor speed allows the thrust response to be considered as quasi-instantaneous and is modeled by Nped [12]. However, the unmodeled dynamics present do still have some effect,

Results

In this chapter the nonparametric methods of system identification presented in Chapter 2 are evaluated on data gathered from flight tests with the Skeldar. All Matlab routines referenced in this chapter are run with the R2013a version and are called from the System Identification Toolbox [9].

In Section 4.1 the experimental setup and preprocessing of data is described. In Section 4.2 the results of the yaw dynamics is presented and in Section 4.3 the results of applying the methodology to a simulated system are presented. Finally, in Section 4.4 the results seen on the cyclic dynamics are presented.

4.1

Experimental setup

For best performance in system identification experiments, the experiment is typi-cally run in open loop [10]. However, in many industrial applications this is often not possible. This is true for helicopters and the Skeldar as well, as it is an inher-ently unstable platform. Due to this, experiments are instead run in closed loop with a feedback controller. Performing system identification on data collected in closed loop is itself a topic of research, and in general most difficulties are depen-dent on the presence of disturbances. In Skeldar, the sensors and measurements are refined to the point that the impact of high frequency noise in the loop is minimal. Disturbances that do have an impact are related to wind. These disturbances are of low frequency and are attenuated by the feedback controller.

To further minimize the impact of the feedback in the system, test signals are applied as an input disturbance on the control signal. In Figure 4.1 the experimental setup is illustrated in a block diagram. For the purpose of system identification performed in this thesis, the input signal uc and output signal y are used. During

experiments the reference r is kept constant and the predefined test signals are added as the input disturbance ud. The controller used has sufficiently low bandwidth to

not act on the frequency range defined by the ud signal. This setup allows for

Controller System

r e

ud

u uc y

−

Figure 4.1. Block diagram of the system and involved signals.

4.1.1 Pre-processing of data

Before collected data are used for system identification, some pre-processing is done. Due to the desire to maintain a state of hover, the input signal may be subject to an offset. This offset is removed to focus the mean at 0. Some low frequency variations may also still be present. These are handled by band-pass filtering both input and output signals. By applying the filter to both input and output data no phase information is lost. The upper limit of the band-pass is set to 2 Hz to exclude higher frequency noise and focus the identification data on the frequency band defined by the input signal. The resulting power spectral density of the input signal from one experiment is displayed in Figure 4.2.

10−2 ₁₀−1 ₁₀0 ₁₀1 ₁₀2 −140 −120 −100 −80 −60 −40 −20 0 20 40 Frequency (Hz) Magnitude (dB)

4.2

Yaw dynamics

The first helicopter dynamics to be investigated are the yaw dynamics. As described in Chapter 3, this subsystem can be modelled by relatively simple first order dy-namics. Considering this, a first order time-delayed model is estimated with the Matlab command procest to be used as reference. It should be noted that this model can not be considered as a description of the true dynamics. However, it does give an indication of the expected results from a parametric method of system identification.

4.2.1 Frequency-domain methods

Estimation of the FRF is limited by the excitation in the system caused by the input signal. As there is essentially no information present in the available data in frequency ranges outside the spectrum given by the frequency sweep, the estimated FRF is limited to this band between 0.2 Hz to 2 Hz.

In Figure 4.3 the results of the frequency-domain methods are presented in a Bode diagram. It can be noted that a high degree of consistency is seen amongst the frequency-domain methods. In the amplitude plot, a slight deviation occurs at the limits of the frequency spectrum, which is to be expected considering the transient errors and limited availability of data in these frequencies. The phase curves deviate minimally, at 1 Hz the maximum difference is 2 degrees between two identified models. 0 10 20 Magnitude (dB) procest TRIMM LPM spafdr 0.2 1 2 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

When comparing the nonparametric methods to the parametric procest method it can be seen that the nonparametric models appear to identify the pole present in the system at a higher frequency. This may either be a consequence of variance in the estimation or that actual dynamics are better captured by nonparametric methods. For a first order system with the pole located outside the excited frequency range a linear slope would typically be seen, similar to the procest curve.

4.2.2 Time-domain methods

In the time-domain, the nonparametric method of system identification used is the estimation of a regularized FIR as described in Chapter 2. This method has been implemented in later Matlab versions as the impulseest command. In this thesis the “Tuned/Correlated” regularization kernel has been used for this method [2]. Apart from the choice of regularization kernel, a second choice is whether to use a base-line AR model or not. Here, both types of models are estimated and compared. An estimated FIR model will be denoted FIRARna, where na corresponds to the

order of the base-line model and na = 0 implies that the model was estimated

without a base-line model.

The previously estimated procest model suggests that the impulse response of the system has subsided in about 10 seconds. Considering the sample time of 0.02 s, an order of n = 500 is chosen to capture the full impulse response. The resulting fit to validation data, as caculated by the compare command, for the three models is summarized in the table below:

Model procest FIRAR0 FIRAR5

Fit 83.20 81.15 79.28

It can be noted that all methods have a similar degree of explanation to the valida-tion data of about 80%. However, the actual impulse responses actually vary greatly although the accuracy of the models appear similar. This is a first indication that the frequencies excited in the experiment data are insufficient to fully model the impulse response. It is also clear that the exact appearance of the impulse response is not of great importance when modelling behaviour in this frequency range.

In Figure 4.4 the impulse responses of the procest model and the FIRAR5

model are illustrated. It is immediately visible that there is a great uncertainty in the estimated response, as indicated by the confidence region. The first samples of the FIRAR5 response show a non-minimum phase behaviour. It is clear from

the experiment that such behaviour does not present itself in practice. This result is attributed to the time-delay present in the system, as time-delays are typically approximated by this type of behaviour [14].

The FIRAR5 impulse response appears to follow the response of the procest

model reasonably well. Although there is uncertainty present, this indicates that such a response can indeed explain the characteristics present in data.

In Figure 4.5 the impulse response of the FIRAR0 model is presented. This

0 1 2 3 4 5 6 7 8 9 10 −40 −20 0 20 40 Time (seconds) Amplitude

Figure 4.4. Impulse response of FIRAR5. The procest model is represented by the solid line and the FIRAR5model by the circle markers. The shaded region indicates a confidence interval of three standard deviations.

0 1 2 3 4 5 6 7 8 9 10 −5 0 5 10 15 20 Time (seconds) Amplitude

why there was such a great uncertainty in the previous response. This response is very different, compared to the response of the FIRAR5 model. The size of the

response indicates very different static properties of this model. The beginning of the response is more similar to the FIRAR0 and procest models.

It is interesting to note that this response appears much less uncertain than that of the previously identified model. It is more likely that this would be a consequence of how the method has been implemented, rather than that this model is more accurate.

The regularized FIR method of identification allows for some estimation of sys-tem properties even outside the region of excited frequencies. It is thus intreresting to look at the Bode diagrams of these models in an extended frequency range. In Figure 4.6 and Figure 4.7 the Bode diagrams of the FIRAR5 and FIRAR0 models

are presented, respectively. The procest model is included for reference.

It can be seen in Figure 4.6 that for lower frequencies, the FIRAR5model tends to

first order characteristics. However, at about 1.5 Hz problems arise and the response deviates. This is somewhat unexpected, as the regularization should penalize such deviations. This appears to be a limitation in the estimation method and will be investigated further in Section 4.3.

While the FIRAR5 model displayed consistent behaviour for lower frequencies,

the opposite can be seen for the FIRAR0 model in Figure 4.7. Here, the static

properties of the response disagree greatly with the expectation of the procest model and higher order characteristics are present in the response. This is in line

−20 0 20 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

−20 0 20 40 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

Figure 4.7. Frequency response of FIRAR0. The solid line shows the frequency response of the procest model and the dashed line is the FIRAR0model. The shaded region indicates a confidence interval of three standard deviations.

with what was seen in the impulse response. However, while there is some variance present in the response at higher frequencies, the amplitude curve is suppressed here, contrary to what could be seen in the FIRAR5model.

These results, along with those seen in the previous section of the FRF identi-fication methods, indicate a high dependency on the input signal. It appears that a model of reasonable accuracy can be estimated in the frequency range given by the input signal spectrum. This is further strengthened by the fact that a good fit is seen on validation data. The validation data is the same type of experiment as the one which the models has been estimated from. To attempt to explain the results seen in this section, some results gained from tests on a simulated system are presented in the next section.

4.3

Simulated first order time-delayed system

the nominal system

G0(s) =

K s+ 2πfe

−Tds (4.1)

is used. The parameters are chosen as K = 10, f = 0.1 Hz and Td= 0.1 s. Thus

the system has a stable pole in 0.1 Hz and a time delay of 100 ms.

The system is simulated in open loop using the same type of input signal used in the experiments. This signal u excites frequencies in the band 0.2 − 2 Hz during a 45 s test period. To attempt to model the behaviour seen in real data, the pole is specifically chosen to be outside the range of frequencies excited by the input. In the following simulations no noise on the output is generated.

4.3.1 Frequency-domain methods

When applying the frequency-domain methods to flight data, the responses of the nonparametric models differed from the procest response. One expectation of using nonparametric frequency domain methods is that they should agree with one another to high degree. This due to the fact that they are derived on the same principle, namely the relation of the Fourier transformed output and input data.

The frequency responses of the investigated methods are illustrated in Figure 4.8 together with the response of G0. It can be seen that the estimated responses follow

the true response almost perfectly and with good certainty. The visible deviations from the true response bear some resemblance to the results seen in the experimental data of the previous section. This indicates that similar difficulties as encountered in flight data may be present in the simulated data set.

4.3.2 Time-domain methods

Great variance could be seen in the impulse responses estimated from the flight data. Here, similar steps are taken to estimate FIR models of the simulated system described by (4.1). For reference, a parametric model is estimated wtih the Matlab command procest and the initial guess that the system is a first order system with a time delay. As expected, this model perfectly explains the behaviour of the true system.

Two separate FIR models are generated for comparison. Both models are esti-mated with n = 250 which implies that at most a 5 s long impulse response can be generated. One model, denoted by FIRAR5, is estimated with a fifth order

base-line AR part while the second model, denoted by FIRAR0, is estimated without a

base-line model. The fit to validation data, as calculated by the Matlab compare command, of these models is summarized in the table below:

Model procest FIRAR0 FIRAR5

Fit 99.99 98.76 95.38

−5 0 5 10 15 20 Magnitude (dB) G0 TRIMM LPM spafdr 0.2 1 2 -135 -90 -45 0 Frequency (Hz) Phase (deg)

Figure 4.8. Results of the estimated FRF models.

impulse response of FIRAR5 in Figure 4.9 it can be seen that the accuracy of this

model is limited only by the choice of order. The impulse response of FIRAR0 is

displayed in Figure 4.10. This response differs greatly from the response of FIRAR5

and displays some characteristics also visible in the FIR models estimated on real data. A non-minimum phase characteristic is present in this model as well as some higher frequency gain. From the confidence interval it is clear that the impulse response can not be estimated to any great certainty.

The difference in the two responses can be fully explained by the presence of the AR base-line model. As described in Chapter 2.5, the base-line model is a parametric model structure. In this case this base-line model can fully explain the system behaviour, as expected considering how well the procest model performed. This leaves no residual output to be explained by a FIR model, i.e. the FIR part of the estimation is superfluous.

The frequency response of FIRAR0is displayed in Figure 4.11. The true system is

0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 Time (seconds) Amplitude

Figure 4.9. Impulse response of FIRAR5. The true system is represented by the solid line and the estimated system is by the circle markers.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 −200 −150 −100 −50 0 50 100 150 200 Time (seconds) Amplitude

−20 0 20 40 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

Figure 4.11. Frequency response of FIRAR0. The solid line shows the frequency response of the true system and the dashed line is the estimated system. The shaded region indicates a confidence interval of three standard deviations.

not covered by this interval. Similar results could be seen in the estimates of real data. This is contrary to what is expected, as the concept of regularization implies a correlation among the impulse response coefficients. Thus, the expectation is that the regularization should suppress the response at higher frequencies, in line with the behaviour seen in the lower frequency range.

4.3.3 Analysis

The simulations done on the known system G0 display some of the same behaviour

seen in the real data. The hypothesis that problems with nonparametric identifica-tion methods arise when an insufficient range of frequencies is excited is tested by choosing G0 to have a slow pole outside the frequency range excited by the input

signal.

However, the regularization present in the time-domain methods is intended to allow for some extrapolation over other frequencies. As verified in the simulations, a regularized FIR model without an AR base-line model is not capable of estimating the system in other frequencies other than those excited by the input. This needs to be taken into consideration when evaluating models estimated by nonparametric methods.

4.4

Cyclic dynamics

The dynamics of the main rotor offer different properties compared to the yaw dynamics previously investigated. This subsystem is considered as a two input, two output MIMO system, which is the simplest case of a MIMO system. All the investigated methods of nonparametric system identification can be applied in a MIMO scenario. As discussed in Chapter 3, the main rotor dynamics can be extended to include several other inputs and outputs such as collective pitch δcol

and translational velocities u, v and w. However, this thesis is limited to roll rate

p, pitch rate q and the cyclic inputs δlat and δlon.

Due to how the helicopter fuselage is shaped, the system is not symmetric de-spite the actuator dynamics being the same for both roll and pitch channels. This asymmetric behaviour is especially visible in the roll channel where the system is less damped. The implication of this is that higher frequency pitch excitation cause roll movement as well. This behaviour is illustrated in Figure 4.12. A frequency

−10 −5 0 5 10 q (deg/s) 0 5 10 15 20 25 30 35 40 45 −10 0 10 Time (seconds) p (deg/s)

sweep in the δloncommand has been conducted and it is clearly visible how the roll

channel is excited in the second part of the experiment. The greater amplitude seen in the roll rate implies that although a pitch sweep is conducted, the helicopter ends up with more roll movement than pitch movement.

While the roll dynamics are excited by a pitch command, the opposite does not apply to the same extent. The low excitation in the δlat to pitch rate q channel

has the implication that identification of this cross coupling is difficult and results become uncertain. Therefore analysis of this channel will be omitted in the following sections.

To understand which results are to be expected from a parametric system iden-tification method a model has been estimated as reference. This model has been chosen to be a 4th order model estimated with the Matlab command ssest. This is a black-box model with all parameters free. The delay of this model is estimated by shifting the output data and choosing the model with best fit corresponding to the given delay. By this method the delay has been estimated to 5 samples or 0.1 seconds in both input channels. No shifting of data is done when applying the other methods of identification, as an important aspect is to test how well nonparametric system identification account for time-delays.

4.4.1 Frequency-domain methods

The results are presented in a similar way as the yaw dynamics with a Bode diagram for each interesting input-output channel. The frequency respone of the δlonto pitch

rate q channel is illustrated in Figure 4.13. It is apparent that uncertainty exists in the parametric model as indicated by the confidence interval. The nonparametric methods of identification tend to display similar behaviour as the reference model but there is some disagreement in the magnitude plots. The LPM and TRIMM models do not appear as smooth as the spafdr model, indicating a greater vari-ance. However, all models appear to be in general agreement of the main system properties.

The results across the other channels show a similar or greater level of agree-ment. The Bode diagrams of the δlon to roll rate p and δlat to roll rate p channels

are displayed in Figure 4.14 and Figure 4.15, respectively. Here, a second order characteristic can be identified in the roll rate p output by the shallow resonance peak at 0.5 Hz. This peak is present from both input commands. While such be-haviour is expected, it is clear that by these methods, it can not be identified with any great certainty.

4.4.2 Time-domain methods

−20 −10 0 10 Magnitude (dB) ssest LPM TRIMM spafdr 0.2 1 2 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

Figure 4.13. Results of the estimated FRF models in the δlonto pitch rate q channel.

−5 0 5 10 15 Magnitude (dB) ssest LPM TRIMM spafdr 0.2 1 2 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

0 5 10 15 Magnitude (dB) ssest LPM TRIMM spafdr 0.2 1 2 −180 −135 −90 −45 0 Frequency (Hz) Phase (deg)

Figure 4.15. Results of the estimated FRF models in the δlatto roll rate p channel.

impulse response. The following models have been estimated with an order of

n= 220.

As separate experiments are carried out in the pitch and roll channels, it is interesting to validate the identified models on both a pitch and a roll experiment. In the table shown below, the calculated fit of the identified models is displayed when validated against a roll test.

Model ssest FIRAR0 FIRAR5

Fit in pitch rate q 29.45 48.47 −6.57 Fit in roll rate p 85.47 87.09 85.99

It is clear that none of the models are capable of producing a fit better than 50 % in the pitch channel. This is likely due to the low excitation δlat commands cause

in pitch. However, the fit to roll data can be considered to be very good for all models.

In the next table the calculated fit on a pitch test is displayed. In this case the fit in the pitch channel is significantly better due to the higher level of excitation of the δlon signal. Judging by the fit seen in the roll channel, it is clear that roll

dynamics are well modeled for both δlat and δlon excitation.

Model ssest FIRAR0 FIRAR5

It appears that the use of a base-line model when estimating a regularized FIR has a large impact on the modelling of the pitch channel, where the fit differs greatly from the other models.

The differences between the regularized FIR models estimated with and without base-line model can be further investigated by examining the impulse responses. In Figure 4.16 the impulse response of the δlon to pitch rate q channel is displayed.

A high level of agreement can be seen between the ssest model and the FIRAR0

model. The impulse response shows first order characteristics with a time delay present. Similar non-minimum phase behaviour as seen in the yaw dynamics is present in the beginning of the response. The length of the undershoot can be interpreted as the time-delay present in the control loop.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 −2 −1 0 1 2 3 4 5 Time (seconds) Amplitude

Figure 4.16. Impulse response of FIRAR0 in the δlonto pitch rate q channel. The ssest model is represented by the solid line and the FIRAR0 model by the circle markers.

Interestingly, the impulse response of the FIRAR5model shown in Figure 4.17 is

completely different. This implies that the addition of a base-line AR model causes issues in the estimation of impulse response coefficients. This is also reflected in the lower fit values produced by this model. A reasonable explanation for this result is that the low-order AR model is insufficient when the time-delay is present. A model estimated where data is shifted nk = 3 samples to take the time-delay into account

confirms this assumption, as the impulse response has an appearance like that of the FIRAR0 model. In the remaining input-output channels of the FIRAR5 model the

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 ·10 7 Time (seconds) Amplitude

Figure 4.17. Impulse response of FIRAR5 in the δlonto pitch rate q channel. The ssest model is represented by the solid line and the FIRAR5 model by the circle markers.

From the high fit seen in the roll rate p output of this model it can be concluded that the properties of the impulse response does not have a great impact on excitation in the frequency band defined by the input.

The Bode diagrams of the FIR models are presented in Figure 4.18 and Fig-ure 4.19. As expected from the impulse response plots, the FIRAR0 very closely

resemble the ssest model. The FIRAR5appears much more uncertain and displays

the high gain at high frequencies also appearing in estimates of the yaw dynamics. The low fit of the FIRAR5in the pitch rate q output can be attributed to the higher

variance seen in the frequency region of interest.

As expected from the similar model fit values in the roll rate p, the Bode dia-grams of the two models show a more consistent behaviour. Bode diadia-grams of the

δlon to roll rate p channel are displayed in Figure 4.20 and Figure 4.21. Both FIR

−20 0 20 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

Figure 4.18. Frequency response of the δlon to pitch rate q channel. The solid line shows the frequency response of the ssest model and the dashed line is the FIRAR0 model. The shaded region indicates a confidence interval of three standard deviations.

−20 0 20 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

−20 −10 0 10 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

Figure 4.20. Frequency response of the δlon to roll rate p channel. The solid line shows the frequency response of the ssest model and the dashed line is the FIRAR0 model. The shaded region indicates a confidence interval of three standard deviations.

−20 −10 0 10 Magnitude (dB) 10−2 ₁₀−1 ₁₀0 ₁₀1 −180 −90 0 90 180 Frequency (Hz) Phase (deg)

Discussion and conclusion

5.1

Discussion of results

This thesis has examined how current methods of nonparametric system identifica-tion can be applied to helicopter dynamics. A major goal has been to investigate whether the inclusion of nonparametric models in the system identification process at Saab can be of value. Nonparametric system identification is of interest due to the “unbiased” models it provides, in contrast to the grey-box models typically used. If well-performing nonparametric models can be estimated, their characteristcs may be used as an additional means to decide how a parametric model structure should be chosen.

5.1.1 Frequency-domain approach

One difficulty in determining the usefulness of a nonparametric model is the consid-eration of how such a model should be validated. Especially the FRF models suffer from this shortcoming. Results have shown that while all FRF models display a high level of agreement amongst eachother, it is not easy to validate them in the time-domain. A transformation of frequency response coefficients may be done to construct the impulse response, but is cumbersome. However, considering the re-sults seen on the simulated first-order system, where magnitude and phase curves are reconstructed almost perfectly, the methodology appears viable. If there is a particular interest for the system characteristics in the frequency range excited by the input, a FRF model may be used as a complement. Based on the analysis and for the sake of simplicity, spafdr which has a Matlab implementation can be used for this purpose.

5.1.2 Time-domain approach

The time-domain approach to nonparametric system identification considered in this thesis is the regularized FIR estimate. This method has the benefit of the possibility to be validated on data and in that manner be directly compared to an identified grey-box model. Additionally, time-delays may be easily estimated and visualized by impulse and step responses. The use of this method has proved that some tuning must be done to achieve easily interpreted results.

The undamped yaw dynamics require a high-order model to capture the full length of the impulse response. The slow first-order characteristics are better cap-tured by the inclusion of a base-line AR model in the estimation process. However, it should be clearly stated that while a regularized FIR model attempts to explain dynamics in frequencies outside those excited by the input signal, the model can only be considered trustworthy in the excited frequency range. Current Matlab implementation does not appear to display this uncertainty on all occasions though. In the identification results of the more damped main rotor cyclic dynamics, it can be seen that better results are achieved by omitting the base-line AR model. This result implies that there may exist a point where it is detrimental to add a base-line model to the FIR estimate. Exactly where this point of inflection is located is difficult to assess, but as a rule of thumb a base-line model can be included for a system that is considered undamped.

The results of the cyclic dynamics also show that a regularized FIR model pro-vides a better fit to validation data than the parametric ssest model. This suggests that the nonparametric model may model certain characteristics better and could thus potentially be used for comparison with identified parametric models. Again, the frequency spectrum where the nonparametric model can be considered accurate is defined by the input signal.

5.2

Conclusion

In conclusion, there is definite potential for use of nonparametric models in the system identification process at Saab. While nonparametric methods do in fact also provide quite a few parameters to be chosen, they may be set to a value which appear to work well and be left at that. Choice of parameters have a rather small impact on the properties of the estimated models. When analyzing identification results, care must be taken by the user to consider what frequency range is of interest. In this thesis, dynamics have been excited in a very specific range of frequencies. Should behaviour in other frequencies be of interest, experiments must be designed with that in mind.

useful results and the regularization does not always bring the results that could be expected.

5.3

Future work

Nonparametric methods investigated in this thesis can be implemented in the cur-rent set of functions used for system identification at Saab. This allows for a direct comparison with the grey-box model structures currently used for identification pur-poses. As this thesis has focused on a specific type of experiment conducted in a known frequency range it should be an interesting test to perform the analysis on data covering a broader frequency range.

To generalize the results of the nonparametric methods, identification can be performed on larger data sets. The impact varying control systems and weather conditions has on identification results can be minimized if all available relevant data is included in the identification process.