[PPT] - Models for LTI systems LTI system stands for linear time invariant PowerPoint Presentation

SLIDE 1

Models for LTI systems

LTI system stands for linear time invariant system Model describing LTI system using an input-output relation yk = G(q)uk + H(q)ek with output (time) series y = {yk : k = 0, 1, . . . , N} input (time) series u = {uk : k = 0, 1, . . . , N} and disturbances (e.g. uncontrollable input, noise with distribution fe(·), . . . ) e = {ek : k = 0, 1, . . . , N}

Lecture 3: Elementary Models for LTI Systems 1 / 91

SLIDE 2

Models for LTI systems

Forward delay operator q quk = uk+1 q−1uk = uk−1 Transfer functions G(q) and H(q) associate with impulse responses gk and hk G(q) =

∞

k=1

gkq−k H(q) = 1 +

∞

k=1

hkq−k

Lecture 3: Elementary Models for LTI Systems 2 / 91

SLIDE 3

Modeling disturbances for LTI systems

matlab demo

H(q) = 1 +

∞

k=1

hkq−k Example: Butterworth lowpass filter of degree 2 and cutoff frequency 1/10 of the Nyquist frequency, with transfer function H(q) = 0.0201 + 0.0402q−1 + 0.0201q−2 1.0000 − 1.5610q−1 + 0.6414q−2 Consider two cases 1. ek =

with probability 1 − γ

r with probability γ with γ small (0.02) and r normally distributed

2. ek normally distributed

Lecture 3: Elementary Models for LTI Systems 3 / 91

SLIDE 4

Modeling disturbances for LTI systems

matlab demo 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −100 −80 −60 −40 −20 normalized frequency (Ts = 1) magnitude (dB/Hz) 5 10 15 20 25 30 35 −0.05 0.05 0.1 0.15 time (samples) impulse response Lecture 3: Elementary Models for LTI Systems 4 / 91

SLIDE 5

Modeling disturbances for LTI systems

matlab demo

Case 1:

500 1000 1500 2000 −4 −2 2 4 time (samples) 500 1000 1500 2000 −0.4 −0.2 0.2 0.4 0.6 time (samples) −2 2 500 1000 1500 2000 −2 2 500 1000 1500 2000

For? Disturbances such as peaks, trends and offsets1

1For trends and offsets use lowpass filter with very low cutoff frequencies Lecture 3: Elementary Models for LTI Systems 5 / 91

SLIDE 6

Modeling disturbances for LTI systems

matlab demo

Case 2:

500 1000 1500 2000 −4 −2 2 4 time (samples) 500 1000 1500 2000 −1.5 −1 −0.5 0.5 1 1.5 time (samples) −4 −2 2 4 200 400 600 −4 −2 2 4 200 400 600

For? Measurement noise

Lecture 3: Elementary Models for LTI Systems 6 / 91

SLIDE 7

Modeling disturbances for LTI systems

matlab demo

Do we need the noise distribution function fe(·)? Often we work with some statistical characteristics such as the first and second

rder moments

E [ek] =

xfe(x)dx

σ2 = E (ek − E [ek])2 =

(x − E [ek])2fe(x)dx

Often we assume ek to be white noise with mean zero and unknown variance σ.

Lecture 3: Elementary Models for LTI Systems 7 / 91

SLIDE 8

System identification for LTI systems

How to estimate G(q) and H(q)? Parametrized model with unknown (to be estimated) parameters θ yk = G(q, θ)uk + H(q, θ)ek Structure choice for transfer functions G(q) and H(q)? A(q)yk = B(q) F(q)uk + C(q) D(q)ek Very general, simplifications possible, why? Occam’s razor, avoid overfitting, easier to solve, . . .

Lecture 3: Elementary Models for LTI Systems 8 / 91

SLIDE 9

Short overview of commonly used models

“equation error” (ARX) model structure Input-output relation described by a linear difference equation yk + a1yk−1 + . . . + anayk−na = b1uk + . . . + bnbuk−nb+1 + ek White noise term ek serves to quantify the direct error in the equations. The parameter vector is θ = a1 a2 . . . ana b1 b2 . . . bnb

t

Transfer function A(q)yk = B(q)uk + ek with A(q) = 1 + a1q−1 + . . . + anaq−na B(q) = b1 + . . . + bnbq−nb+1

r with G(q) = B(q)

A(q) and H(q) = 1 A(q):

yk = B(q) A(q) uk + 1 A(q)ek ARX models possess an autoregressive (AR) part A(q)yk, ‘X’ refers to the input term B(q)uk (eXogeneous variable)

Lecture 3: Elementary Models for LTI Systems 9 / 91

SLIDE 10

Short overview of commonly used models

ARMAX model structure ARX models do not adequately describe disturbances, a moving average (MA) model for the equation error provides more flexibility: yk + a1yk−1 + . . . + anayk−na = b1uk + . . . + bnbuk−nb+1 + ek + c1ek−1 + cnc ek−nc

moving average

where ek is white noise. The parameter vector is θ = a1 . . . ana b1 . . . bnb c1 . . . cnc

t

Now, C(q) = 1 + c1q−1 + . . . + cnc q−nc , and yk = B(q) A(q) uk + C(q) A(q) ek G(q) = B(q) A(q) H(q) = C(q) A(q)

Lecture 3: Elementary Models for LTI Systems 10 / 91

SLIDE 11

Short overview of commonly used models

“output error” (OE) model structure Transfer functions G(q) and H(q) in ARX and ARMAX models share the same denominator A(q). From a physical point of view, it would more sense to decouple and parameterize the functions separately. In “output error” models we model relation the relation between input and

utput as a linear difference equation, and assume additive white noise is added:

yk = B(q) F(q)uk + ek Note that this is an ARMAX model where A(q) = C(q) = F(q)

Lecture 3: Elementary Models for LTI Systems 11 / 91

SLIDE 12

Short overview of commonly used models

