[PPT] - Rational bases for system identification Adhemar Bultheel PowerPoint Presentation

SLIDE 1

Rational bases for system identification

Adhemar Bultheel

Department Computer Science Numerical Approximation and Linear Algebra Group (NALAG) K.U.Leuven, Belgium adhemar.bultheel@cs.kuleuven.ac.be http://www.cs.kuleuven.ac.be/∼ade/

September 2002

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

SLIDE 2

1/63

Summary

1. Definitions from signals, systems
2. Vector orthogonal polynomials (Jacobi/unitary Hessenberg matrices)
3. Orthogonal rational functions in L2(R, µ), L2(T, µ)
4. ORF and inverse (Hessenberg) eigenvalue problem (semiseparable

matrices)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 3

2/63

Definitions

A signal is a complex function of a time variable.
discrete (fn), n ∈ Z (discrete signal) or
continuous f(t), t ∈ R (analog signal)
Define a delay operator (Df)n = fn−1

(Df(t) = f(t − 1))

A pulse δ is the Kronecker delta

(Dirac delta function)

A pulse decomposition:

f = (

m fmDm)δ

(f =

f(s)Dsδ(·)ds)
An inner product: f, g = fngn

(f, g =

f(t)g(t)dt)
The convolution h = f ∗ g:

hn =

m fmgn−m

(h(t) =

f(s)g(s − t)ds)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 4

3/63

The (auto)correlation function:

(rfg)n =

m fmgn+m

(rfg(t) =

f(s)g(s + t)ds)

(rff)n = rn =

m fmfn+m

(rff(t) = r(t) =

f(s)f(s + t)ds)
The Fourier transform: F = Ff

F(eiω) =

m fme−imω

(F(ω) =

f(t)e−itωdt)
The Z-transform:

F(z) =

m fmz−m

(F(z) =

f(t)e−itzdt)

(compare with (2-sided) Laplace transform F(p) =

f(t)e−tpdt)

F(z)|z=eiω∈T = (Ff)(eiω) F(z)|z=ω∈R = (Ff)(ω)

power spectrum: R = Fr = |F|2

white noise: R = constant; normalized constant = 1 colored noise: g = Hf, f white noise.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 5

4/63

A (linear) system H is a (linear) operator that maps a signal (input) to

another signal (output) Hu = y

It is stationary or time invariant if HD = DH
The impulse response: h = Hδ

y = Hu = H ( umDm) δ = umDm(Hδ) = ( umhn−m) = h ∗ u y = Hu = H

u(s)Ds

δ =

u(s)Ds(Hδ) =
u(s)h(· − s) = h ∗ u

y = h ∗ u = u ∗ h = hmDmu = H(D−1)u hence H = H(D−1).

The transfer function: H(z) = (Zh)(z) = |H(z)|eiΦ(z)

amplitude: |H(z)| and phase: Φ(z).

A causal system: if u = 0 before time 0, then y = 0 before time 0.
A (BIBO) stable system: if u < ∞ ⇒ Hu < ∞

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 6

5/63

For a causal stable system: H ∈ H2(E)

(H ∈ H2(L))

The inverse system: if H(z) = 0, z ∈ T (z ∈ R), inverse = 1/H(z).

Stable & causal then inverse if poles of H in D (in U) Inverse is table & causal if zeros of H in D (in U) the latter is called minimal phase

A system is all pass if |H(z)| = 1 (hence Y = HU = U)
FIR/MA: H is polynomial in z−1 (otherwise IIR)

AR: H = cnst/polynomial in z−1 ARMA: H = rational in z−1 proper: deg(numerator) ≤ deg(denominator) (⇐ stable & causal)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 7

6/63

State Space

If H is proper rational then partial fraction expansion: H(z) = D + c1 z − α1 + · · · + cn z − αn = C(zI − A)−1B + D where A = diag(α1, . . . , αn), B = [1, . . . , 1]T, C = [c1, . . . , cn], D = H(∞). This is called a state space description of the system. It is characterized by the 4 matrices (A, B, C, D). Such a state space description is not unique.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 8

7/63

If T is an invertible matrix, then (A, B, C, D) ≡ (A′, B′, C′, D′) if

A′

B′ C′ D′

=
T

1 A B C D T 1 −1 =

TAT −1

TB CT −1 D

since C(zI − A)−1B + D = C′(zI − A′)−1B′ + D′.

A minimal realization: size A is minimal. Then σ(A) = poles of the transfer function. H(z) = hk

