ZRS

CTU in Prague, CZ - 2009

Model Identification (ARX, (ARX, ARMAX, ARMAX, OE OE model) model) Jirka Roubal [email protected]

hhttp://support.dce.felk.cvut.cz/pub/roubalj/i Department of Control Engineering, Faculty of Electrical Engineering Czech Technical University in Prague

hwww.cvut.czi 2009

J. Roubal, CTU in Prague

–1/19–

ZRS

CTU in Prague, CZ - 2009

Contents and References

Contents 1. Least Squares Method 2. Identification of Dynamical Systems (ARX, ARMAX, OE model) 3. Recursive Identification (forgetting) 4. Examples

References ˇ ˇ TECHA , J. (2000), Modern´ı teorie rˇ´ızen´ı, Praha: Vydavatelstv´ı CVUT. H AVLENA , V. & S ISBN 80-01-02095-9. L JUNG , L. (1999), System Identification: Theory for the User, 2nd edn, Prentice Hall, Upper Saddle River, N.J. ISBN 0-13-656695-2. L UENBERGER , D. G. (1996), Optimization by Vector Space Methods, Wiley-Interscience. ISBN 0-471-18117-X. ´ ı a identifikace system ´ u, N OSKIEVI Cˇ , P. (1999), Modelovan´ ˚ Montanex, a.s., Ostrava. ISBN 80-7225-030-2. J. Roubal, CTU in Prague

–2/19–

ZRS

CTU in Prague, CZ - 2009

Least Squares (LS) Method

1 Least Squares (LS) Method • equation system y = Zθ +e –

y (m × 1) vector of the left sides



Z (m × n), m >> n data matrix exactly known



θ (n × 1) unknown vector



e (m × 1) error – random variable

(1)



E {e} = 0

T

cov {e} = E ee



= σ2 I

• Least Squares (LS) Method (L UENBERGER , D. G. 1996) n o n T o 2 ∗ T θ = arg min e e = arg min y − Z θ y−Z θ = arg min ky − Z θ k2 θ θ θ • minimization =⇒ the best linear nondeviated estimation 

θˆ = Z TZ

J. Roubal, CTU in Prague

−1

Z Ty



P θ = σ 2 Z TZ

−1

(2)

–3/19–

ZRS

CTU in Prague, CZ - 2009

Least Squares (LS) Method

1.1 Example Using the Least Squares method approximate the points in Figure 1 by a linear function. We are looking for parameters k and q of equation

y = kx + q.

3

%% parameters k = 0.2; q = 1; sigma2 = 0.5;

k = theta(1) q = theta(2)

2 y

%% experiments X = [1 2 3 4 5 6]’; y = k*X + q + ... sigma2*randn(size(X)); Z = [X, ones(size(X))]; theta = inv(Z’*Z)*Z’*y

2.5

1.5

1

0.5

points errors line

0 −1

0

1

2

3

4

5

6

7

8

x

Figure 1: Approximation of points by linear function J. Roubal, CTU in Prague

–4/19–

ZRS

CTU in Prague, CZ - 2009

Identification of Dynamical Systems

2 Identification of Dynamical Systems • ARX model (AutoRegresive model with eXternal input) (L JUNG , L. 1999) a(d) y(t) = b(d) u(t) + e(t)

1 b(d) u(t) + e(t) y(t) = a(d) a(d)

– prediction of mean value y ˆ(t|t − 1) is linear function of measurable data – linear regression can be used for model parameters estimation

• ARMAX model (Average model with eXternal input) (L JUNG , L. 1999) a(d) y(t) = b(d) u(t) + c(d)e(t)

b(d) c(d) y(t) = u(t) + e(t) a(d) a(d)

– enables us to model deterministic and stochastic parts of the system independently – linear regression cannot be used for model parameters estimation → pseudolinear reg.

• OE model (Output Error model) (L JUNG , L. 1999) a(d) y(t) = b(d) u(t) + a(d) e(t) J. Roubal, CTU in Prague

b(d) y(t) = u(t) + e(t) a(d) –5/19–

ZRS

CTU in Prague, CZ - 2009

Identification of Dynamical Systems