Box-Jenkins model structure Just as ARX models have been extended to ARMAX models by modeling the “equation error” more extensively, so can OE models be extended to yk = B(q) F(q)uk + C(q) D(q)ek These are called Box Jenkins models. General model structure A(q)yk = B(q) F(q)uk + C(q) D(q)ek While this model unifies all previous ones, we often use the more compact special cases, for which specialized and more efficient algorithms exist. Important: ARMAX, OE and Box-Jenkins models are nonlinear in the parameters!

Lecture 3: Elementary Models for LTI Systems 12 / 91

SLIDE 13

AutoRegressive modeling (AR)

Autoregressive model for a sequence of scalar data yk : k = 0, 1, . . . , N: yk + a1yk−1 + . . . + anayk−na = ek Assumption: sequence ek : k = 0, 1, . . . , N is zero mean white noise E [ek] = 0 ∀k E [ekel] = δklσ2 ∀k, l

◮ AR-model can be considered a whitening filter, converting yk into white

noise ek

◮ or yk generated by internal feedback mechanism combined with white noise

−a1yk−1 − . . . − anayk−na = yk − ek

Lecture 3: Elementary Models for LTI Systems 13 / 91

SLIDE 14

AR transfer function

Transfer function given by (all-pole representation, G(q) = 0) H(q) = 1 A(q, a) = 1 1 + a1q−1 + . . . + anaq−na where θ = a. Solve using linear system

 

yna−1 yna−2 . . . y1 y0 yna yna−1 . . . y2 y1 . . . . . . yN−1 yN−2 . . . yN−na+1 yN−na

    

−a1 −a2 . . . −ana−1 −ana

   =   

yna yna+1 . . . yN−1 yN

   −   

ena ena+1 . . . eN−1 eN

  

Two problems

◮ Consistency not guaranteed ⇒ least squares solution ◮ Real time measurements yk allow updates to a ⇒ recursive least squares

solution

Lecture 3: Elementary Models for LTI Systems 14 / 91

SLIDE 15

AR error covariance

No knowledge of residuals e and noise variance σ2 known beforehand, but can be estimated using E ete = (N − na + 1)σ2 (1) Solve the least squares problem min

ver x inRnb = Ax2

(2) In the assumption that rank (A) = na, we find ˆ x = (AtA)−1Atb (3) Estimate of the variance can be obtained from the norm of the residuals σ2

est = b − A(AtA)−1Atb2 2/(N − na + 1)

(4) while the error covariance of the estimate x is given by E (x − ˆ x)(x − ˆ x)t = σ2(AtA)−1 (5)

Lecture 3: Elementary Models for LTI Systems 15 / 91

SLIDE 16

AR modeling

matlab demo

How to interpret vector of residuals e? Vector of residuals e considered as a perturbation of the right hand side may not be interpreted as measurement noise: noise data yk in the right hand side side also appears as noisefree data in the matrix A:

 

¯ yna−1 ¯ yna−2 . . . ¯ y1 ¯ y0 ¯ yna ¯ yna−1 . . . ¯ y2 ¯ y1 . . . . . . ¯ yN−1 ¯ yN−2 . . . ¯ yN−na+1 ¯ yN−na

    

−a1 −a2 . . . −ana−1 −ana

   =   

¯ yna ¯ yna+1 . . . ¯ yN−1 ¯ yN

  

The residuals should be considered as an ‘input’ to the system instead, consider two signals described by

1. Data yk are observed exactly without noise but are generated by driving a

white noise sequence through an all-pole system yk + a1yk−1 + . . . + anayk−na = ek

2. Data yk are generated by an autonomous linear deterministic system but
bserved with measurement errors ek

yk + a1yk−1 + . . . + anayk−na = 0 ¯ yk = yk + ek

Lecture 3: Elementary Models for LTI Systems 16 / 91

SLIDE 17

AR modeling

matlab demo 50 100 150 200 250 300 −5 5 50 100 150 200 250 300 −5 5 Lecture 3: Elementary Models for LTI Systems 17 / 91

SLIDE 18

Asymptotic properties

What if we increase the number of data samples as opposed to the number of experiments?2 x − ˆ x = (AtA)−1Ate As m → ∞, the number of rows in A increases, when ek is a white noise sequence lim

m→∞

Ate m = 0 Therefore lim

m→∞ x = ˆ

x Example: yk − 2.4yk−1 + 2.7225yk−2 − 1.4865yk−3 + 0.3741yk−4 = ek with poles z = 0.7 ± 0.65i and z = 0.5 ± 0.4i.

2Keep in mind that in many practical applications only a limited amount of experiments is

allowed!

Lecture 3: Elementary Models for LTI Systems 18 / 91

SLIDE 19

Asymptotic properties

matlab demo 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −100 −50 50 100 frequency

ne−sided PSD (dB/Hz)

PSD estimate of x PSD of model output 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −60 −40 −20 20 40 frequency

ne−sided PSD (dB/Hz)

Figure : (above) Fourier transform of output data together with magnitude of identified transfer function. (below) Fourier transform of residuals. These are clearly white (flat spectrum).

Lecture 3: Elementary Models for LTI Systems 19 / 91

SLIDE 20

Asymptotic properties

matlab demo 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0.7 0.8 0.9 1 1.1 1.2 variance estimate of white noise sequence

Figure : Estimate σest of the variance σ of the white noise sequence e for an increasing number of samples N.

Lecture 3: Elementary Models for LTI Systems 20 / 91

SLIDE 21

Asymptotic properties

matlab demo 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 −2.6 −2.4 −2.2 estimate and error bounds (66%) for first coefficient 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 2.2 2.4 2.6 2.8 3 3.2 estimate and error bounds (66%) for second coefficient

Figure : Estimates and 66 % error bounds of a1 (above) and a2 (below) for an increasing number of samples N.

Lecture 3: Elementary Models for LTI Systems 21 / 91

SLIDE 22

Asymptotic properties

matlab demo

500 1000 1500 2000 2500 3000 3500 4000 4500 5000 −2 −1.8 −1.6 −1.4 −1.2 −1 estimate and error bounds (66%) for third coefficient 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0.1 0.2 0.3 0.4 0.5 0.6 estimate and error bounds (66%) for fourth coefficient