zk = D + C(zI − A)−1B = D + CB z + · · · + CAk−1B zk

+ · · ·. Y (z) = [C(zI − A)−1B + D]U(z) = CX(z) + DU(z) with X(z) = (zI − A)−1BU(z) or zX(z) = AX(z) + BU(z). This zX(z) = AX(z) + BU(z) Y (z) = CX(z) + DU(z)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 9

8/63

is the z-domain formulation of xn+1 = Axn + Bun yn = Cxn + Dun ˙ x(t) = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t) ˙ x(t) = 1 i dx(t) dt if system is causal and stable and state is zero at time zero. x is the state. SS description very powerful because of

linear algebra techniques applicable
generalization of SISO to MIMO is relatively simple

n nu n A B ny C D

many other advantages

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 10

9/63

Identification problem

U U N U G N Y Y0 Y

Y = G(U − NU) + NY = GU + (NY − GNU) = GU + V , V = H0E, E white noise, H0 some coloring filter. It is assumed we know Y and U (or G) in a number of frequency points: zk ∈ T or R and some information about the noise (covariances)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 11

10/63

Frequency Domain Identification

Estimate G as ˆ G such that we minimize Y − ˆ Y wy or G − ˆ Gwg where F2

w = N k=1 w2 k|F(zk)|2.

Note with wg = wyU: Y − ˆ Y wy = GU − ˆ GUwy = (G − ˆ G)Uwy = G − ˆ Gwg, E.g. measurements {G(zj)}N

j=1, variance {σj}N j=1, then set wg = 1/σ.

G − ˆ Gwg =

Y

U − B A

wg

=

Y A − BU

AU

wg

= Y A − BUw, w = wg AU In ML estimation: w−2 = Cov(AY − BU) = σ2

Y |A|2 + σ2 U|B|2 − 2Re(σ2 Y UAB).

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 12

11/63

Linear/Nonlinear problem

Suppose ˆ G parametrised by parameter set θ. Cost function: C(θ) = C(L(θ), N(θ)), L linear, N nonlinear. If N does not depend on θ, the problem becomes a linear LSQ problem. Solving for min C(θ) directly is costly and not convex so convergence to global min is not guaranteed, unless we use a good starting point. The method of solution: ITERATE θ(k+1) = argmin C(L(θ), N(θ(k))). May work in some cases to find global solution. Requires a stable and efficient method to solve the linear LSQ problem. Need constraints. E.g. system stability needs system poles in D. So have to solve a constrained and weightes NLSQ problem in C. ITERATE converges, then use as starting point for min C(θ) with standard routine.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 13

12/63

Polynomial method

min C = min Y A − UB2

w w may depend upon (estimated) solution.

Set W = w[Y − U], P = [A B]T ⇒ w(Y A − UB) = WP C = WP2 = P2

W = N

j=1

P(zj)HW H

j WjP(zj),

P ∈ P2×1

n

. f, gW = f(zj)HW H

j Wjg(zj) a p.d. inner product in P2×1 n

if n ≤ N. min PW needs constraint: e.g. P “monic”. How to parametrize P?

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 14

13/63

dim P2×1 = 2(n + 1):

1 z−1 · · · z−n · · · · · · 1 z−1 · · · z−n

r

1 z−1 · · · z−n 1 z−1 · · · z−n

Simpler to orthogonalize {I, z−1I, . . . , z−nI} w.r.t. matrix-valued SP

f, gW = f(zj)HW H

j Wjg(zj) in P2×2

Gives orthogonal block polynomials φk: φk, φlW = δk−lI to make unique choose e.g. l.c. upper triangular. Write P = n

k=0 φkλk = Φ(z)λ, φk ∈ P2×2 k

, λk ∈ C2×1, λn = 0. P2

W = ΦλW = λHΦHWHWΦλ, Φ = [φj(zi)], W = diag(Wk).

φk orthonormal ⇒ ΦHWHWΦ = QHQ = I, Q = WΦ is unitary min P2

W = min λ2 = λH n λn ⇒ P = φnλn.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 15

14/63

Recurrence relation