2.1 Example - Model Identification Using ARX Model

• system identification (L JUNG , L. 1999, N OSKIEVI Cˇ , P. 1999) • consider the following model b3 z 3 + b 2 z 2 + b 1 z + b 0 Y (z) = 4 G(z) = U (z) z + a3 z 3 + a2 z 2 + a1 z + a0

Figure 2: LTI model

J. Roubal, CTU in Prague

–6/19–

ZRS

CTU in Prague, CZ - 2009

Identification of Dynamical Systems

• difference equation y(k) = −a3 y(k − 1) − a2 y(k − 2) − a1 y(k − 3) − a0 y(k − 4)+ + b3 u(k − 1) + b2 u(k − 2) + b1 u(k − 3) + b0 u(k − 4) • n-measuring with noise y(k)

= −a3 y(k − 1) − · · · − a0 y(k − 4) + + b3 u(k − 1) + · · · + b0 u(k − 4) + e(k)

y(k + 1)

= −a3 y(k) − · · · − a0 y(k − 3) + + b3 u(k) + · · · + b0 u(k − 3) + e(k + 1) .. .

y(k + n − 1) = −a3 y(k + n − 2) − · · · − a0 y(k + n − 5) + + b3 u(k + n − 2) + · · · + b0 u(k + n − 5) + e(k + n − 1) • compact form of the previous equation system y = Zθ + e J. Roubal, CTU in Prague

–7/19–

ZRS

CTU in Prague, CZ - 2009

Identification of Dynamical Systems





   y=  



   Z =  



y(k) y(k + 1) .. .

y(k + n − 1)

−y(k − 1) −y(k)

   ,  

··· ···

.. .

−y(k + n − 2)

···



   e=  



e(k)

   ,  

e(k + 1) .. .

e(k + n − 1)

−y(k − 4)

u(k − 1)

−y(k − 3)

u(k)

.. .

.. .

−y(k + n − 5) u(k + n − 2)

      θ=       ··· ···

a3 .. .

a0 b3 .. .

b0



      ,      

u(k − 4) u(k − 3) .. .

···

u(k + n − 5)

      

• measured data are in Figure 3 and the result of the identification experiment is in Figure 4

J. Roubal, CTU in Prague

–8/19–

CTU in Prague, CZ - 2009

Identification of Dynamical Systems

1

u

0.5 0 −0.5 −1 0

5

10

15

10

15

t [s] 1

0.5 y

ZRS

0

−0.5 0

5 t [s]

Figure 3: Data for identification – system input and system output

J. Roubal, CTU in Prague

–9/19–

CTU in Prague, CZ - 2009

Identification of Dynamical Systems

1.8 System Model (n=15) Model (n=30) Model (n=100)

1.6 1.4 1.2 1 y

ZRS

0.8 0.6 0.4 0.2 0 −0.2 0

0.5

1

1.5

t [s]

Figure 4: Result of the identification experiment – step response

J. Roubal, CTU in Prague

–10/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

3 Recursive Identification • one new measured data y(t) = z T(t) θ + e(t)

e ∼ N (0, σe2 )

we know

E {θθ (t)|t − 1} = θˆ θ(t|t − 1)

cov {θθ (t)|t − 1} = σe2 P (t|t − 1)

• recursive estimation of constant parameters (H AVLENA , V. & Sˇ TECHA , J. 2000)   P (t|t − 1)z(t) T θˆ θ(t|t) = θˆ θ(t|t − 1) + y(t) − z (t) θˆ (3) θ(t|t − 1) T 1 + z (t)P (t|t − 1)z(t) P (t|t − 1)z(t)z T (t)P (t|t − 1) P (t|t) = P (t|t − 1) − 1 + z T (t)P (t|t − 1)z(t)

J. Roubal, CTU in Prague

(4)

–11/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

– equations (3) and (4) do not develop in time

θˆ θ(t + 1|t) = θˆ θ(t|t)

P (t + 1|t) = P (t|t)

because the constant parameters are supposed – recursive identification with linear forgetting (H AVLENA , V. & Sˇ TECHA , J. 2000)

P (t + 1|t) = P (t|t) + P 0

(5)