Figure : Estimates and 66 % error bounds of a3 (above) and a4 (below) for an increasing number of samples N.

Lecture 3: Elementary Models for LTI Systems 22 / 91

SLIDE 23

Asymptotic properties

What if ek is colored? First element of limm→∞(Ate)/m is lim

m→∞

1 m

N−1

k=na−1

ykek−na+1 If ek is a white noise sequence, E [ykek−n+1] = 0, if not yk = ...some function of previous yk... + ek E [ekek−na+1] = 0 Generalized to all coefficients, we get that limm→∞(Ate)/m = 0, or lim

m→∞ x = ˆ

x In addition, the error covariance matrix decays to zero as 1/m, implying better estimates as m → ∞

Lecture 3: Elementary Models for LTI Systems 23 / 91

SLIDE 24

AutoRegressive modeling with eXogeneous input (ARX)

An AutoRegressive model with eXogeneous input can be described by yk + a1yk−1 + . . . + anayk−na = b1uk + . . . + bnbuk−nb + ek Assumption: sequence ek : k = 0, 1, . . . , N is zero mean white noise with standard deviation σ

◮ Easy to calculate (least squares) ◮ Fairly good results even with violated assumptions ◮ Inaccurate noise models: 1/A(q), pushes towards high order ARX models

What about the model orders? Not a trivial problem. Try out different combinations of model orders and inspect the errors in function of their model orders.

Lecture 3: Elementary Models for LTI Systems 24 / 91

SLIDE 25

ARX transfer function

A(q)yk = B(q)uk + ek Because of white noise assumptions we obtain the linear least squares problem

 

yna−1 yna−2 . . . y1 y0 una . . . una−nb+2 una−nb+1 yna yna−1 . . . y2 y1 una+1 . . . una−nb+3 una−nb+2 . . . . . . . . . . . . yN−1 yN−2 . . . yN−na+1 yN−na uN . . . uN−nb+2 uN−nb+1

          

−a1 −a2 . . . −ana−1 −ana b1 . . . bnb−1 bnb

        

=

  

yna yna+1 . . . yN−1 yN

   −   

ena ena+1 . . . eN−1 eN

  

(6) Lecture 3: Elementary Models for LTI Systems 25 / 91

SLIDE 26

ARX properties

Error covariance of the estimate is given by E (x − ˆ x)(x − ˆ x)t = σ2AtA−1 Asymptotic properties are similar to AR model, notice that one can validate the ARX model based upon input and error correlation E [ukek−n+1] = 0 In case of a colored noise sequence, the result will be biased, even for an infinite amount of data.

Lecture 3: Elementary Models for LTI Systems 26 / 91

SLIDE 27

ARX properties and unbiased estimates

Special case 1: What if we have information on the noise properties? yk = B(q) A(q) uk + κ(q) A(q)ek Filtering sequences uk and yk using inverse of κ(q) uF

k = κ−1(q)uk

y F

k = κ−1(q)yk

gives the ARX identification problem y F

k = B(q)

A(q)uF

k +

1 A(q)ek

Lecture 3: Elementary Models for LTI Systems 27 / 91

SLIDE 28

ARX properties and unbiased estimates

matlab demo

Consider the fourth order ARX system G(q) = 1 − 1.0011q−1 + 1.3699q−2 − 1.1227q−3 1 + 0.8060q−4 H(q) = 1 1 + 0.8060q−4 with na = 4 and nb = 4. Two experiments are conducted:

1. 100 Monte Carlo simulations of the output yk are simulated with ek white
noise. From each of the 100 input-output pairs, an ARX model with

na = 4, nb = 4 is determined. From these 100 models, the average transfer function is then determined.

2. The same 100 experiments are repeated (na = 4, nb = 4) but with ek

colored noise (coloring by a fourth order low pass Butterworth filter).

Lecture 3: Elementary Models for LTI Systems 28 / 91

SLIDE 29

ARX properties and unbiased estimates

matlab demo 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −20 −10 10 20 30 40 normalized frequency (Ts = 1) magnitude (dB/Hz) ARX system ARX model (white) ARX model (colored) Lecture 3: Elementary Models for LTI Systems 29 / 91

SLIDE 30

ARX properties and unbiased estimates

Output-error model (OE) yk = B(q) A(q) uk + ek in this case κ(q) = A(q), but A(q) is unknown! Two possibilities

1. rough estimate of A(q)
2. ARX model fitting (A1(q), B1(q)), filter uk, yk with 1/A1(q) and obtain

second ARX model

1

A1(q)yk

= B2(q)

A2(q)

1

A1(q)uk

+

1 A2(q)ek

Lecture 3: Elementary Models for LTI Systems 30 / 91

SLIDE 31

ARX properties and unbiased estimates

Special case 2: noise term well described by high order AR model 1/D(q) (order r) A(q)yk = B(q)uk + 1 D(q)ek Equivalently A(q)D(q)yk = B(q)D(q)uk + ek which is an ARX model with orders na + r and nb + r and transfer function B(q)D(q) A(q)D(q) = B(q) A(q) High order ARX models are often fairly good models. Can we reduce the model order?

1. Estimate a high order ARX model
2. Do model reduction on this model. Most model reduction techniques give

an indication of the number of important modes that have to be retained.

Lecture 3: Elementary Models for LTI Systems 31 / 91

SLIDE 32

Prediction Error Methods (PEM)

Basic input-output equation for a single input-output linear system with additive disturbance yk = G0(q)uk + H0(q)ek Parametrization (e.g. ARX, ARMAX, . . . ) defines model set M M = {G(q, θ), H(q, θ) | θ ∈ Rm×1} Postulate that (G0(q), H0(q)) ∈ M Given a set of input-output data uk, yk, k = 1, . . . , N and a parametrization, find θopt such that the given data is predicted as “wel” as possible

Lecture 3: Elementary Models for LTI Systems 32 / 91

SLIDE 33

Prediction

What is the best prediction ˆ yk with G0(q) and H0(q) using ul, yl, l = 1, . . . , k − 1 and the next input to the system uk? ˆ y pred

k

= H−1