z−1φk−1(z) = φk(z)ηkk + · · · + φ0(z)η0k, k = 1, . . . , n + 1 z−1[φ0 · · · φn−1 φn] = [φ0 · · · φnφn+1]       η01 · · · η0n η0,n+1 η11 · · · η1n η1,n+1 ... . . . . . . ηnn ηn,n+1 ηn+1,n+1       z−1Φ(z) = Φ(z)H + [0 · · · 0φn+1ηn+1,n+1] Now write this down for z = [z1, . . . , zN]T. ˜ Z−1Φ = ΦH + [0 · · · 0 φn+1(z)ηn+1,n+1] ˜ Z = diag(zI2), Φ = [φj(zi)] = Φ(z)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 16

15/63

˜ Z−1Φ = ΦH + [0 · · · 0 φn+1(z)ηn+1,n+1] multiply with W and set Q = WΦ (W˜ Z−1 = Z−1W) Z−1Q = QH + [0 · · · 0 Wφn+1(z)ηn+1,n+1] multiply from the left with QH and because φn+1 ⊥W Pn QHZ−1Q = H. For N = n + 1, Q is square unitary matrix and thus QHZ−1Q = H ⇔ Z−1 = QHQH

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 17

16/63

In the analog case: Z is real ⇒ H = QHZ−1Q = HH. ⇒ H block tridiagonal ⇒ scalar pentadiagonal (block 3-term recurrence relation) In the digital case: Z−HZ−1 = I ⇒ HHH = QHZ−1QQHZ−HQ = I. ⇒ H is unitary (block) upper Hessenberg. (Szeg˝

-type recurrence relation)

Z−1 = QHQH is Inverse eigenvalue problem. Not completely defined. Need some condition on the W Set w =   W1 . . . WN   ⇒ φ0, φ0W = φH W H

k Wkφ0 = φH 0 wHwφ0 = I2.

With η00 = φ−1 = sqrt(wHw): φ0, I = η00. So QHw = ΦHWHw = ΦHWHW[I, . . . , I]H = [ηH

00, 0, . . . , 0]H ⇒

QHw = e1η00

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 18

17/63

Inverse Eigenvalue Problem

QHZ−1Q = H & QHw = e1η00 & ηH

00η00 = wHw

r

QH   W1 . . . Z−1 WN   I2 Q

=

  η00 H  

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 19

18/63

Hessenberg matrix

    w1 z−1

1

w2 z−1

2

... wN z−1

N

    unitary similarity ⇒ transf. Block Upper Hessenberg × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × | × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 20

19/63

Recursively: continuous time

× × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ⇒ × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 21

20/63

Recursively: discrete time

Szeg˝

-type recurrence:

φnσn = z−1φn−1 + φ∗

n−1γn

φ∗

nσH n

= z−1φn−1γH

n + φ∗ n−1

σn, γn ∈ C2×2 (block Schur parameters) Store the block upper Hessenberg matrix in factored form H = G1G2 · · · Gm Gk =     I2(k−1) −γk ˆ σk σk ˆ γk I...     Chasing the elements to block upper Hessenberg form by similarity transformations percolating through the product requires operations on 3 × 3 or 5 × 5 blocks ⇒ fast algorithm.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 22

21/63

Remarks

Fast algorithm (O(N 2) operations)
Economize if n < N (operate on N × n matrix)
can do in real arithmetic for complex data
Works for MIMO systems (OVP are s × 1, s = (nu + ny)ny)
Several degree structures are possible (column permutations)
Representation of ˆ

G is optimal: representation = orthogonal basis If the weights change then basis changes and is optimal. The condition number of systems to solve = 1.

nonlinear (ML) solution:

fix last basis, condition number grows moderately

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 23

22/63

real computations for T

We can preprocess complex data to do real arithmetic. 2 data points eiωk, e−iωk ∈ T with corresponding weights Wk and W k. Define Q =

1 √ 2

1

1 1 −1 1 −i

Then QH

Wk W k

=

1 √ 2

1

i 1 1 1 −1 Wk W k

=

√ 2

Re(Wk)

Im(Wk)

and QH
eiωk

e−iωk

Q =
cos ωk

sin ωk − sin ωk cos ωk

M´

etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 24

23/63

MIMO system

M= stack columns in long vector

model A−1B: A is ny × ny, B is ny × nu

P =
A
B
this is ny(ny + nu) long

E = AY − BU (ny × 1) ⇒ min

k

EH(zk)C−1

E (zk)

E(zk)

E = ([Y T

− U T] ⊗ Iny) P C

E = Cov(

E) = C1+C2−2C3    C1 = (Inu ⊗ A)C

Y (Inu ⊗ AH)

C2 = (Inu ⊗ B)C

U(Inu ⊗ BH)

