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–