0 (q)G0(q)uk + (1 − H−1 0 (q))yk

= Wuuk + Wyyk Prediction error is the difference between measured and predicted output epred

k

= yk − ˆ y pred

k

= H−1

0 (q)(yk − G0(q)uk)

= ek

Lecture 3: Elementary Models for LTI Systems 33 / 91

SLIDE 34

ARX prediction

G0(q) = b1 + b2q−1 + . . . + bnbq−nb+1 1 + a1q−1 + . . . + anaq−na H0(q) = 1 1 + a1q−1 + . . . + anaq−na Prediction of ARX ˆ y pred

k

= H−1

0 (q)G0(q)uk + (1 − H−1 0 (q))yk

= (b1 + . . . + bnbq−nb+1)uk − (a1q−1 + . . . + anaq−na)yk = b1uk + . . . + bnbuk−nb+1 − a1yk−1 − . . . − anayk−na = yk − ek If ek a white noise sequence, we can compute using linear least squares

Lecture 3: Elementary Models for LTI Systems 34 / 91

SLIDE 35

ARMAX prediction

G0(q) = b1 + b2q−1 + . . . + bnbq−nb+1 1 + a1q−1 + . . . + anaq−na H0(q) = 1 + c1q−1 + . . . + cnc q−nc 1 + a1q−1 + . . . + anaq−na Prediction of ARMAX ˆ y pred

k

= H−1

0 (q)G0(q)uk + (1 − H−1 0 (q))yk

= b1 + b2q−1 + . . . + bnbq−nb+1 1 + c1q−1 + . . . + cnc q−nc uk + (c1 − a1)q−1 + . . . + (nnc − anc ) + . . . − anaq−na 1 + c1q−1 + . . . + cnc q−nc yk Non-linear in the parameters means non-linear optimization problems!

Lecture 3: Elementary Models for LTI Systems 35 / 91

SLIDE 36

Prediction error and model sets

Prediction errors play an important role in search for the optimal θopt, need a way of expressing in terms of parametrization of the system epred

k

(θ) = yk − y pred

k

(θ) = H−1(q, θ)(yk − G(q, θ)uk) = ek For a time series data set u1, u2, . . . , uN and y1, y2, . . . , yN and a specific choice for θ, we get a sequence of prediction errors e1(θ), e2(θ), . . . , eN(θ) Based on the data uk, yk and the parametrization G(q, θ), H(q, θ), we can compute the prediction error ek(θ). Select your estimate of θ (θopt) in such a way that the series of prediction errors is minimized

Lecture 3: Elementary Models for LTI Systems 36 / 91

SLIDE 37

Minimizing prediction errors

Prediction error sequence ek as a vector in RN

◮ many norms available in RN ◮ filtering before applying norm (frequency weighting)

ef

k(θ) = L(q)ek(θ) ◮ norm as an objective function to be minimized

V (θ) = 1 N

N

k=1

l(ef

k(θ))

with l(·) a positive scalar function

◮ optimal parameter θopt

θopt = argmin

θ

V (θ)

Lecture 3: Elementary Models for LTI Systems 37 / 91

SLIDE 38

Minimizing prediction errors

Example: quadratic function used in least squares V (θ) = 1 N

N

k=1

1 2(ef

k(θ))2

Example: l1 norm where l(·) = | · | is more robust against outliers Example: using the l∞ norm, defined as l∞(·) = maxdomain |l(·)|, outliers will determine the resulting model

Lecture 3: Elementary Models for LTI Systems 38 / 91

SLIDE 39

Frequency domain interpretation of the PEM

ef

k(θ) = L(q)(yk − ˆ

y pred

k

(θ)) = L(q)H−1(q, θ)(G0(q) − G(q, θ))uk + L(q)H−1(q, θ)H0(q)ek Using Parseval’s theorem V (θ) = 1 2π

π

−π

1 2|G0(ejω) − G(ejω, θ)|2 |L(ejω)|2|U(ω)|2 |H(ejω)|2 dω + higher order terms (7) where

◮ G0(q) is the unknown (real) transfer function from uk to yk ◮ G(q, θ) and H(q, θ) form the model ◮ U(ω) is the spectrum of the input sequence xk ◮ L(q) the filter used to filter the prediction errors

Lecture 3: Elementary Models for LTI Systems 39 / 91

SLIDE 40

Frequency domain interpretation of the PEM

PEM is equivalent to optimal fitting in frequency domain of G(q, θ) to G0(q), using the weighting function Q(ω) = |L(ejω)|2|U(ω)|2 |H(ejω)|2 Weighting function is signal to noise ratio weighted by |L(ejω)|2

◮ In case of undermodeling, more emphasis will be put on improving the fit

at frequencies where signal to noise ratio is highest.

◮ When the model order is sufficient to model the transfer function (or in

case of overmodeling), more or less emphasis can be put onto certain frequencies or frequency bands. Good fit is obtained in frequency bands that are excited with a high signal to noise ratio and vice versa.

Lecture 3: Elementary Models for LTI Systems 40 / 91

SLIDE 41

Frequency domain interpretation of the PEM

What are the effects of a wrong parameter class? Example: true OE model yk = G0(q)uk + ek (8) to which we fit an ARX model G(q, θ) = B(q, θ) A(q, θ) H(q, θ) = 1 A(q, θ) (9) Frequency fit: V (θ) ≈ 1 2π

π

−π

1 2|G0(ejω) − B(ejω, θ) A(ejω, θ) , θ)|2|L(ejω)|2|U(ω)|2|A(ejω, θ)|2dω (10)

Lecture 3: Elementary Models for LTI Systems 41 / 91

SLIDE 42

Frequency domain interpretation of the PEM

For many practical systems the filter 1/A(q) is a low pass filter, thus A(q) will be a high pass filter. When estimating an ARX model to an OE system, the resulting fit will be good for high frequencies but poor for low frequencies. One solution consists of filtering the data with an estimate of 1 A(q) ≃ 1 ˜ A(q) (11)

r equivalently use a filter L(q) = 1/˜

A(q): V (θ) ≈ 1 2π

π

−π