C3 = Herm[(Inu ⊗ A)C

Y U(Inu ⊗ BH)]

min PW, W = ([Y T − U T] ⊗ Iny)C−1/2

E

has size 1 × (ny + nu)ny

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 25

24/63

Robot arm

rder [6/6], 100 data points, 200Hz, discrete time

0.05 0.1 −60 −50 −40 −30 −20 −10 10 20 amplitude Robot arm (dB) 0.02 0.04 0.06 0.08 0.1 0.12 −50 −40 −30 −20 −10 10 rel mag Cmplx Err (dB) 0.02 0.04 0.06 0.08 0.1 0.12 0.1 0.2 0.3 0.4 0.5 0.6 0.7 mag Cmplx Err Vml=109.943

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 26

25/63

Band pass filter

rder [6/6], 50 data points, 20KHz, discrete time

0.05 0.1 0.15 0.2 −80 −70 −60 −50 −40 −30 −20 −10 amplitude frf (dB) 0.05 0.1 0.15 0.2 0.25 −80 −60 −40 −20 rel mag Cmplx Err (dB) 0.05 0.1 0.15 0.2 0.25 1 2 3 4 5 x 10−4 mag Cmplx Err Vml=18701.9636

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 27

26/63

CD’s radial servo system

sampling frequency 9.7 kHz, order [5/5]

Relative Freq 0.1 0.2 0.3

60
40
20

20 40 Magnitude FRF (dB) error frf

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 28

27/63

Mechanical structure

rder [120/120]

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 29

28/63

References OPV

[1] A. Bultheel, M. Van Barel, Y. Rolain, J. Schoukens. Robust rational approximation for identification, Proc. CDC conference, 2001 [2] R. Pintelon, Y. Rolain, A. Bultheel, M. Van Barel. Numerically robust frequency domain identification of multivariable systems, (2002). [3] M. Van Barel, A. Bultheel, Discrete linearized least squares approximation on the unit circle, J. Comput. Appl. Math. 50 (1994) 965-972. [4] A. Bultheel, M. Van Barel. Vector orthogonal polynomials and least squares approximation, SIAM J. Matrix Anal. Appl. 16 (1995) 863-885. [5] M. Van Barel, A. Bultheel, Orthogonal polynomial vectors and least squares approximation for a discrete inner product, Electron. Trans. Numer. Anal. 3 (1995) 1-23.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 30

29/63

Identification and ORF

Identification problem: find rational function ˆ G such that min

N

k=1

|wk(G(zk) − ˆ G(zk))|2 The wk can possibly depend on the solution, but in the present situation we use an estimate of the solution so that we can consider wk as known. Looking for a rational ˆ G: ˆ G ∈ Ln =

pn(z)

n

j=1(1 − αj z ) : pk ∈ Pn

Note that if all αk = 0 ⇒ Ln = Pn.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 31

30/63

So why not represent it w.r.t. a basis of rational functions: Ln = span{B0, . . . , Bn} with B0 = 1, Bk =

k

j=1

1 z − αj , k ≥ 1 For numerical stability reasons: use orthogonalize the basis {Bk} → {φk} ˆ G(z) =

n

k=0

λkφk(z), φk orthogonal basis functions with respect to an appropriate inner product. The inner product can be discrete or continuous, but is in general with respect to a weight (measure) on T or R. f, gµ =

f(z)g(z)dµ(z)
r

f, gw =

N

j=1

f(zj)g(zj)wj.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 32

31/63

Takenaka-Malmquist basis

Work by Ninnes, Van den Hof, Heuberger, Bokor, de Hoog, and others: Use ORF with respect to Lebesgue measure on T. These have explicit expressions φn(z) =

1 − |αn|2

z − αn

n−1

k=1

1 − αkz z − αk but that does not help for the condition number of ΦHWHWΦ in the linear least squares problem. And/or use a finite number of αk that are cyclically repeated.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 33

32/63

State space formulation

parahermitian conjugate f∗(z) = [f(1/z)]H or f∗(z) = [f(z)]H
All-pass/inner function: G(z)G∗(z) = I/ analytic
balanced state space: G(z) = C(zI − A)−1B + D

Xc = ∞

k=0(AkB)(AkB)H = Σ = Xo = ∞ k=0(CAk)H(CAk)

all-pass+balanced ⇒ Σ = I and
A

B C D A B C D H = I =

A

B C D H A B C D

matrix valued inner product: F, G =

1 2π

F(eiω)[G(eiω)]Hdω

