Spectral methods to compute a solution to some H interpolation - PowerPoint PPT Presentation

Spectral methods to compute a solution to some H ∞ interpolation problems A. E. Frazho The talk is mainly my approach on using FFT and spectral methods to compute solutions to H ∞ problems. It bypasses many state space methods. Some of the methods are standard and nothing new is claimed.

A review of FFT Consider a trigonometric polynomial � a k e − i ω k p ( e i ω ) = k The DFT is the Vandermonde matrix F ν on C ν containing { e − i ω j }         p ( e i ω 0 ) p ( e i ω 0 ) a 0 a 0         p ( e i ω 1 ) p ( e i ω 1 ) a 1 a 1          .   .      . . . . . = F − 1 .   = F ν   and     . . . . ν                 p ( e i ω ν − 2 ) p ( e i ω ν − 2 ) a − 2 a − 2 p ( e i ω ν − 1 ) p ( e i ω ν − 1 ) a − 1 a − 1 is ν evenly spaces points in [0 , 2 π ). If ν = 2 n the Here { ω j } ν − 1 0 computation is ν log( ν ) and called FFT. The inverse DFT is F − 1 = 1 ν F ∗ ν ν

Approximating Fourier series and H 2 functions + → L 2 be the Fourier transform: Let F : ℓ 2 � � � ⊤ α k e − ik ω = F f ( e i ω ) = · · · α − 1 α 0 α 1 · · · For continuous f evaluate f ( e i ω ) on 2 n points of [0 , 2 π ): take ifft     f ( e i ω 0 ) a 0     f ( e i ω 1 ) a 1         . . = F − 1 . . a =     . . ν         f ( e i ω ν − 2 ) a − 2 f ( e i ω ν − 1 ) a − 1 � 2 π a k ≈ α k = 1 e ik ω f ( e i ω ) d ω ifft computes Reiman intergal 2 π 0 � f � ∞ ≈ � f ( e i ω j ) � ∞ , C ν � f � 2 ≈ � a � C ν and f ∈ H 2 ⇐ ⇒ the bottom half of a equals zero Operator theory uses e ik ω . Matlab, engineers e − ik ω ( λ = z − 1 ).