1 2|G0(ejω) − B(ejω, θ) A(ejω, θ) , θ)|2|U(ω)|2dω (12)

Lecture 3: Elementary Models for LTI Systems 42 / 91

SLIDE 43

Frequency domain interpretation of the PEM

matlab demo

Example: effect of the spectrum of the input signal |U(jω)|2 G0(q) = 1 − 1.0011q−1 + 1.3699q−2 − 1.1227q−3 1 + 0.8060q−4 (13) and no disturbance H0(q) = 0. Identify an OE model of the form G(q, θ) = b1 + b2q−1 + b3q−2 1 + f1q−1 + f2q−2 H(q, θ) = 1 with 4 different input signals applied to the system

1. A white noise signal
2. A white noise signal filtered through a low pass filter
3. A white noise signal filtered through a high pass filter
4. A block wave

Lecture 3: Elementary Models for LTI Systems 43 / 91

SLIDE 44

Frequency domain interpretation of the PEM

matlab demo 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −20 −10 10 20 30 40 50 normalized frequency (Ts = 1) magnitude (dB/Hz) true model white noise input low freq. input high freq.input block wave input Lecture 3: Elementary Models for LTI Systems 44 / 91

SLIDE 45

Frequency domain interpretation of the PEM

matlab demo

Example: effect of the wrong parameterization class Consider the system yk = G0(q)uk (14) with G0(q) = 0.001 10q−2 + 7.4q−3 + 0.924q−4 + 0.1764q−5 1 − 2.0236q−1 + 1.4022q−2 − 0.3932q−3 + 0.0372q−4 (15) The system suffers no disturbances and white noise is used as input or |U(jω)|2 = 1. The filter L(q) is also equal to 1. We generate 1000 samples and generate two different models from the same data (noisefree), an OE and an ARX model: G(q, θ) = b1 + b2q−1 + b3q−2 1 + f1q−1 + f2q−2 G(q, θ) = b1 + b2q−1 + b3q−2 1 + a1q−1 + a2q−2 (16) H(q, θ) = 1 H(q, θ) = 1 1 + a1q−1 + a2q−2 (17)

Lecture 3: Elementary Models for LTI Systems 45 / 91

SLIDE 46

Frequency domain interpretation of the PEM

matlab demo 10

−2

10

−1

10 −90 −80 −70 −60 −50 −40 −30 −20 −10 10 20 frequency (fmax = π) magnitude (dB/Hz)

Figure : True model (full line), OE estimate (dashed line) and error of the estimate (dotted line). OE orders: nb = 3, nf = 2.

Lecture 3: Elementary Models for LTI Systems 46 / 91

SLIDE 47

Frequency domain interpretation of the PEM

matlab demo 10

−2

10

−1

10 −90 −80 −70 −60 −50 −40 −30 −20 −10 10 20 frequency (fmax = π) magnitude (dB/Hz)

Figure : True model (full line), ARX estimate (dashed line) and error of the estimate (dotted line). The error is clearly smaller for high frequencies. The weighting |ˆ A(q)| is also depicted (dashed dotted line). ARX orders: na = 2, nb = 3.

Lecture 3: Elementary Models for LTI Systems 47 / 91

SLIDE 48

Frequency domain interpretation of the PEM

matlab demo 10

−2

10

−1

10 −90 −80 −70 −60 −50 −40 −30 −20 −10 10 20 frequency (fmax = π) magnitude (dB/Hz)

Figure : Second iterate of the ARX estimate, after filtering both input and output data with 1/ˆ A(q) with ˆ A(q) the estimated numerator in the first step. The improvement if clearly visible.

Lecture 3: Elementary Models for LTI Systems 48 / 91

SLIDE 49

Asymptotic distribution

Case 1: true system θ0 in selected model class M θopt = argmin

θ

V (θ) V (θ) = 1 N

N

k=1

1 2(ek(θ))2 If N → ∞ then E [(θopt − θ0)] = 0 E (θopt − θ0)(θopt − θ0)t = Pθ N where Pθ = E ∇2V (θ0)−1 Qθ

E

∇2V (θ0)−1 Qθ = lim

N→∞ NE

(∇V (θ0))(∇V (θ0))t

Lecture 3: Elementary Models for LTI Systems 49 / 91

SLIDE 50

Asymptotic distribution

Example: ARX using least squares V (θ) = 1 2Aθ − b2 Here θ corresponds to x in least squares framework V (θ) = 1 2θtAtAθ − btAθ + 1 2bbt Derivatives ∇V (θ) = At(Aθ − b) ∇2V (θ) = AtA

Lecture 3: Elementary Models for LTI Systems 50 / 91

SLIDE 51

Asymptotic distribution

Let ∇V (θ)|θ0 = e(θ0), the residuals of the least squares equations for θ0 Pθ N = (AtA)−1 Qθ N (AtA)−1 = (AtA)−1 1 N lim

N→∞ NE

Ate(θ0)e(θ0)tA (AtA)−1 = (AtA)−1(σ2AtA)(AtA)−1 = σ2(AtA)−1 Using first order approximation of some f (ˆ θ) and Gauss’ approximation formula, the variance of parameter estimate θopt can be used to compute variances of

◮ transfer function ◮ poles and zeros of model ◮ impulse response

Lecture 3: Elementary Models for LTI Systems 51 / 91

SLIDE 52

Asymptotic distribution

Case 2: true system θ0 not in model class There is no θ ∈ Rm such that G0(q) = G(q, θ) H0(q) = H(q, θ) Parameter estimates θopt will be Gaussian distributed around some ˆ θ0 different from θ0

◮ Cannot compare ˆ

θ0 and θ0 directly (dimension, model, . . . )

◮ Can compare in frequency domain

G0(ejω) − G(ejω, ˆ θ0)

◮ Variance of θopt around ˆ

θ0 can be described

◮ More parameters will reduce the error but increase variance

Lecture 3: Elementary Models for LTI Systems 52 / 91

SLIDE 53

Asymptotic distribution