Xc =

(zI − A)−1B, (zI − A)−1B
,

Xo =

[C(zI − A)−1]H, [C(zI − A)−1]H

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 34

33/63

Thus: take balanced SS realization (An, Bn, Cn, Dn) of Blaschke product
f degree n, and set

(zI − An)−1Bn = [φ1(z), φ2(z), . . . , φn(z)]T and the φk are ORF.

Product of Blaschke factors ⇒

recursive SS (Ak, Bk, Ck, Dk) ⇒ Takenaka-Malmquist ORF. SS for (1 − αz)/(z − α) is

A

B C D

=
α
1 − |α|2
1 − |α|2

−α

Hambo basis: inner fct = (Blaschke prod)k ⇒ cyclic repetition of finite

number of poles.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 35

34/63

Special case (1)

Laguerre basis: all αk = α: φn(x) = √

1−|α|2 z−α

1−αz

z−α

n−1 . Laguerre polynomials Ln(x): ∞

0 e−xLn(x)Lm(x)dx = δn−m.

Define ˆ φn(x) = e−x/2Ln(x) Take Fourier transform: cφn(x) =

√ 2α jx+α

jx−α

jx+α

n . This is the continuous time Laguerre basis (orthogonal on R). A bilinear (Jukowky) transform to T gives the discrete time Laguerre basis.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 36

35/63

Special case (2)

Kautz basis: pairs of complex conjugate poles. φ2k−1(z) = Ck(1 − βkz)ψk(z), φ2k(z) = C′

k(1 − β′ kz)ψk(z),

ψk(z) =

k−1

j=1(1−αjz)(1−αjz)

k

j=1(z−αj)(z−αj)

Obtained from Legendre polynomials: 1

−1 Pn(x)Pm(x)dx = δn−m.

Define ˆ φ(t) = e−γtPn(2e−γt − 1) ⇒⊥ on [0, ∞) Taking Fourier transform: cnφn(x) =

cn jx+(n+1/2)γ

n−1

k=0 jx−(k+1/2)γ jx+(k+1/2)γ.

This is the continuous time Kautz basis (orthogonal on R). Set αk = 2−γ(2k+1)

2+γ(2k+1) and bilinear transform to T gives discrete time Kautz

system for a specific choice of the αk.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 37

36/63

Recursive state space

Φn(z)T = (zI − An)−1Bn with (An, Bn, Cn, Dn) a balanced state space realization of a general Blaschke product n

k=1(1 − ¯

αkz)/(z − αk) If Gi ↔ (Ai, Bi, Ci, Di), i = 1, 2 (inner ↔ balanced) then G = G1G2 ↔ (A, B, C, D) with A B C D

=

  A1 B1 B2C1 A2 B2D1 D2C1 C2 D2D1   Thus with G1 = n−1

k=1 1−¯ αkz z−αk ↔ (An−1, Bn−1, Cn−1, Dn−1) and

G2 = 1−¯

αnz z−αn ↔ (αn,

1 − |αn|2,
1 − |αn|2, −¯

αn) An =

An−1

Xn αn

,

Bn =

Bn−1

bn

M´

etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 38

37/63

An =

An−1

Xn αn

,

Bn =

Bn−1

bn

Xn =
1 − |αn|2Cn−1,

bT

n =

1 − |αn|2Dn−1,

Cn = [(−¯ αn)Cn−1 |

1 − |αn|2],

Dn = (−¯ αn)Dn−1.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 39

38/63

Hambo basis

Take G(z) = d

j=1 1−¯ αjz z−αj

with balanced state space realization (A, B, C, D). Define Vk(z) = (zI − A)−1BGk−1(z), set Vk(z) = [φk1, . . . , φkd]T. Then the {φk,i : i = 1, . . . , d; k = 0, 1, . . .} form an orthonormal set ORF for f, g =

1 2π

2π f(eiω)g(eiω)dω. Note that vk(j) = Vk(D)δj is a vector in the state space. Orthogonal state space basis in the sense that [Vk, Vk] = I ⇔ {vk, vl} = I with [f, g] =

1 2π

2π f(eiω)[g(eiω)]Hdω and {f, g} =

j f(j)[g(j)]H.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 40

39/63

Identification

Make row vector of ORF: Φ(z) = [(zI − A)−1B]T. Model: ˆ G = Φλ, set Φ = Φ(z), G = G(z), W = diag(w), E = G − Φλ. minimize G− ˆ Gw =