where P 0 can include a prior information about the direction of the parameters variance – recursive identification with exponential forgetting (H AVLENA , V. & Sˇ TECHA , J. 2000)

P (t + 1|t) = where ϕ

1 P (t|t) ϕ

(6)

∈ (0, 1i is forgetting coefficient

J. Roubal, CTU in Prague

–12/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

3.1 Example in Matlab Experiments:

1.125 G(z) = 2 z + 0.25z − 0.125 (sampling time is Ts

t=0.5 s



1.16 G(z) = 2 z + 0.16

= 0.01 s)

• σe = 0.01 – without forgetting –

P0 = 0



ϕ = 0.95

P 0 = 0, ϕ = 1 influence of σe to estimation

P 0 = 0.01I ϕ = 0.9

P 0 = 0.1I

P0 = I

ϕ = 0.85

• σe = 0.1

J. Roubal, CTU in Prague

–13/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

2

P0 = 0.1, φ = 0.85, σe = 0.01 1

1

u

u

0.5

0 −1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

estimated exp 0.1 0.2 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

−1 0

θ1

−0.5

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.2 0.1 0 −0.1 0

real estimated estimatedlin

0.2 θ2

1

0 −0.2

0

0

−0.5

1

0

5

−1

θ

y

0.5

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.5 0 0

Figure 5: Parameters identification σe

J. Roubal, CTU in Prague

= 0.01

–14/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

2

P0 = 0.1, φ = 0.85, σe = 0.1 1

1 u

0.5

−1 0

0

u

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

estimated exp 0.1 0.2 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

−1 0

θ1

−0.5

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.2 0.1 0 −0.1 0

real estimated estimatedlin

0.2 θ2

1

0

0 1 5

−1

θ

y

−0.2 0

0

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.5 0 0

Figure 6: Parameters identification σe

J. Roubal, CTU in Prague

= 0.1

–15/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

2

P0 = 0.1, φ = 0.85, σe = 0.01 1

1 u

0.5

−1 0

0

u

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

estimated exp 0.1 0.2 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

−1 0

θ1

−0.5

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.2 0.1 0 −0.1 0

real estimated estimatedlin

0.2 θ2

1

0

0 1 5

−1

θ

y

−0.2 0

0

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.5 0 0

Figure 7: Parameters identification σe

J. Roubal, CTU in Prague

= 0.01

–16/19–

ZRS

CTU in Prague, CZ - 2009

Recursive Identification

2

P0 = 0.1, φ = 0.85, σe = 0.1 1

1 u

0.5

−1 0

0

u

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

estimated exp 0.1 0.2 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

−1 0

θ1

−0.5

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.2 0.1 0 −0.1 0

real estimated estimatedlin

0.2 θ2

1

0

0 1 5

−1

θ

y

−0.2 0

0

0.1

0.2

0.3

0.4

0.5 t [s]

0.6

0.7

0.8

0.9

1

0.5 0 0

Figure 8: Parameters identification σe

J. Roubal, CTU in Prague

= 0.1

–17/19–

ZRS

CTU in Prague, CZ - 2009

Problems

Problems • Convergency of the Least Squares method. ˇ TECHA , J. 2000) – sufficient excitation condition (L JUNG , L. 1999, H AVLENA , V. & S

• In (1), the unknown error is also in the left side. – TLS Total Least Squares methods

• What type of forgetting to use? • How to set the foregetting parameters? • For the system identification the system must be excited which is an opposite to control that try to stabilize the system.

J. Roubal, CTU in Prague

–18/19–

ZRS

CTU in Prague, CZ - 2009

Questions ...

Questions ... • How the Least Squares method is defined? • What is the condition for LS solution (2)? • What is the convergency condition for dynamical system identification? • Consider that we measured the distances (in the plane) of the stars in the ring of the Orion constellation (see Figure 9). Using the Least Squares method compute the distances a ˆ b from measured data: and ˆ – a) a

= 1.1, b = 1.2, c

– b) a

= 1.1, b = 1.2, c = 2.3,

– c) a

= 1.1, b = 1.2, c = 2.2.

b a

Figure 9: Orion constellation J. Roubal, CTU in Prague

–19/19–