Total mean square error in frequency domain E |G0(ejω)| − |G(ejω, θopt)|2 Independence of bias and variance E |G0(ejω)| − |G(ejω, θopt)|2 = |G0(ejω)| − |G(ejω, ˆ θ0)|2

bias

+ E |G(ejω, ˆ θ0)| − |G(ejω, θopt)|2

variance

(18) Minimal mean square error by selecting an optimal number of parameters

Lecture 3: Elementary Models for LTI Systems 53 / 91

SLIDE 54

Correlation methods

What is a “good” model? prediction errors ek(θ) independent (i.e. not correlated) with past data uk, uk−1, . . . and yk−1, yk−2, . . .

utput signal yk contains all information from input signal uk and ek is a

remainder

Lecture 3: Elementary Models for LTI Systems 54 / 91

SLIDE 55

Correlation vectors

Choose a vector sequence ξk (correlation vector) derived from uk, uk−1, . . . and yk−1, yk−2, . . . and impose non-correlation constraint3 1 N ξt

ke(θ) = 0

r using frequency weighting

ef

k(θ) = L(q)ek(θ)

to find the optimal parameters θopt = solution

1 N ξt

kef (θ) = 0

3Note that ξ is a matrix of correlation vectors

Lecture 3: Elementary Models for LTI Systems 55 / 91

SLIDE 56

Instrumental variables

Example: ARX model yk = −a1yk−1 − . . . − anayk−na + b1uk + . . . + bnbuk−nb+1 + vk with vk colored noise (unknown covariance) How to obtain consistent estimates? A(q)yk = B(q)uk + vk vk = C(q)ek ARMAX model in disguise, but not interested in C(q)

Lecture 3: Elementary Models for LTI Systems 56 / 91

SLIDE 57

Instrumental variables and linear least squares

ARX model b = Aθ + v with solution θleastsquares = solution

1 N At(Aθ − b) = 0

= solution

1 N Atv(θ) = 0

Linear least squares is an instrumental variable technique with

ξ = A L(q) = 1 Instruments ξ (= A) badly chosen, may be correlated with error vk, expect biased result

Lecture 3: Elementary Models for LTI Systems 57 / 91

SLIDE 58

Instrumental variables and linear least squares

Choose general correlation vectors ξ and formulate the problem θopt = solution

1 N ξt(Aθ − b) = 0

= solution

1 N ξtv(θ) = 0

which can be solved by

ξtb = (ξtA)θopt ⇒ θopt = (ξtA)−1ξtb

Lecture 3: Elementary Models for LTI Systems 58 / 91

SLIDE 59

Instrumental variables and linear least squares

Typical choice for ξ

 

xna−1 xna−2 . . . x1 x0 una . . . una−nb+2 una−nb+1 xna xna−1 . . . x2 x1 una+1 . . . una−nb+3 una−nb+2 . . . . . . . . . . . . xN−1 xN−2 . . . xN−na+1 xN−na uN . . . uN−nb+2 uN−nb+1

 

Variables xk generated from input uk by a linear system N(q)xk = M(q)uk If we knew the real system, then N(q) = A(q), M(q) = B(q), and xk = B(q) A(q) uk = yk − vk so xk is the undisturbed output of the system!

Lecture 3: Elementary Models for LTI Systems 59 / 91

SLIDE 60

Instrumental variables and linear least squares

matlab demo

Closed loop approach

1. Compute A(q, θopt) and B(q, θopt) using linear least squares
2. Use estimates as N(q) and M(q) in the instrumental variables method

Example: estimate ARX and IVX model of A(q)yk = B(q)uk + vk (19) vk = C(q)ek (20) where A(q) = 1 + 0.8060q−4 B(q) = 1 − 1.0011q−1 + 1.3699q−2 − 1.1227q−3 and vk is obtained by filtering white noise through a Butterworth bandpass filter from normalized frequencies 0.2 to 0.4. The estimated ˆ A(q) and ˆ B(q) are used to generate xk from uk

Lecture 3: Elementary Models for LTI Systems 60 / 91

SLIDE 61

Instrumental variables and linear least squares

matlab demo 10

−3

10

−2

10

−1

−20 −15 −10 −5 5 10 15 20 25 30 normalized frequency (Ts = 1) magnitude (dB/Hz) true system ARX model

Figure : True system (full line) and estimated ARX model (dashed line) of orders na = 4, nb = 4.

Lecture 3: Elementary Models for LTI Systems 61 / 91

SLIDE 62

Instrumental variables and linear least squares

matlab demo 10

−3

10

−2

10

−1

−20 −15 −10 −5 5 10 15 20 25 30 normalized frequency (Ts = 1) magnitude (dB/Hz) true system IVX model

Figure : True system (full line) and estimated IVX model (dashed line) from an ARX model of orders na = 4, nb = 4.

Lecture 3: Elementary Models for LTI Systems 62 / 91

SLIDE 63

Preprocessing

Before you start, always clean your data! Disturbances present in raw data

◮ spikes ◮ drifts ◮ offsets ◮ significant differences in power of various input and output signals ◮ time delays in measured transfers

For time delay, there are two approaches

1. More parameters
2. Correct signals for time delay first, then system identification (Occam’s

razor)

Lecture 3: Elementary Models for LTI Systems 63 / 91

SLIDE 64

Trend determination and correction

Not all process variables contributing to changes found at the output may be available in a dataset As a result, these disturbances are considered as colored output noise

◮ Characteristics of contributions are often known ◮ In many cases slow variations of the output: trends

Slow variations imply we only see a period to a few periods of the trend signal in the data interval of a time series dataset

Lecture 3: Elementary Models for LTI Systems 64 / 91

SLIDE 65

Trend determination and correction

Using the instrumental variables method, it is required that

N

k=1

ξkvk = 0 For vk a large disturbance with few periods, uk will be correlated to vk

N

k=1

ukvk = 0 for higher frequency trends, the sum averages out faster Slower disturbances have a larger impact on the estimated model!

Lecture 3: Elementary Models for LTI Systems 65 / 91

SLIDE 66

Trend determination and correction