k |wk[G(zk)−Φ(zk)λ]|2 = EHWHWE = E2 w.

Hence solve Φλ = G in weighted least squares sense: λ = Φ†G, Φ† = [ΦHWHWΦ]−1[ΦHWHW] However δk−l = φk, φl =

1 2π

π

−π φk(eiω)φl(eiω)dω = j |wj|2φk(zj)φl(zj).

Hence ΦHWHWΦ = I. Condition number of WΦ can be very large. Loss of correct digits in λ (! iterate on the nonlinear parameters α)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 41

40/63

ORF recurrence

OPS: z−1φk(z) ∈ span{φl : l = k + 1, k, . . . , 0} ORF:(generic case = regular system: φk(αk−1) = 0) z−1φk(z) ∈ span{[1 − αl/z]φl(z) : l = k + 1, . . . , 0} (α0 = 0) φk(z) ∈ span{[z − αl]φl(z) : l = k + 1, . . . , 0} (α0 = 0) Φ(z) = Φ(zI − A)H + [0 . . . 0 (z − αn+1)φn+1bn] Φ = [φ0, . . . , φn], A = diag(α0, α1, . . . , αn), H =     η01 · · · η0n η0,n+1 η11 · · · η1n η1,n+1 ... . . . . . . ηnn ηn,n+1     n = N: φN+1 does not exist (φN+1w = 0): choose ΦN+1(zk) = 0, k = 1, . . . , N.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 42

41/63

ORF recurrence (continuous time)

In case of real line: 3-term recurrence for OPS. Similar for ORF (if αk ∈ R). φk(z) ∈ span{[z − αl]φl(z) : l = k − 1, k, k + 1}, (φ−1 = 0) φn(z) = bn−1(z − αn−1)φn−1 + an(z − αn)φn + bn(z − αn+1)φn+1. Φ = Φ(zI − A)T + [0...0 (z − αn+1)φn+1bn] Φ = [φ0, . . . , φn], A = diag(α0, α1, . . . , αn), T =       a0 b0 b0 a1 b1 b1 ... ... ... ... bn−1 bn−1 an       If αk are complex, T not tridiagonal but . . . (see later).

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 43

42/63

ORF recurrence (discrete time)

Forward recursion

φn(z)

φ∗

n(z)

= en

z − αn−1 z − αn

1

Ln Ln 1 ζn−1(z) 1 φn−1(z) φ∗

n−1(z)

Ln = −
φk, 1−αn−1z

z−αn φn−1

φk, z−αn−1

z−αn φ∗ n−1

, en = 1 − |αn|2 1 − |αn−1|2 1 1 − |Ln|2 1/2 ζk(z) = 1 − αkz z − αk , Bn = ζ1 · · · ζn, φ∗

n(z) = Bn(z)φn(1/z).

Backward recursion = Nevanlinna-Pick algorithm

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 44

43/63

Recurrence (discrete time)

Eliminate the φ∗

k ⇒ φk(z) = k+1

j=0

ηjk(z − αj)φj(z) Φ(z) = Φ(z)(zI − A)H + [0 . . . 0 (z − αn+1)φn+1(z)bn] H is not unitary (see later)

Conclusion

continuous time: H almost tridiagonal. discrete time: H almost unitary Hessenberg

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 45

44/63

Recursive state space

Φn(z)T = (zI − An)−1Bn with (An, Bn, Cn, Dn) a state space realization for the ORF An =

An−1

Xn αn

,

Bn =

Bn−1

bn

If Xn = [xn1, . . . , xnn], then

ρnn = 1, ρn,n−k−1 = en−k(−¯ αn−k−1)(ρn,n−k + Ln−kσn,n−k) σnn = 0, σn,n−k−1 = en−k(¯ Ln−kρn,n−k + σn,n−k), xn,n−k = en−k(1 − |αn−k−1|2)(ρn,n−k + Ln−kσn,n−k) and bn = e1(L1ρn,1 + σn,1) 1, 1−1/2

w

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 46

45/63

Inner product evaluation

If data are available in z = {zj}N

j=1 ⊂ T with corresponding weights

w = {wj}N

j=1 then choose the discrete inner product.

For a continuous inner product, choose points on the circle, e.g. equidistant zj = exp{2ijπ/N}, j = 1, . . . , N and evaluate weight wj = w(zj) and use discrete inner product. For example w(z) = |U(z)|2, and zj equidistant on T, then given time domain data uk, this can be computed very efficiently by FFT. We can evaluate the moments cl in |U(z)|2 =