d ( z ) = � a k z − k Transfer functions G ( z ) = n ( z ) For example, 0 . 84 z 2 + 1 . 88 z + 0 . 66 G ( z ) = z 4 + 0 . 44 z 3 − 0 . 43 z 2 − 0 . 056 z + 0 . 014 = 0 . 84 + 1 . 5 z 3 + 0 . 35 + 0 . 53 + · · · z 2 z 5 z 6 Matlab commands: 2 16 = 65536 is overkill. g = fft([0 , 0 , 0 . 84 , 1 . 88 , 0 . 66] , 2 16 ) ./ fft([1 , 0 . 44 , − 0 . 43 , − 0 . 056 , 0 . 014] , 2 16 ); m = real(ifft( g )); � 0 0 . 53 � m (1: 6) = 0 0 . 84 1 . 5 0 . 35 norm( m (2 15 : 2 16 ) = 3 . 64 × 10 − 16 ; Therefore G ( z ) is stable norm( g , inf) = 3 . 47 = � G � ∞ and norm( m ) = 1 . 87 = � G � 2

The Nyquist and Magnitude plot Figure: The plot of G ( e i ω ) and | G ( e i ω ) | ; winding number = 3.

Stability of polynomials z n p ( z ) is in H 2 . The The polynomial p ( z ) is stable if and only if polynomial � � 5 � z − k p ( z ) = ( z − 2) 10 k =1 = z 6 − 3 . 5 z 5 + 3 . 85 z 4 − 1 . 925 z 3 + 0 . 4774 z 2 − 0 . 056 z + 0 . 0024 is unstable. Matlab commands: (2 20 = 1048576 is way overkill) p = poly([(1 : 5) / 10 , 2]); f = 1 ./ fft( p , 2 20 ); m = ifft( f ); norm( m (2 19 : 2 20 )) = 1 . 324 ( unstable ) norm(ifft( f )) = 2 . 3625 = � f � 2 (ifft( f ) , inf) = 1 . 2933 = � f � ∞

Hankel approximation: Classical Kalman-Ho Consider a rational stable ∞ � G ( z ) = n ( z ) a n z − n d ( z ) = n =0 The corresponding Hankel matrix     a 1 a 2 a 3 · · · C     · · · a 2 a 3 a 4 CA � �     A 2 B H =   =    B AB · · · CA 2 a 3 a 4 a 5 · · ·    . . . . . . . . . . . . . . . rank H = dim of minimal realization of G ( z ) = D + C ( zI − A ) − 1 B Using the finite partition of H with the svd the Kalman-Ho algorithm computes { A , B , C , D } .

Hankel approximation: example Consider the non rational function 2 z + 1 f ( z ) = e z 2 Using the FFT with Kalman-Ho f ( z ) ≈ g = z 4 + 0 . 5039 z 3 + 0 . 1091 z 2 + 0 . 0123 z + 0 . 0006292 z 4 − 0 . 4961 z 3 + 0 . 1052 z 2 − 0 . 01152 z + 0 . 0005635 � f � ∞ = 20 . 0855 and � f � 2 = 6 . 9435 and � f − g � ∞ = 0 . 0082 As expected, g ( z ) is outer. Matlab commands: f = exp(fft([0 , 2 , 1] , 2 16 )); m = real(ifft( f )); [ a , b , c , d ] = kalho( m (1 : 700) , . 01); [ n , dd ] = ss2tf( a , b , c , d ); g = fft( n , 2 16 ) ./ fft( dd , 2 16 ); norm( f , inf ) = 20 . 085 = � f � ∞ norm( m ) = 6 . 94 = � f � 2 , � f − g � ∞ = 0 . 0082

Toeplitz matrices If r = � ∞ −∞ e − ik ω ∈ L ∞ , then T r is the Toeplitz matrix:   r 0 r − 1 r − 2 · · ·   r 1 r 0 r − 1 · · ·    on ℓ 2 T r = and � T r � = � r � ∞   r 2 r 1 r 0 · · · +  . . . ... . . . . . . If θ = � ∞ 0 θ k z − k is in H ∞ , then   θ 0 0 0 · · ·   θ 1 θ 0 0 · · ·    on ℓ 2 T θ =   θ 2 θ 1 θ 0 · · · +  . . . ... . . . . . . ◮ θ ∈ H ∞ is outer if ran( T θ ) = ℓ 2 + . ◮ θ is an outer spectral factor for T r if θ is outer and T r = T ∗ θ T θ ◮ If T r ≫ 0, then T r admits a unique outer spectral factor θ .

Outer spectral factorization Let T r ≫ 0 be Toeplitz. Solve √ e � � ⊤ = � � ⊤ and θ ( z ) = � ∞ T r 1 a 1 a 2 · · · e 0 0 · · · 0 a k z − k θ is outer, T r = T ∗ θ T θ . Let T n the upper n × n corner of T r . Solve � � ⊤ = � � ⊤ 1 · · · 0 0 · · · 0 T n a 1 a 2 a n − 1 e √ ez n − 1 θ n ( z ) = z n − 1 + a 1 z n − 2 · · · + · · · a n − 1 θ ( z ) = lim n →∞ θ n ( z ) � � ◮ T n = T ∗ n and θ n ( z ) is outer θ n T θ n ◮ θ n ( z ) can be computed fast by Levinson, n = 10000 is easy. ◮ If r is rational, then M deg ( θ ) < ∞ . ◮ M deg ( θ n ) = n − 1 → ∞ . The Kalman-Ho finds θ from θ n . ◮ The Matrix case is similar.

Inner outer factorization g = g i g o : an example The matrix case is similar. Pole-zero cancelation unstable. ◮ Form the Toeplitz T r from | g | 2 . Use Levinson to find θ n ≈ g o . ◮ Compute g i = g g o . ◮ Apply Kalman-Ho to find state space realizations. 0 . 84 z 2 + 1 . 88 z + 0 . 66 g ( z ) = z 4 + 0 . 44 z 3 − 0 . 43 z 2 − 0 . 056 z + 0 . 014 1 . 5 z 4 + 1 . 5 z 3 + 0 . 37 z 2 0 . 557 z + 1 g o ( z ) = and g i ( z ) = z 4 + 0 . 44 z 3 − 0 . 43 z 2 − 0 . 0555 z + 0 . 014 z 2 (1 + 0 . 557 z ) g = fft(num , 2 16 ) ./ fft(den , 2 16 ); r = real(ifft(abs( g ) . 2 )); [ a , e ] = levinson( r (1 : 1000)); go = sqrt( e ) ./ fft( a , 2 16 ); gi = g ./ go ; mo = real(ifft( go )); mi = real( ifft ( gi )); [ ao , bo , co , dko ] = kalho( mo (1 : 700)); [ no , do ] = ss2tf( ao , bo , co , dko ); Go = tf( no , do ) [ ai , bi , ci , dki ] = kalho( mi (1 : 700)); [ ni , di ] = ss2tf( ai , bi , ci , dki ); Gi = tf( ni , di )

A classical outer factorization formula: scalar case � � � 2 π z + e i ω 1 z − e i ω ln( | g ( e i ω ) | d ω g o ( z ) = exp 2 π 0 Change the form for FFT. Notice that � � ∞ ∞ a n e − in ω = h + h and h = a 0 2 ln | g ( e i ω ) | = a n e − in ω 2 + n = −∞ n =1 Claim: g o ( z ) = e h which is suitable for FFT | g ( e i ω ) | 2 = e 2 ln | g ( e i ω ) | = e h e h = g o g o Same example: 0 . 84 z 2+1 . 88 z +0 . 66 g ( z ) = z 4+0 . 44 z 3 − 0 . 43 z 2 − 0 . 056 z +0 . 014 g = fft(num , 2 16 ) ./ fft(den , 2 16 ); a = ifft(2 log(abs( g ))); gc = exp(fft([ a (1) / 2 , a (2 : 2 15 )] , 2 16 )); % this ∼ e h norm( gc − go , inf ) = 8 . 9543 × 10 − 15 = � g c − g lev � ∞ 1 . 5 z 4 + 1 . 5 z 3 + 0 . 37 z 2 g o ( z ) = z 4 + 0 . 44 z 3 − 0 . 43 z 2 − 0 . 0555 z + 0 . 014

The band formula: Toeplitz extension Many matrix H ∞ interpolation problems similar calculations. Let T n +1 ≫ 0 be Toeplitz on C n +1 with entries { r j } n 0 . Carath´ eodory-Toeplitz interpolation: Find all strictly positive Toeplitz extensions T r ∼ { r j } of T n +1 . The band method solution: � � ⊤ = � � ⊤ 1 · · · 0 0 · · · 0 T n +1 a 1 a 2 a n e √ eu = 1 + a 1 z − 1 + a 2 z − 2 + · · · a n z − n All strictly positive Toeplitz extensions T r ∼ { r j } of T n +1 are ∞ � 1 − | g | 2 r k z − k = | θ | 2 r = | z − n ug + u | 2 = ( g ∈ S ◦ ) k = −∞ How to compute { r k } and outer θ such that T r = T ∗ θ T θ .

A band method example: Let T 4 the Toeplitz from { 10 , 9 , 8 , 6 } and 1 . 275 z − 1 . 797 g ( z ) = and � g � ∞ = . 9 0 . 9 z 3 − 1 . 068 z 2 − 0 . 3655 z + 0 . 5393 Some Matlab commands: 2 17 = 131072 is overkill. [ a , e ] = levinson([10 , 9 , 8 , 6]); u = fft( a , 2 17 ) ./ sqrt( e ); r = (1 − abs( g ) . 2 ) ./ (abs(conj( u ) . ∗ g . ∗ fft([0 , 0 , 0 , 1] , 2 17 ) + u ) . 2 ); � � r j = real ( ifft ( r )); r j (1 : 5) = 10 9 8 6 4 . 25 [ q , e ] = levinson( r j (1 : 2000)); θ = sqrt( e ) ./ fft( q , 2 17 ); norm ( r − abs ( θ ) . 2 , inf) = 2 . 2021 × 10 − 10 = � r − | θ | 2 � ∞ Applying Kalman-Ho to ifft( θ ) yields the outer factor 1 . 106 z 6 − 1 . 335 z 5 − 0 . 4449 z 4 + 0 . 677 z 3 θ ( z ) = z 6 − 2 . 103 z 5 + 0 . 1902 z 4 + 2 . 128 z 3 − 1 . 041 z 2 − 0 . 5014 z + 0 . 328