50 100 150 200 250 300 350 400 450 500 −0.1 −0.05 0.05 0.1 0.15 50 100 150 200 250 300 350 400 450 500 −1 −0.5 0.5 1 Lecture 3: Elementary Models for LTI Systems 66 / 91

SLIDE 67

Detrending

matlab demo

Causal or anti-causal filtering: how to combat phase shift error? ⇒ combine two filters in parallel (causal + anti-causal) and average out

50 100 150 200 250 300 350 400 450 500 −5 5 50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10

Figure : Original signal (above) and signal with trend (below) showing forward and backward filtering (dashed lines) and combination (full black line).

Lecture 3: Elementary Models for LTI Systems 67 / 91

SLIDE 68

Detrending

matlab demo

Causal or anti-causal filtering: how to combat phase shift error? ⇒ combine two filters in parallel (causal + anti-causal) and average out

50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10 50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10

Figure : Original trend (above) and estimated trend (below). The quality of detrending depends a lot on prior knowledge of the frequency band in which the trend

ccurs and its overlap with the original signal.

Lecture 3: Elementary Models for LTI Systems 68 / 91

SLIDE 69

Detrending

◮ Initial conditions of low pass filters important to separate trend as

accurately as possible

◮ No applications possible to online processing ◮ Highpass filtering to data, but better to focus only on output signals ◮ Finish with offset correction for input and output signals separately, most

industrial processes have high offset levels and small variations Alternative detrending schemes

◮ Constant/linear detrending using detrend ◮ Polynomial detrending using polyfit ◮ Locally weighted scatterplot smoothing (LOESS) using smooth

Lecture 3: Elementary Models for LTI Systems 69 / 91

SLIDE 70

Detrending

Polynomial fitting on centered and normalized time series

50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10 50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10

Figure : Original trend (above) and estimated trend (below). The fitted polynomial is

f order 10.

Lecture 3: Elementary Models for LTI Systems 70 / 91

SLIDE 71

Detrending

Local regression using weighted linear least squares and a 2nd degree polynomial model (LOESS)

50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10 50 100 150 200 250 300 350 400 450 500 −15 −10 −5 5 10

Figure : Original trend (above) and estimated trend (below). The span of the LOESS algorithm has been set to 50 data points.

Lecture 3: Elementary Models for LTI Systems 71 / 91

SLIDE 72

Peak shaving

Disturbance of measures signals with spikes introduced in sensors propagating to measuring equipment

◮ Amplitudes of spikes large compared to process signal ◮ Typically last one to tens of samples ◮ Make up large chunk of the noise energy

Removing spikes by ‘peak shaving’ using four step method

1. Clip signal amplitudes to values never reached by the real process signals
2. Compute a trend signal of the clipped signal
3. Compute a standard deviation of the trend-corrected, clipped signal
4. Interpolate all samples or the original signal that are outside a band

defined by the trend signal

Lecture 3: Elementary Models for LTI Systems 72 / 91

SLIDE 73

Peak shaving

matlab demo 50 100 150 200 250 300 350 400 450 500 −5 5 50 100 150 200 250 300 350 400 450 500 −10 10 50 100 150 200 250 300 350 400 450 500 −20 20 50 100 150 200 250 300 350 400 450 500 −20 20

Figure : Test setup: signal from ARX model (top) is corrupted with a trend signal and peaks into raw data (bottom).

Lecture 3: Elementary Models for LTI Systems 73 / 91

SLIDE 74

Peak shaving

matlab demo

1 Clip signal amplitudes to values never reached by the real process signals ˜ sk = f cl

k sk

with ˜ sk the clipped signal, sk the measured signal and f cl

k =

smax(1/sk)

sk ≥ smax 1 smin < sk < smax smin(1/sk) sk ≤ smin 2 Compute a trend signal ¯ sk of the clipped signal ˜ sk

Lecture 3: Elementary Models for LTI Systems 74 / 91

SLIDE 75

Peak shaving

matlab demo 50 100 150 200 250 300 350 400 450 500 −20 20 50 100 150 200 250 300 350 400 450 500 −5 5 50 100 150 200 250 300 350 400 450 500 −5 5

Figure : The clipped version (middle) of the starting signal (top) preserves values between the range [−7.0, 5.5]. After clipping the trend is determined (bottom).

Lecture 3: Elementary Models for LTI Systems 75 / 91

SLIDE 76

Peak shaving

matlab demo

3 Compute a standard deviation of the trend-corrected, clipped signal σ =

N
k=1

(˜ sk − ¯ sk)2 3 Interpolate all samples or the original signal that are outside a band defined by the trend signal plus and minus α times the standard deviation obtained from the third step. The permitted signal is given by xk =

¯

sk + ασ upper limit ¯ sk − ασ lower limit All consecutive values outside the permitted band are replaced by values from a linear interpolation from the last accepted value to the first sample value within the permitted band after the spike Typically, α = 3

Lecture 3: Elementary Models for LTI Systems 76 / 91

SLIDE 77

Peak shaving

matlab demo 50 100 150 200 250 300 350 400 450 500 −5 5 50 100 150 200 250 300 350 400 450 500 −10 10 50 100 150 200 250 300 350 400 450 500 −10 10

Figure : Based on the trend of the clipped signal (top), a permitted band (middle) is determined with α = 3.0. Values of the original signal not lying in the band are interpolated (bottom).

Lecture 3: Elementary Models for LTI Systems 77 / 91

SLIDE 78

Peak shaving

matlab demo 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −10 −5 5 10 15 normalized frequency (Ts = 1) magnitude (dB/Hz) true system estimated model

Figure : The peak shaved signal is detrended and used to estimate an ARX model (dashed line) of the same order as the true system (full line).

Lecture 3: Elementary Models for LTI Systems 78 / 91

SLIDE 79

Delay estimation

◮ If signals are not corrected for time delays before system identification,

they become part of the model (more parameters)

◮ Time delays do not contribute to the further description of the dynamic

behavior of the process

◮ Signals can be compensated if identification is performed offline

Estimation can be carried out using correlation techniques Main idea is that cross-correlation values start to differ from zero when the shift between input and output signals reaches the sought after time delay

Lecture 3: Elementary Models for LTI Systems 79 / 91

