Rational bases for system identification Adhemar Bultheel - PowerPoint PPT Presentation

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´ on en teor´ ıa de sistemas, Laredo, September 2002

1/63 Summary 1. Definitions from signals, systems 2. Vector orthogonal polynomials (Jacobi/unitary Hessenberg matrices) 3. Orthogonal rational functions in L 2 ( R , µ ) , L 2 ( T , µ ) 4. ORF and inverse (Hessenberg) eigenvalue problem (semiseparable matrices) M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

2/63 Definitions • A signal is a complex function of a time variable. -discrete ( f n ) , n ∈ Z (discrete signal) or -continuous f ( t ) , t ∈ R (analog signal) • Define a delay operator ( D f ) n = f n − 1 ( D f ( t ) = f ( t − 1) ) • A pulse δ is the Kronecker delta (Dirac delta function) • A pulse decomposition: m f m D m ) δ f ( s ) D s δ ( · ) ds ) � f = ( � ( f = • An inner product: � f, g � = � f n g n � ( � f, g � = f ( t ) g ( t ) dt ) • The convolution h = f ∗ g : � h n = � m f m g n − m ( h ( t ) = f ( s ) g ( s − t ) ds ) M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

3/63 • The (auto)correlation function: � ( r fg ) n = � m f m g n + m ( r fg ( t ) = f ( s ) g ( s + t ) ds ) � ( r ff ) n = r n = � m f m f n + m ( r ff ( t ) = r ( t ) = f ( s ) f ( s + t ) ds ) • The Fourier transform: F = F f F ( e iω ) = � m f m e − imω f ( t ) e − itω dt ) � ( F ( ω ) = • The Z -transform: m f m z − m f ( t ) e − itz dt ) � F ( z ) = � ( F ( z ) = f ( t ) e − tp dt ) � (compare with (2-sided) Laplace transform F ( p ) = F ( z ) | z = e iω ∈ T = ( F f )( e iω ) F ( z ) | z = ω ∈ R = ( F f )( ω ) • power spectrum: R = F r = | F | 2 white noise: R = constant; normalized constant = 1 colored noise: g = H f , f white noise. M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

4/63 • A (linear) system H is a (linear) operator that maps a signal (input) to another signal (output) H u = y • It is stationary or time invariant if HD = DH • The impulse response: h = H δ y = H u = H ( � u m D m ) δ = � u m D m ( H δ ) = ( � u m h n − m ) = h ∗ u u ( s ) D s � u ( s ) D s ( H δ ) = �� � � y = H u = H δ = u ( s ) h ( · − s ) = h ∗ u y = h ∗ u = u ∗ h = � h m D m u = H ( D − 1 ) u hence H = H ( D − 1 ) . • The transfer function: H ( z ) = ( Z h )( z ) = | H ( z ) | e i Φ( 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 � < ∞ ⇒ �H u � < ∞ M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

5/63 • For a causal stable system: H ∈ H 2 ( E ) ( H ∈ H 2 ( 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´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

6/63 State Space If H is proper rational then partial fraction expansion: c 1 c n = C ( zI − A ) − 1 B + D H ( z ) = D + + · · · + z − α 1 z − α n where B = [1 , . . . , 1] T , A = diag( α 1 , . . . , α n ) , C = [ c 1 , . . . , c n ] , 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´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

7/63 If T is an invertible matrix, then ( A, B, C, D ) ≡ ( A ′ , B ′ , C ′ , D ′ ) if � − 1 A ′ B ′ TAT − 1 � � � � � � � � � T 0 A B T 0 TB = = C ′ D ′ CT − 1 0 1 C D 0 1 D since C ( zI − A ) − 1 B + D = C ′ ( zI − A ′ ) − 1 B ′ + D ′ . A minimal realization: size A is minimal. Then σ ( A ) = poles of the transfer function. H ( z ) = � h k z + · · · + CA k − 1 B z k = D + C ( zI − A ) − 1 B = D + CB + · · · . z k Y ( z ) = [ C ( zI − A ) − 1 B + D ] U ( z ) = CX ( z ) + DU ( z ) with X ( z ) = ( zI − A ) − 1 BU ( 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´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

8/63 is the z-domain formulation of x ( t ) = 1 dx ( t ) x n +1 = Ax n + Bu n x ( t ) = Ax ( t ) + Bu ( t ) ˙ y ( t ) = Cx ( t ) + Du ( t ) ˙ y n = Cx n + Du n i 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 n u n A B n y C D - many other advantages M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

9/63 Identification problem U 0 Y0 G N N U Y Y U Y = G ( U − N U ) + N Y = GU + ( N Y − GN U ) = GU + V , V = H 0 E , E white noise, H 0 some coloring filter. It is assumed we know Y and U (or G ) in a number of frequency points: z k ∈ T or R and some information about the noise (covariances) M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

10/63 Frequency Domain Identification Estimate G as ˆ G such that we minimize � Y − ˆ Y � w y or � G − ˆ G � w g where w = � N � F � 2 k =1 w 2 k | F ( z k ) | 2 . Note with w g = w y U : � Y − ˆ Y � w y = � GU − ˆ GU � w y = � ( G − ˆ G ) U � w y = � G − ˆ G � w g , E.g. measurements { G ( z j ) } N j =1 , variance { σ j } N j =1 , then set w g = 1 /σ . � � � � Y U − B Y A − BU w = w g � G − ˆ � � � � G � w g = = = � Y A − BU � w , � � � � A AU AU � � � � w g w g In ML estimation: w − 2 = Cov( AY − BU ) = σ 2 Y | A | 2 + σ 2 U | B | 2 − 2Re( σ 2 Y U AB ) . M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

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´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

12/63 Polynomial method min C = min � Y A − UB � 2 w w may depend upon (estimated) solution . Set W = w [ Y − U ] , P = [ A B ] T ⇒ w ( Y A − UB ) = WP N C = � WP � 2 = � P � 2 � P ( z j ) H W H P ∈ P 2 × 1 W = j W j P ( z j ) , . n j =1 � f, g � W = � f ( z j ) H W H j W j g ( z j ) a p.d. inner product in P 2 × 1 if n ≤ N . n min � P � W needs constraint: e.g. P “monic”. How to parametrize P ? M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel

13/63 � 1 z − 1 z − n � · · · 0 0 · · · 0 z − 1 z − n 0 0 · · · 0 1 · · · dim P 2 × 1 = 2( n + 1) : or � 1 z − 1 z − n � 0 0 · · · 0 z − 1 z − n 0 1 0 · · · 0 Simpler to orthogonalize { I, z − 1 I, . . . , z − n I } w.r.t. matrix-valued SP � f, g � W = � f ( z j ) H W H j W j g ( z j ) in P 2 × 2 Gives orthogonal block polynomials φ k : � φ k , φ l � W = δ k − l I to make unique choose e.g. l.c. upper triangular. Write P = � n k =0 φ k λ k = Φ( z ) λ , φ k ∈ P 2 × 2 , λ k ∈ C 2 × 1 , λ n � = 0 . k W = � Φ λ � W = λ H Φ H W H WΦ λ , Φ = [ φ j ( z i )] , W = diag( W k ) . � P � 2 φ k orthonormal ⇒ Φ H W H WΦ = Q H Q = I , Q = WΦ is unitary W = min � λ � 2 = λ H min � P � 2 n λ n ⇒ P = φ n λ n . M´ etodos de aproximaci´ on en teor´ ıa de sistemas, Laredo, September 2002 A. Bultheel