l∈Z clzl by convolution

and use a Nevanlinna-Pick type algorithm on the PR function Ω(z) = c0/2 + ∞

l=1 clzl. (approximately)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 47

46/63

The linear least squares problem

So we can compute Φ(α) = Φ(z; α) = [φ0(z) | φ1(z) | · · · | φn(z)] ∈ CN×n+1 W = diag(w), λ = {λk}n

k=0,

G = G(z) Setting ˆ Gn = n

k=0 λkφk, solve in least squares sense

WΦ(α)λ(α) = WG ≡ Q(α)λ(α) = WG, Q(α) = WΦ(α) by orthogonality however, λ(α) = [Q(α)HQ(α)]−1Q(α)HWG = Q(α)HWG

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 48

47/63

Nonlinear problem

K(α) =

j

|G(zj) −

n

k=0

λk(α)φk(zj; α)|2wj ⇒ min

α⊂D K(α)

= (I − Φ(α)Q(α)HW)G2

W = (I − Q(α)Q(α)H)WG2

constraint: set αj = rjeiωj with −1 ≤ rj ≤ 1, j = 1, . . . , n. αj can be forced to be real or complex conjugate. Therefore use parameters bi, ci where (z − αi)(z − αi) = z2 + biz + ci with appropriate constraints.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 49

48/63

Robot arm

100 data points, 200 Hz, condensed in the beginning, order [6/6] discrete time

0.2 0.4 0.6 0.8 1 10

−4

10

−2

10 10

2

ampl frequency response fct 1645 0.2 0.4 0.6 0.8 1 −4 −2 2 4 phase frequency response fct 1645 0.2 0.4 0.6 0.8 1 10

−6

10

−4

10

−2

10 10

2

variance 0.2 0.4 0.6 0.8 1 10

−2

10

−1

10 10

1

relative error 1645

−1 −0.5 0.5 1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 poles and zeros 1645

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 50

49/63

Bandpass filter

rder [6/6], 50 data points, 20KHz, discrete time

0.5 1 1.5 2 10

−6

10

−4

10

−2

10 10

2

ampl frequency response fct 623 0.5 1 1.5 2 −3 −2 −1 1 2 3 4 phase frequency response fct 623 0.5 1 1.5 2 10

−12

10

−11

10

−10

10

−9

variance 0.5 1 1.5 2 10

−6

10

−4

10

−2

10 relative error 623

−1 −0.5 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 poles and zeros 637

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 51

50/63

Electrical Machine

110 data points, 4000 Hz, condensed in the beginning, order [5/5], discrete time

0.1 0.2 0.3 0.4 0.5 10

−2

10

−1

10 ampl frequency response fct 594 0.1 0.2 0.3 0.4 0.5 −0.5 0.5 1 1.5 phase frequency response fct 594 0.1 0.2 0.3 0.4 0.5 10

−11

10

−10

10

−9

10

−8

10

−7

variance 0.1 0.2 0.3 0.4 0.5 10

−4

10

−3

10

−2

10

−1

relative error 594

−1 −0.5 0.5 1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 poles and zeros 594

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 52

51/63

Sensitivity, Electrical Machine

−1 −0.5 0.5 1 10

−3

10

−2

10

−1

10 relative error, variable parameters 5

poles: [-0.3999

0.3989

0.9437 0.9935 0.9971]

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 53

52/63

Conclusion

Rational approximation in weighted discrete least squares sense.
Linear and nonlinear parameters. Solve for linear in terms of nonlinear.

Use orthogonal basis to minimize the condition number.

Either direct orthogonal rational basis or linearized vector polynomial as

a combination of orthogonal block polynomials.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 54

53/63

References ORF

[1] A. Bultheel, P. Gonz´ alez-Vera, E. Hendriksen, and O. Nj˚

astad. Orthogonal rational
functions. Cambridge University Press, 1999.

[2] Th. de Hoog. Rational orthonormal bases and related transforms in linear system modeling PhD, T.U. Delft, 2001 [3] P. Van gucht, A. Bultheel. Using orthogonal rational functions for system identification, Report TW314, Dept. Computer Science, K.U.Leuven, September 2000. [4] P. Van gucht and A. Bultheel. State-space realization and orthogonal rational

functions. Report TW292, Dept. Computer Science, K.U.Leuven, September 1999.