SLIDE 80

Delay estimation

Input signal: stationary, white, (inter-channel independent for MIMO), zero mean white noise sequence with standard deviation σ E [ukul] = σ2δkl Process response to an arbitrary input signal as a convolution using impulse reponse hk yk =

∞

i=0

hiuk−i + nk Here we assume output noise signal uncorrelated with input signal Cross-correlation between process input and output ψyu(τ) = lim

N→∞

1 N + 1

N

i=0

yi+rui

Lecture 3: Elementary Models for LTI Systems 80 / 91

SLIDE 81

Delay estimation

Substituting convolution expression gives ψyu(τ) =

∞

j=0

hjψuu(j − τ) Here, ψuu is the autocorrelation function of the input and an impulse at zero if the assumption of white noise is valid, hence ψuu(τ) = σ2hr Cross-correlation will have shape of the impulse response, the time delay can be found by looking for the beginning of a sequence of values that differ significantly from zero

Lecture 3: Elementary Models for LTI Systems 81 / 91

SLIDE 82

Delay estimation

Confidence bands upper = 3σ[ψyu(τ < 0)] lower = −3σ[ψyu(τ > 0)] If the time delay is ambiguous, it’s safer to underestimate and extra parameters during modeling to prevent incorrect models

Lecture 3: Elementary Models for LTI Systems 82 / 91

SLIDE 83

Filtering

Filtering has to prohibit aliasing effects and reduce influence of disturbances on the measured signals Importance filter design: prevent loss of any relevant process data

1. Filter out all irrelevant information during recording of process data
2. Record extra information of process and use to increase ratio relevant

process information and blurring disturbances

Lecture 3: Elementary Models for LTI Systems 83 / 91

SLIDE 84

Filtering

Filters must be designed so that they will only pass information in the frequency range that is relevant for the identification of the system dynamics + anti aliasing

◮ In the pass band, filters should not be allowed to have noticeable influence

n the signal characteristics (flat band pass)

◮ Cut-off frequency of filters must be chosen high enough such that no

damping is introduced at frequencies covered by process dynamics

◮ Disturbances at frequencies beyond bandwidth of process have to be

removed as much as possible

Lecture 3: Elementary Models for LTI Systems 84 / 91

SLIDE 85

Sampling frequency

Advised to sample process signals at a frequency higher than the sampling frequency essentially needed for identification of process dynamics Choose sampling rate 5-10 higher than bandwidth of process Summarized

1. Sample system at high sampling frequency
2. Preprocess signals
3. Power at frequencies associated to disturbances (high frequencies) will be
low. Filter and decimate (under sample) to bring sampling frequency down

Lecture 3: Elementary Models for LTI Systems 85 / 91

SLIDE 86

User choices

Choice of model structure A(q)yk = B(q) F(q)uk + C(q) D(q)ek Keep in mind that more complex models

◮ . . . are less efficient to compute ◮ . . . are less efficient to use ◮ . . . requires more parameters, therefore exhibits higher variance of θopt

around ˆ θ0, meaning estimates are less accurate Aside from opting for simple model structures, the number of parameters can be penalized using various criterions θopt = argmin

θ

V (θ) + dimθ

N

Lecture 3: Elementary Models for LTI Systems

86 / 91

SLIDE 87

User choices

However, keep in mind that more complex models span larger classes and increase the chance of containing the true model

◮ There is always a trade off in place between flexibility and simplicity ◮ A priori knowledge is valuable ◮ Rule of thumb: first try simple models, and work your way up more

complex models if deemed necessary

◮ Number of parameters can be visually estimated by computing an

experimental frequency response using sinusoidal inputs sequences4

◮ Idea behind information criterions to penalize more complex models is to

asses whether an addition of parameter is capable of improving the fit substantially enough

4this strategy corresponds to sampling in the frequency domain Lecture 3: Elementary Models for LTI Systems 87 / 91

SLIDE 88

User focus

Prediction error methods as optimal (weighted) fitting in the frequency domain V (θ) ≈ 1 2π

π

−π

1 2|G0(ejω) − G(ejω, θ)|2 |L(ejω)|2|U(ω)|2 |H(ejω)|2 dω User can influence the misfit by

◮ Selecting the interesting frequency content of input signals, for example a

frequency bandwidth

◮ Choice of noise model H(q), either empirically or by assumption ◮ Modification of filter L(q); can be interpreted as either a modified input

spectrum L(ejω)U(ω) or a modified noise model L−1(q)H(q)

Lecture 3: Elementary Models for LTI Systems 88 / 91

SLIDE 89

Model validation

How do we validate a model?

1. If model parameters represent physical parameters, assess whether they are

physically meaningful

2. Compare transfer function of estimated models to empirically determined

transfer function

3. Relative error of output signal to simulated output signal

ǫ = 100

N

k=1(yk − y sim k

)2

N

k=1 y 2 k

%

4. Prediction error ǫ can be computed for independent validation datasets; if

the error is much higher compared than the error for the dataset used to estimate the model, it is a strong indication of overmodeling

Lecture 3: Elementary Models for LTI Systems 89 / 91

SLIDE 90

Model validation

5 White noise properties of errors ek can be put to the test, different ‘whiteness tests’ exist in literature 6 Correlation of prediction errors epred

k

to previous input values up to uk. Correlation implies model output has not fully extracted the information present in the input5 Re,u(τ) = 1 N

N

k=τ

epred

k

uk−τ 7 Other mode properties

◮ transfer function ◮ impulse response ◮ poles and zeros ◮ stability ◮ sign (polarity) of the DC gain 5Note this will never be the case for the least squares method Lecture 3: Elementary Models for LTI Systems 90 / 91

SLIDE 91

Model validation

8 Do not forget to inspect the variance of the estimated parameters, too large means too useless! 9 Variance of model parameters, transfer function, . . . is important to determine the working domain of your model such as the frequency range input sequences are allowed to cover. Pole-zero pairs can annihilate one another if their confidence regions overlap significantly, thereby implying a potential model reduction

Lecture 3: Elementary Models for LTI Systems 91 / 91