[5] P. Van gucht, A. Bultheel. Matlab routines for system identification using orthogonal rational functions. http://www.cs.kuleuven.ac.be/˜nalag/research/software/ORF/ORFidentification.html.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 55

54/63

ORF as inverse eigenvalue problem

Φ(z) = Φ(z)(zI − A)H + [0 . . . 0 (z − αn+1)φn+1bn] write for all zj: (recall φN+1(z) = 0): Φ = (ZΦ − ΦA)H multiply with W with Q = WΦ: Q = (ZQ − QA)H assuming S = H−1 exists: Q(S + A) = ZQ multiply with QH: S + A = QHZQ as before: QHw = e1η00, w = [W1, . . . , WN]T, η00 = w QH   w Z     1 Q   =   w S + A  

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 56

55/63

Semiseparable matrices

A semiseparable matrix is of the form S = tril(uvH, 0) + R, e.g. S =     u1¯ v1 r12 r13 r14 u2¯ v1 u2¯ v2 r23 r24 u3¯ v1 u3¯ v2 u3¯ v3 r34 u4¯ v1 u4¯ v2 u4¯ v3 u4¯ v4     its inverse is upper Hessenberg H. In our case S = H−1 and from the solution of the inverse eigenvalue problem it follows that for zi all on R or all on T, R is also “rank 1”: R = triu(stH, 1). Hence fast algorithms exist of O(n2).

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 57

56/63

Special cases

continuous time case: zk ∈ R S + A = QHZQ, ⇒ S + A = (S + A)H SH − S = A − AH = 2iIm(A), α’s are real or come in pairs: (a + ib, −a + ib). Hence the off-diagonal part of S is hermitian ⇒ R = triu(vuH, 1) S = tril(uvH, 0) + triu(vuH, 1) and Im(diag(S)) = −Im(A). Such an H = S−1 = tridiagonal (only tridiagonal when αk’s real)

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 58

57/63

discrete time case: zk ∈ T QHZQ unitary ⇒ S + A unitary. α’s are real or come in complex conjugate pairs: (a + bi, a − bi) This implies some structure for S, hence for H This structure is a consequence of the solution of the inverse eigenvalue

problem. It is sufficient to know that the upper triu(S−1, 1) has rank 1.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 59

58/63

Algorithm

Update can be done recursively in O(n) operations ⇒ total of O(n2)

perations = fast algorithm.

Possibly O(nN) if degree n with N data points.     1 QH         × × × × × × × ×            1 1 Q        =     × × × S + A     S + A =

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 60

59/63

     w4 z4 ∗ M      =     w4 z4 ∗ u1¯ v1 + α0 v1¯ u2 v1¯ u3 u2¯ v1 u2¯ v2 + α1 v2¯ u3 u3¯ v1 u3¯ v2 u3¯ v3 + α2     Q =

c

s −s c

,

QH

w4

∗

=
∗
updates:

u′

4

    similar for the v’s and α’s. Applying QH I z4 M Q I

gives

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

d = u′′

2¯

v′′

2 + α′ 1 ⇒ α′ 1 = d − u′′ 2¯

v′′

2

Is SS+D, but (2,2) diagonal element α′

1 should be α1.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 62

61/63

Needs another unitary similarity transform on rows/columns 2 and 3. We can use a matrix of the form Q =

1

√

1+|t|2

1

¯ t −t 1

where t =

α1−α′

1

¯ v2u′′

2−¯

u2v′′

2

gives updates: QH

u′′

2

u2

=
u′

2

u′′

3

,

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 63

62/63

First versions implemented in matlab (M. Van Barel). Fast and numerically stable. Algorithms with complex poles yet to be designed. We may expect that

complex “conjugate” poles can be combined into 2 × 2 real blocks which

means that we can construct real basis functions, which will guarantee that the poles and zeros will have symmetry.

generalization to MIMO systems become possible.

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel

SLIDE 64

63/63

References SS matrices

[1] D. Fasino, L. Gemignani. Direct and inverse eigenvalue problems for diagonal-plus- semiseparable matrices (2001) Submitted [2] D. Fasino, L. Gemignani. A Lanczos-type algorithm for the QR factorization of Cauchy-like matrices (2002) Submitted [3] M. Van Barel, D. Fasino, L. Gemignani, N. Mastronardi. Orthogonal rational functions and diagonal-plus-semiseparable matrices, (2002) Submitted

M´ etodos de aproximaci´

n en teor´

ıa de sistemas, Laredo, September 2002

A. Bultheel