In this paper, we present a novel non-iterative algorithm to identify linear time-invariant systems from frequency response data

Full text

(1)Ecient Construction of Transfer Functions from Frequency Response Data Tomas McKelvey, Huseyin Akcay, and Lennart Ljung Department of Electrical Engineering Linkoping University S-581 83 Linkoping, Sweden, Tel: +(46) 13 282461, Fax: +46 13 282622 Email: tomas@isy.liu.se September 15, 1994 Report: LiTH-ISY-R-1609. Submitted to IEEE Trans. of Automatic Control for review. Abstract. In this paper, we present a novel non-iterative algorithm to identify linear time-invariant systems from frequency response data. The algorithm is related to the recent time-domain subspace identication techniques. Promising results are obtained when the algorithm is applied to the real frequency data originating from a large exible structure. A robustness analysis is performed and under weak conditions on the measurement noise, it is shown that the algorithm is consistent.. 1 Introduction Recently, identi cation and control of large space structures has received considerable attention. 18, 19, 20, 4, 14, 13]. This type of systems are also frequently encountered in the modal analysis area of mechanical engineering. Typically such systems are lightly damped and quite often as in the system analysis and control design of mechanical structures, high order models with many outputs are needed. For structural design purposes, the nite element method provides accurate enough models. Then static and dynamic tests on the structure are performed to re ne the nite element model. However, this traditional approach to model development may not be accurate enough if the intended use of the analytical model is to design a control system since most modern multivariable control design techniques are based on the state-space models of the systems. A direct method is Author to whom all correspondence should be sent to.. 1.

(2) then to realize the model from the experimental results. Identi cation of multi-input multi-output (MIMO) systems of high orders is a challenging problem. If time-domain measurements are available, a great number of algorithms exist in the literature. These algorithms can be classi ed either as iterative or noniterative. Among the iterative algorithms, we nd the prediction error methods 21, 31] and among the noniteratives, we nd the more recent subspace based algorithms 23, 34, 26]. Noniterative methods do not involve nonlinear parametric optimization. In particular, subspace based algorithms deliver state-space models without the need for an explicit parametrization. Essentially, there is no dierence between MIMO system identi cation and single-input single-output (SISO) system identi cation for a subspace based algorithm. They are more robust to inaccuracies than the parametrized canonical models due to the fact that there is no need for explicit parametrization. In 35, 24], subspace based algorithms were shown to be consistent and asymptotic expressions for the quality of the estimates were derived in. 35]. In practice, information about a system is often characterized in terms of the frequency response of the system at some discrete set of frequencies. Each transfer function measurement is compiled from an enormous number of time-domain measurements using high performing sophisticated data analyzers and data acquisition equipment and is thus a condensed version of the crude measurements. In this step a signi cant noise reduction can be expected if the excitation of the system is well designed. (See 2, 5, 28] for a discussion on the data acquisition and the statistical properties of the frequency response data). Some experimental studies 3] suggest performing identi cation in the frequency domain for large exible structures due to the poor t in the time-domain. The problem of tting a real-rational model to a given frequency response has been addressed by many authors, e.g. 30, 22, 28]. In the traditional way, a system is modeled as a fraction of two polynomials with real coecients and a nonlinear least-squares t to the frequency response data as in the prediction error methods for the time-domain problem is sought. In an early result 16], and later further re ned in 29], a sequence of linear least-squares problems called SK-iterations are solved. However, SK-iterations are not guaranteed to terminate at the global minimum as indicated in 36]. A second drawback is the sensitivity of the poles and zeros of the system to polynomial factoring if the system order is high. This de ciency can be reduced by introducing other parametrizations e.g. orthogonal Chebyshev polynomials 9, 2], zero-pole-gain form, or the related RPM-parametrization 25]. The algorithm in 14] is noniterative and based on the famous Ho and Kalman realization algorithm 12] or Kung's smoothed version 15]. In 14], the impulse response coecients of the system also called Markov parameters are estimated applying the inverse discrete Fourier transform (DFT) to the frequency response data. The realization algorithms 12, 15] nd a minimal statespace realization given a nite sequence of Markov parameters. In 14], a recursive scheme to calculate the estimates of the Markov parameters is proposed. This approach is exact only if the impulse response dies out completely within the number of given frequency points, in other words if the system has a nite impulse response and therefore for lightly damped systems yields very poor estimates. In order to perform the inverse DFT, the frequency data must also be uniformly spaced. In 4], Bayard suggests rst tting a high order rational model to the data using the SK-iterations and then calculating Markov parameters of the high order model. Next, a reduced order model in the state-space form is obtained utilizing the realization algorithm of Ho and Kalman 12]. A new frequency domain approach proposed by Liu and co-workers 18, 19] is a frequency domain counterpart of the time-domain subspace methods by Moonen 23] and Liu 20]. This 2.

(3) approach does not require the data to be uniformly spaced in frequencies and also oers some frequency weighting capabilities. We will now outline the contents of this paper. In Section 2, we formulate the general problem and present formulas to compute state-space realizations from frequency domain data. We consider the subspace formulation of 18], where frequency-domain data are samples of the Fourier transforms of input and output data or spectrum samples when wide-sense stationary inputs are applied. In this section, we show that all these formulations are dierent implementations of the familiar timedomain subspace algorithms. In Section 3, we present the main result of the paper: The inverse DFT using only a nite number of uniformly spaced frequency response samples combined with a subspace identi cation step indeed yields the true system in spite of the bias, introduced by the inverse DFT technique and possessed by the estimated impulse response coecients. As a result, it will follow that the approach of Liu et.al. 18] is closely related to our algorithm if the frequencies are uniformly spaced and equal if no noise is present. To illustrate the main result, in Section 4, we test our algorithm on the real data originating from a frequency response experiment on a exible structure testbed at the Jet Propulsion Laboratory, Pasadena, California. The example clearly indicates that the subspace methods are competitive compared with the iterative least-squares methods. In Section 5, we perform a sensitivity and robustness analysis of the proposed algorithm to model imperfections, which provides some clear insights on how user choices aect the resulting estimates. An asymptotic noise analysis is pursued in Section 6 showing the consistency of the estimates. Section 7 contains the conclusions.. 1.1 Notation. j. AT A AH Ay tr(A) jjAjjF (A) i (A) AB. ;i (G). H1. kGk1. p ;1. the transpose of A the complex conjugate of A the complex conjugate transpose of A (AT A);1 AT the Moore-Penrose pseudo inverse of the full-rank matrix A. P a , the trace of A q i ii H tr(A A), the Frobenius norm of A the spectral radius of A the ordered singular values of A, 1 : : : n the Kronecker product 2 a11 B a12 B a1cB 3 64 ... .. .. .. 75 . . . ar1 B ar2 B arc B when A is a r c matrix the ordered Hankel singular values of G, ;1 : : : ;n Hardy space of matrix-valued bounded analytic functions in the open unit disc of the complex plane sup norm of G 2 H1 , = sup! (G(ej! )) 3.

(4) 2 State Space Model Identication in Frequency Domain Assume that the true system is a stable, nite-dimensional, linear time-invariant, discrete-time system. Then its input/output properties can be described by a pair of state-space equations. x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k). (2.1). where u(t) 2 IRm , y(t) 2 IRp and x(t) 2 IRn . Note that all such pairs describing the same input/output behavior of the system are equivalent under a similarity transformation. In this paper, we consider the problem of nding a state-space realization (2.1), given a nite number of noisy frequency response data of the system. Gk = G(ej!k ) + ek !k 2 0 ] k = 1 M:. (2.2). 2.1 Observability Range Space Extraction. In this section, we present formulas to compute observability range space from frequency domain data. This section provides background for the next section where we specialize to frequency response data and introduce an ecient and robust algorithm for uniformly spaced data. Most formulas in this section were presented in a recent paper 18]. Our contribution is Subsection 2.2, where a connection between the algorithm presented in 18] and the more familiar time-domain subspace identi cation algorithms is made. Assume that the unknown system is described by the state-space equations (2.1). We also assume that the system is minimal and thus both observable and controllable. Concatenating outputs and inputs as. 2 3 2 3 y (k ) u(k) 66 y(k + 1) 77 66 u(k + 1) 77 7 66 77 yq (k) = 66 u ( k ) = .. .. 75 q . . 4 4 5 y(k + q ; 1) u(k + q ; 1). (2.3). and introducing the extended observability and the block Toeplitz matrices. 2 C 66 CA Oq = 66 .. 4 .. CAq;1. 3 2 D 77 66 CB 77 ;q = 66 .. 5 4 .. ::: ::: ... :::. 0 D .. . CAq;2B CAq;3B. we obtain the relation. yq (k) = Oq x(k) + ;q uq (k): Furthermore, let us construct the following block Hankel matrices. h. 3 77 77 5. (2.4). (2.5). i. YqN = yq (1) yq (2) yq (N ; q) 4. 0 0 .. . D. (2.6).

(5) h. UqN = uq (1) uq (2) uq (N ; q) and concatenate the states as follows. h. XN ;q = x(1) x(2) x(N ; q). i i. (2.7) (2.8). which together with (2.5) results in the matrix equation. YqN = Oq XN ;q + ;q UqN : (2.9) The subspace identi cation algorithms 23, 34] and others are based on the equality (2.9). After certain projections, a matrix whose column range space is the same as the range space of the columns of extended observability matrix is obtained. Then utilizing a singular value decomposition, a state-space model is extracted. It turns out that the same idea also works in the frequency domain given samples of the Fourier transform of the input and output data at some set of frequencies. By taking the Fourier transform of both sides of (2.5), we obtain (2.10) W (!) Y (!) = Oq X (!) + ;q W (!) U (!) where denotes the Kronecker product and. h. W (!) = 1 ej! ej2! ej!(q;1). iT. :. (2.11). Assume that Y (!) and U (!) are given at a discrete set of frequencies !k 2 0 ], k = 1 : : : M . Let us now use these samples to de ne the following matrices. h. i. WqM = W (!1 ) W (!2) W (!M ) . 2 3 Y (!1 ) 0 0 66 .. 77 0 Y ( ! . 77 6 2) : : : YqM = (WqM Ip) 66 . . .. 7 . .. .. . 5 4 .. 0 Y (!M ) 2 3 U (!1 ) 0 0 66 .. 77 0 U ( ! ) : : : . 77 6 2 UqM = (WqM Im ) 66 . . . . .. .. .. 75 4 .. 0 U (!M ) h i XqM = X (!1 ) X (!2 ) X (!M ) :. (2.12). (2.13). (2.14) (2.15). With the matrices above, we arrive at the complex matrix equation. YqM = Oq XqM + ;q UqM : 5. (2.16).

(6) This equation is the frequency domain counterpart of the time domain equation (2.9). Since Y (!) and U (!) are computed from the input/output data by a numerical procedure, e.g., the fast Fourier transform (FFT), it is reasonable to assume that equidistant frequencies are used in the computation of Y (!k ) U (!k ) k = 1 : : : M . The equidistance property will turn out to be quite fruitful in the next section. However in this section, we do not impose this condition on !k k = 1 : : : M . Instead of samples of the Fourier transform of input/output data, assuming u is a wide-sense stationary process, we can use spectrum samples de ned by. Syq u (!) = E yq (!)u(!) ] Suq u(!) = E uq (!)u(!) ] Sxu(!) = E x(!)u(!) ]. (2.17). where E denotes expected value. Then one obtains the following analogue of (2.10). W (!) Syq u(!) = Oq Sxu(!) + ;q W (!) Suq u(!):. (2.18). as described in 18]. Therefore in the discussion of observability range space, without loss of generality, we will consider only (2.10) whether the Fourier transform or spectrum samples are used.. 2.1.1 Estimates of B and D. In the frequency response of (2.1), computed as. G(ej! ) = C (ej! In ; A);1 B + D. (2.19). we notice that both B and D appear linearly assuming A and C given. Therefore, given state-space parameters A and C , we can calculate B and D by solving the following linear least-squares problem M . ^ D^ = arg min X Y (!k ) ; D + C (ej!k I ; A);1 B U (!k )2 : (2.20) B F BD k=1. Since input/output cross spectrum and input autospectrum are related as. h. i. Syu(!) = C (ej! I ; A);1 B + D Suu(!). (2.21). B and D can be calculated from an analogous equation to (2.20).. 2.2 Observability Range Space Extraction in Frequency Domain. There are several methods to extract the observability range space from (2.16). In (2.16), the real matrices Oq , ;q and complex matrix X are unknown. We begin with the construction of the following three matrices. U = UqM UqM ] X = XqM XqM ] Y = YqM YqM ]: Then, from (2.16), we obtain. Y = Oq X + ;q U : 6. (2.22) (2.23).

(7) The rst method we employ is closely related to the time domain work in 23] and we form the following matrix U ? = I ; U H (UU H );1 U (2.24) assuming that U has full rank and multiply (2.23) from right by U ? to obtain. Y U ? = Oq X U ?: (2.25) Notice that Rank (U ?) = 2M ; qm since U ? is a projection matrix on IR2M . Hence it has maximal rank among all annihilators of U . Multiplication of (2.24) from right by U H yields. U ?U H = U H ; U H (UU H );1 UU H = U H I ; (UU H );1 UU H = 0 (2.26) The relation (2.26) is used in 18] to arrive at. YY H ; YU H (UU H );1 (YU H )H = Oq XX H ; XU H (UU H );1 (XU H )H OqT . which also can be written as. Y U ?Y H = Oq X U ? X H OqT. (2.27) (2.28). since U ?(U ?)H = U ?. The imaginary parts of the above equation are zero. Since (2.25) and (2.27) are equivalent, we could have stopped at (2.25). However, typically M q which provides data compression prior to the extraction of A and C from the observability range space. The second method we consider was described for the time domain problem by Verhaegen 34]. Applying an LQ factorization as follows. ". # " #" H # U = L11 0 Q1 Y L21 L22 QH2. (2.29). and using (2.25), we can obtain the observability range space from. L22 = Oq X Q2 :. (2.30). Thus all the methods employed are equivalent but dierent implementations. One might be tempted to say that one of the implementations above is to be preferred over the others for noisy data. However, this is not the case since the state-space parameters A and C are obtained by means of an SVD, which smoothes and compresses the data. It seems that the implementation (2.27) has less storage requirements than the others if M q. In order that YU ? has the same column range space as Oq , certain conditions on M and UqM must be imposed which we will not dwell on. We have already assumed that rank U = qm . However, in the next section, we will state a condition on the number of data when frequency response data are available. We have yet to discuss how to compute A and C given observability range space. This will be the subject of the next section.. 7.

(8) 3 State Space Identication from Frequency Response Data In this section, we assume that frequency response data G(ej!k ) on a set of of frequencies !k 2 0 ] k = 1 : : : M are given. Since G is a real-rational transfer function, frequency response data on 0 ] can be extended to 0 2] as follows. G(e;j! ) = G (ej! ) 0 ! : (3.1) and without loss of generality to simplify the notation, we may assume that frequency response data for w0 = 0 and wM = are also available and let wM +k = wM ;k k = 1 : : : M ; 1: (3.2). 3.1 Generation of frequency response data. There are mainly two methods to generate frequency response data. The rst method is to apply persistently exciting inputs to the system and sample the Fourier transforms of the input/output data. The second method is to apply wide-sense stationary inputs to the system and sample input autospectrum and input/output crosspectrum. Let B and D in (2.1) be partitioned columnwise as B = b1 b2 bm] D = d1 d2 dm ] (3.3) and let Gi (!) denote the frequency response of (A bi C di ): Suppose that at each identi cation experiment, only one input to a particular channel is applied and outputs at all channels are measured and experiments are performed exhausting all input channels. Let U i (!k ) and Y i (!k ) denote the Fourier transforms of the ith components of the input and the corresponding output data, respectively. Then the frequency response data are computed from i Gi (ej!k ) = UY i((!!k )) !k 2 0 ] k = 0 1 : : : M i = 1 : : : m: (3.4) k Similarly, assuming that inputs are spatially white, we can calculate frequency response data using spectrum samples as follows (3.5) Gi(ej!k ) = SSyui((!!)) !k 2 0 ] k = 0 1 : : : M i = 1 : : : m: u i ui k The frequency response data above are readily computed by microprocessor-based spectrum analyzers.. 3.2 Observability range space identication. Given frequency response data construct the following matrix. 2 i j!0 66 j!G0 (ei j!) 0 Giq2M := 666 e G. (e ) .. 4. ej!0(q;1) Gi (ej!0 ). Gi (ej!1 ) ej!1 Gi(ej!1 ) .. .. 8. . Gi(ej!2M ;1 ) .. . .. ... . ej!2M ;1(q;1) Gi(ej!2M ;1 ). 3 77 77 75. (3.6).

(9) for i = 1 : : : m , where the entries of Giq2M comply with (3.1). Using the notation Xqi 2M for the sampled state frequency response matrix, we can derive the following equation as a special case of (2.16) Giq2M = Oq Xqi 2M + ;iq Wq2M i = 1 : : : m (3.7) where the construction of ;iq bi and di are used. Since a Vandermonde matrix has full rank if frequencies are distinct and qr < M , Wq2M has a right annihilator Wq?2M and multiplying (3.7) from right by Wq?2M we get. Giq2M Wq?2M = Oq Xqi 2M Wq?2M i = 1 : : : m: From (3.8), we get the following relation between the range spaces Range. h. Gq2M W 1. ? q2M. . Gmq2M Wq?2M. In (3.9), the equality is attained if Rank. m i=1. m i . =. Range. i=1. (3.8). Range Giq2M Wq?2M Range Oq :. Giq2M Wq?2M. !. = n:. (3.9). (3.10). Thush if (3.10) is satis ed, we can extract i the observability range space from the column range space 1 ? m ? of Gq2M Wq2M Gq2M Wq2M : The major dierence from the the algorithm in 18] is the use of a low rank annihilator Wq?2M instead of the full rank annihilator (2.24). This dierence will be signi cant in the noise analysis in Section 5.. 3.3 Uniformly Spaced Frequency Response Data. Let us assume that equidistant frequency response data G(ejk=M ) k = 0 : : : M are available. Then, we can derive a simple yet powerful algorithm resulting from the equation (3.7) as follows. First, extend the data to satisfy (3.1) and observe that if maxfq rg < M , then Wq2M is annihilated by the following matrix. 2 66 1 ? Wqr = 2M 66 4. 1. 1. ej2=2M .. .. ej 2(2M ;1)=2M. . ej4=2M : : : .. ... .. . 1. ej2r=2M .. .. ej2r(2M ;1)=2M. 3 77 77 : 5. (3.11). ? is calculated as Furthermore, (i1 i2 )-th entry of Giq2M Wqr. 1 2MX;1 Gi (ej 2k=2M ) ej 2k (i2 +i1 ;1)=2M ? Giq2M Wqr (i1 i2 ) = 2M k=0. := g^i (i2 + i1 ; 1) i1 = 1 : : : q i2 = 1 : : : r 9. (3.12).

(10) which 2M -point DFT ofi the frequency response data. Since the column range space of h 1 is the Gq2M Wq?2M Gmq2M Wq?2M is invariant under the permutations of its columns, we can use G instead of Gi in the DFT computation above. Thus, we have demonstrated that the column range space of the following Hankel matrix. 2 r 66 gg^^21 gg^^32 :: :: :: g^rg^+1 H^ qr := 66 .. .. . . .. . 4 . . .. 3 77 77 5. (3.13). g^q g^q+1 : : : g^q+r;1 where the rst q + r impulse response coecients of G are estimated from the 2M -point inverse DFT as follows 1 2MX;1 G(ej 2k=2M ) ej 2ik=2M i = 0 : : : q + r ; 1 g^i := 2M (3.14) k=0 is a subspace of the observability range space if maxfq rg < M: Furthermore, if Rank (H^ qr ) = n then the column range space of H^ qr is equal to the observability range space. Now, we present our algorithm. Algorithm 3.1 1. Expand the given frequency data according to (3.1) and perform the 2M -point inverse DFT (3.14) on the expanded data to obtain the q +r ;1 estimates of the impulse response coe cients g^k . 2. Construct the q r block Hankel matrix (3.13). Perform a singular value decomposition of H^ qr as follows "^ #" ^T # 0 V1 1 ^ ^ ^ Hqr = U1 U2] 0 ^ (3.15) V^2T 2 where ^ 1 contains the n dominant singular values on the diagonal. 3. The system matrices are estimated as. h i. y h i. A^ = I(q;1)p 0(q;1)pp U^1 ^ 11=2 0(q;1)pp I(q;1)p U^1 ^ 11=2. (3.16). C^ = Ip 0p(q;1) U^1 ^ 11=2. (3.17). h. i. 2M ;1 ^ D^ = argmin X Gi ; G^ (eji=M )2 B F. where. i=0. ^ G^ (z ) = D^ + C^ (zI ; A^);1 B: 10. (3.18) (3.19).

(11) From (3.14), notice that lim g^ = M !1 k. Z 2 0. G(ej 2 ) ej2k d = gk k = 0 : : : q + r ; 1. (3.20). where gk is the kth impulse response coecient of G. Thus for the noise free data, H^ qr tends to a limit as M tends to in nity and the resulting algorithm is known as Kung's algorithm 15]. In Kung's algorithm, B^ is calculated from. B^ = ^ 11=2 V^1T. ". Im. #. 0m(r;1) :. (3.21). The realization given by (3.15), (3.17), and (3.21) is balanced i.e. the q-block row observability matrix Oq and the r-block column controllability matrix. h i Cr = B^ A^B^ : : : A^r;1B^. satisfy. (3.22). OqT Oq = Cr CrT = ^ 1. (3.23) As q and r jointly tend to in nity, the products (3.23) will converge to the observability and controllability Gramians, respectively, and the diagonal elements of ^ 1 will converge to the Hankel singular values of the system. Thus for nite q and r, the Hankel singular values of the system are underestimated by ^ 1 . Although ^ does not play any role in the construction of a state-space for G^ (except for a selection of a base for the states), it will be essential in the selection of a model order in the presence of unmodeled dynamics and noise. (Notice that A^ and C^ can be calculated from (3.15-3.17) letting ^ 1 = In ). As q and r increase, singular values of H^ qr also increase in the absence of noise. However, small singular values of H^ qr and corresponding left and right singular vectors will be less reliable in the presence of noise. Therefore q and r should be chosen suciently large to obtain ^ 1 as large as possible and model order should be kept lowest possible. These suggestions will be supported by the computations in Section 5. The results in this section are summarized in the following theorem.. Theorem 3.2 Let G be an nth order stable system represented by (2.1). Then n + 2 noiseless. equidistant frequency response measurements of G on 0 ] are su cient to identify G by the above algorithm 3.1.. Proof: See Appendix A ^ C^ can also be calculated by the total least-squares technique. Usually, little The parameters A or nothing is gained over the least-squares method due to the prior smoothing by the DFT and SVD.. 3.4 Nonuniformly Spaced Frequency Response Data and the Identication of Continuous-Time Systems. ? in (3.11). This Given equally spaced frequency response measurements, we easily obtained Wqr particular form is numerically well conditioned and hence robust to measurement noise. With. 11.

(12) ? will be ill conditioned if some frequencies are closely located and nonuniformly spaced data, Wqr therefore sensitive to noise. A direct subspace identi cation algorithm for the continuous-time systems could be developed ? would along the same lines of this paper. However in this case, numerical conditioning of Wqr ? for the discrete-time systems since the entries of the latter is be much worse than that of Wqr magnitude bounded by one. The Tustin and the zero-order hold sampling equivalence are two popular techniques to discretize continuous-time systems. The Tustin method requires measuring the frequency response at very high frequencies if the equivalent discrete-time frequency response is to be equally spaced. Hence the Tustin method is limited to the narrow bandwidth systems. Therefore in the design of identi cation experiments for large bandwidth systems especially for exible structures, the zero-order hold sampling equivalence method seems to be preferable for numerical considerations.. 4 Example: A Flexible Structure To illustrate the eectiveness of Algorithm 3.1, we will use the experimental frequency response data (JPL-data) obtained from the tests on exible structures at the Jet Propulsion Laboratory, Pasadena, California 4]. The Algorithm 3.1 will compete with the nonlinear least-squares algorithm outlined in Subsection 4.1 on the same data. The JPL-data have the following characteristics: 512 number of complex frequency samples Gk , k = 1 : : : 512 Frequency range !kc = 1:23 ; 628 radian/sec. equally spaced. Since the data origin from a exible structure, it contains many lightly damped modes. Our aim is to construct a discrete-time model matching the given frequency response. For the discrete-time models, we will map the continuous-time frequencies to the discrete-time frequencies as c k !k = ! c !N which gives the zero-order hold sampling equivalent of a continuous-time system when !Nc is chosen as the Nyquist frequency. In all estimates, we let models to be proper i.e. a feed-through term D is estimated.. 4.1 Nonlinear least-squares frequency identication. In this subsection, we will briey describe a least-squares method in the identi cation of a rational model G^ (z ) = ^b(z )=a^(z ) given the frequency response data Gk , to compare our algorithm with an existing one. This method is implemented in the Matlab Signal Processing Toolbox, 17], as the command invfreqz. The method is divided in two stages. First, an initial estimate is obtained by a linear least-squares solution, see 16]. X j!k a^0 (ej! ) ^b0 (ej! ) = arg min jb(e ) ; a(ej!k )Gk j2 ab N. with. k=1. a(z) = z na + ana ;1 zna ;1 + : : : + a0 12. (4.1).

(13) and. b(z ) = bnb znb + bnb ;1 znb;1 + : : : + b0 : This step is equivalent to a least-squares minimization of the residuals j!k "k = Gk ; ab((eej!k )) weighted by the estimated denominator a(ej!k ). Only if the frequency data is noise free and the model order is suciently large this step will yield an unbiased estimate. Using real data, a second step is thus necessary. The second step is a nonlinear least-squares minimization N b(ej!k ) X j! j! ^ ; Gk j2 (4.2) a^(e ) b(e ) = arg min j ab k=1 a(ej!k ) solved by a Gauss-Newton procedure initialized with the estimates a0 and b0 from the rst step. The second step removes the weights of the residuals introduced by the rst step. Since the second step involves a nonlinear optimization, the degree of success is completely dependent on the quality of the initial model a0 and b0 provided by the rst step.. 4.2 Model order selection by cross validation. Model validation is one issue that distinguishes the system identi cation from curve tting. In system identi cation, the measured data are assumed to be generated by some nite-dimensional system corrupted by some noise. Finding a good model order selection criteria is then an interesting problem since increasing model order usually decreases the estimation error (unless numerical problems occur or iterative methods are used and local minima are reached). If the model order is increased above the true order, the model will start to t the noise. In time-domain identi cation, a common solution to this problem is cross validation 33]: divide the data set into two disjoint sets and use one set for the identi cation and the other set for the model validation. If the order of the estimated model increases over the correct model order we will see no increase in model quality using the independent validation data set. If the true system however is in nite-dimensional, then the cross validation techniques will not give the same guidance since high order models will always approximate the underlying in nite-dimensional system better than low complexity models. However utilizing the cross validation step will tell us if the data are the noisy measurements of a nite-dimensional system or from a noise free in nite-dimensional system. In the frequency-domain, the cross validation is easily performed by dividing the frequency response measurements in two disjoint sets estimation data and validation data. The most natural division is to take every other frequency point as the estimation data and the rest as validation data. This division also preserves the uniform spacing between the frequencies. A model is then estimated using the estimation data only and the quality of the model is assessed by comparing the estimated transfer function with the validation data set.. 4.3 Model quality assessment. The following two norms on the error between the predicted and measured frequency responses will be used to compare dierent models ^ j!k (i) kG ; G^ km1 := max !k jGk ; G(e )j 13.

(14) Estimation error on estimation data and validation data 8. 7. 6. Error. 5. 4. 3. 2. 1. 0 15. 20. 25. 30 Model order n. 35. 40. 45. Figure 4.1: jjG ; G^ jj2e error on estimation data \+" and validation data \*" for dierent model orders using algorithm (3.1{3.19). (ii) kG ; G^ km2 :=. r1. ^ j!k 2 N jGk ; G(e )j :. 4.4 Model Order Determination. To determine an appropriate model order we divided the JPL-data into two disjoint sets as described previously. In Figure 4.1, the results of estimating models of order 19 to 43 is given using Algorithm 3.1 together with the estimation data and validate using the validation data. From the gure we notice that the transfer function error measured on the validation data keeps decreasing as model order increases. This gives an indication that the JPL-data is fairly noise free and originate from a system of a high order, possibly in nite. Any choice of model order will thus be a trade o between model complexity and amount of unmodeled dynamics we are willing to accept. We will estimate models in a range of model orders to illustrate the trade o.. 4.5 Estimation Results. To obtain best possible models we use all frequency response data in the identi cations. Figure 4.2 shows the magnitudes of the frequency data together with an estimated model of order 24 using Algorithm 3.1. The estimated model accurately picks up all the eleven dominant peaks. Increasing model order further improves the t as expected. In Figure 4.3, jjG ; G^ jjm1 is plotted for the models of order 24 to 62 using Algorithm 3.1 presented in the paper and the nonlinear least-squares method. The jjG ; G^ jjm2 errors have also the same characteristics. For model order 24, our algorithm produced jjG ; G^ jjm1 = 13:2 and for model order 42 we obtained jjG ; G^ jjm1 = 2:3. It is quite clear that, for this data, the algorithm outperforms the least-squares 14.

(15) Measured and estimated transfer function, model order 24 120. 100. Magnitude. 80. 60. 40. 20. 0 0. 100. 200. 300 400 Frequency rad/s. 500. 600. 700. Figure 4.2: Magnitude of measured transfer function (solid line) and estimated model of order 24 using Algorithm 3.1 (dashed line). All frequency response data is used in the identi cation and shown in the graph. method for all model orders. This data have also been used by Gu and Khargonekar in 11]. They presented a model of order 24 with jjG ; G^ jjm1 = 13. This model was obtained by SK iterations, see 29]. By introducing a second step and increasing model order to 42, they reduced the error to jjG ; G^ jjm1 = 6:1. Comparing these results with our results indicates that the presented algorithm is a promising alternative to existing methods. The algorithm proposed in 18] produced almost identical results as Algorithm 3.1, which is not surprising since, as we previously showed, for uniformly spaced frequencies the two methods are closely related and the noise level is low. For model orders larger than 60, the nonlinear least-squares method produced poor models probably due to the numerical problems with the coecients of the high order polynomials. It is also interesting to mention that model reductions of the 60th order model to a 24th order model always produced a model with higher errors than the directly estimated model by the subspace technique. This fact holds for truncated balanced realization, Hankel norm reduction and model reduction by calculating the Markov parameters and estimating a lower order model by Kung's algorithm. The latter result is surprising since this is part of the approach suggested by Bayard 4].. 5 Robustness and Sensitivity In the previous section, we saw that an nth order system can be recovered by a simple algorithm utilizing only the DFT, SVD, and a least-squares solution, given n + 2 noiseless frequency response measurements of the system. The model order is not known a priori and has to be inferred from the data. We also admit the presence of unmodeled dynamics. To be speci c, suppose that we are 15.

(16) Error for NNLS estimates and subspace estimates 40. 35. 30. Infinity error. 25. 20. 15. 10. 5. 0 20. 25. 30. 35. 40 45 Model order n. 50. 55. 60. 65. Figure 4.3: jjG ; G^ jj1e error for dierent model orders. \*" Algorithm 3.1, \+" Nonlinear leastsquares method. All data is used in the estimations. given the frequency response data Gk = G(ejk=M ) + ek k = 0 1 : : : M (5.1) where G represents the transfer function of a stable single-input/single-output (SISO) system of nth order and e captures unmodeled dynamics, assumed only to be norm bounded kek1 := sup jek j

(17) : (5.2) k. In the model structure above, our lack of information or our will to neglect the detailed features of the system such as small time-variations, nonlinearities etc. is reected by a small, norm bounded quantity e. In this section, we show that Algorithm 3.1 is robust to inaccuracies in the data and provide a criterion to estimate a model order when the data is corrupted as in (5.2). Let Hqr denote the Hankel matrix in (3.13) for e = 0 and de ne H~ qr := H^ qr ; Hqr . Let Hqr = U1 1V1T be its SVD as in (3.15). The next lemma provides perturbation bounds on 1 caused by e.. Lemma 5.1 Let frequency response measurements of G be corrupted as in (5.1). Then 2 (H~ ) 2 + log(q + r ; 1)]

(18) + (q + r ; 1)

(19) : 1. qr. 2M. Proof: See Appendix B.1.. (5.3). State-space parameters for the identi ed model are calculated from a set of basis vectors of the estimated observability range space. Hence we need to bound small perturbations of the observability range space caused by unmodeled dynamics. 16.

(20) Lemma 5.2 Let frequency response measurements of G be corrupted as in (5.1). Assume that q r > n and h i p (5.4) 1 > 2 2 n;2 (Hqr ) 2 kgk1 + 1 (H~ qr ) kH~ qr kF := : Then there exists a unique matrix P satisfying. kP kF . (5.5). such that the range space of the columns of U1 + U2 P is the same as the range space of the columns of U^1 .. Proof: See Appendix B.2.. The following theorem establishes robustness of Algorithm 3.1 against small inaccuracies in the Hankel matrix Hqr .. Theorem 5.3 Let frequency response measurements of G be corrupted as in (5.1). Consider Al-. gorithm 3.1 and assume that q r > n for xed q and r. Let the order of G^ be n. Then lim sup kG^ ; Gk1 = 0: !0 kH~ qr kF . (5.6). Proof: See Appendix B.3.. Information on the system order is not a restriction at all since model orders can be increased up to r where r < q while maintaining exactness of the algorithm for the systems in the model set. Therefore q and pr may be xed a priori to certain values and suitable model orders are sought. Since kH~ qr kF q 1 (H~ qr ), we have the following robustness result for xed q and r lim sup kG^ ; Gk1 = 0:. (5.7). !0 kek1 . The following theorem furnishes a simple criterion to determine the model order when estimates of kGk1 (A) and

(21) are available. Recall that n ; 1 th order best approximation of an nth order system in the H1 -norm results in the error ;n 10]. Thus ;n indicates order robustness of the system. Theorem 5.4 Let frequency response measurements of G be corrupted as in (5.1). Let H^ qr denote q by r Hankel matrix in (3.13). Consider the SVDs of H^ qr in (3.15). Then. r ; 1)

(22) jn(H^ qr ) ; ;n(G)j 2 + log(q + r ; 1)]

(23) + (q +2M (5.8) h i +(1 ; );1 kGk1 q+r + 1+minfqrg + 2M +1 (1 ; 2M );1 : 2. Proof: See Appendix B.4.. The logarithmic rate in Theorem 5.4 suggests that q and r can be taken quite large to obtain accurate estimates of ;n. The conclusions in this section apply to the multiple-input multipleoutput systems and only Lemma 5.1 and Theorem 5.4 need trivial modi cations. 17.

(24) 6 Asymptotic Error Analysis In this section, we will investigate asymptotic properties of the subspace identi cation algorithm when the frequency response measurements are corrupted by noise. Under certain weak assumptions on the noise and the system, we will show that the transfer function estimates are consistent. We will assume that in the following frequency response data G~ k = G(ej!k ) + Nk k = 0 M (6.1) the disturbances are independent zero-mean complex random vectors. Each complex variable is assumed to have independent real and imaginary parts with equal distribution. Complex Gaussian noise (see 6] for a de nition of complex noise) is a special case of this general noise model. For a justi cation on the noise model chosen here, we refer the interested reader to 28, 30]. The following theorem demonstrates that the algorithm proposed by Liu et. al. 18] is consistent if the measurement noise is temporally and spatially white at output channels and has equal variances for all frequencies.. Theorem 6.1 Let the true system be described by the state-space equations (2.1). Suppose that the noise is zero mean, has bounded fourth moments and covariance of the form. E Nk NkH = 2 Ip 8k and the data is uniformly spaced in frequencies. Consider the subspace method 18]. Then G^ ! G w.p. 1 as M ! 1:. (6.2) (6.3). Proof: See Appendix C.1.. The uniform covariance requirement for all frequencies on the noise as well as temporal and spatial whiteness limits the practical use of the algorithm 18] to single output systems. However for the algorithm we proposed, we show that no such additional condition on the noise is required:. Theorem 6.2. Consider Algorithm 3.1 and assume that measurement noise is zero mean and has bounded second moments. Then G^ ! G w.p. 1 as M ! 1: (6.4). Proof: See Appendix C.2. This result is quite powerful since consistency is preserved even if noise is correlated between dierent outputs and/or the covariance is nonuniformly distributed over frequencies. The latter case happens if the measurement noise is colored in the time domain.. 7 Conclusions In this paper, we have developed a simple state-space identi cation algorithm to identify systems from equally spaced frequency response measurements. The algorithm was shown to be robust under small norm bounded perturbations. An asymptotic stochastic analysis showed that the algorithm is consistent under the weak conditions on the measurement noise. The algorithm was used to identify 18.

(25) a exible structure in state-space form and a comparison with a nonlinear least-square iterative algorithm was made. The results show that the algorithm outperforms the nonlinear least-squares algorithm (invfreqz) on these almost noise less data and is therefore a viable alternative to classical iterative methods.. Acknowledgments The authors are indebted to Dr. D.S. Bayard at the Jet Propulsion Laboratory, Pasadena, California who provided the experimental data used in this paper and motivated this work. The second author would also like to thank Prof. P.P. Khargonekar of the University of Michigan for introducing him to this problem during his doctoral studies in the University of Michigan.. References. 1] V.M. Adamyan, D.Z. Arov, and M.G. Krein, \Analytic Properties of Schmidt pairs for a Hankel operator and the generalized Schur-Takagi problem," Math. USSR Sbornik, Vol. 15, pp. 31{73, 1971.. 2] J.L. Adcock, \Curve Fitter for Pole-Zero Analysis," Hewlett-Packard Journal, pp. 33{36, 1987.. 3] H.T. Banks and Y. Wang, \Damping Modeling in Timoshenko Beams," Proceedings of the American Control Conference, Chicago, pp. 2139{2143, 1982.. 4] D.S. Bayard, \An Algorithm for State-Space Frequency Domain Identi cation Without Windowing Distortions," Proceedings of the 31st IEEE Conference on Decision and Control, Tucson, Arizona, pp. 1707{1712, 1992.. 5] D.S. Bayard, \Statistical Plant Set Estimation Using Schroeder-Phased Multisinusoidal Input Design," Proceedings of the American Control Conference, Chicago, pp. 2988{2995, 1992.. 6] D.R. Brillinger. Time Series: Data Analysis and Theory. McGraw-Hill, New York, USA, 1981.. 7] E.W. Cheney. Introduction to Approximation Theory. McGraw-Hill Inc., 1966.. 8] K.L. Chung. A Course in Probability Theory, Academic Press, San Diego, 1974.. 9] R.L. Dailey and M.S. Lukich, \MIMO Transfer Function Curve Fitting Using Chebyshev Polynomials," Proceedings of the SIAM 35th Anniversary Meeting, Denver, Colorado, 1987.. 10] K. L. Glover, \All Optimal Hankel-norm approximations of linear multivariable systems and their L1 error bounds," International Journal of Control, Vol. 39, pp. 1115{1193, 1984.. 11] G. Gu and P.P. Khargonekar, \ Frequency Domain Identi cation of Lightly Damped Systems: The JPL Example," Proceedings of the American Control Conference, pp. 3052-3056, 1993.. 12] B.L. Ho and R.E. Kalman, \Eective Construction of Linear State-Variable Models from Input/Output Functions," Regelungstechnik,, Vol. 14, pp. 545{592, 1966. 19.

(26) 13] J.N. Juang and R.S. Pappa, \An Eigensystem Realization Algorithm for Modal Parameter Identi cation and Model Reduction," Journal of Guidance, Control and Dynamics, Vol. 8, pp. 620{687, 1985.. 14] J.N. Juang and H. Suzuki, \An Eigensystem Realization Algorithm in Frequency Domain for Modal Parameter Identi cation," Journal of Vibration, Acoustics, Stress, and Reliability in Design, Vol. 110, pp. 24{29, 1988.. 15] S. Y. Kung, \A New Identi cation and Model Reduction Algorithm via Singular Value Decomposition," Proceedings of the 12th Asilomar Conference on Circuits, Systems and Computers, Paci c Grove, California, pp. 705{714, 1978.. 16] E.C. Levy, \Complex Curve Fitting," IRE Transactions on Automatic Control, AC 4, pp. 37{44, 1959.. 17] J. Little and L. Shure. Signal Processing Toolbox. The Mathworks, Inc., 1988.. 18] K. Liu, R.N. Jacques, and D.W. Miller. Frequency Domain Structural System Identi cation by Observability Range Space Extraction. Technical Report. Space Engineering Research Center, MIT, Cambridge, Massachusetts, 1993.. 19] K. Liu, R.N. Jacques, and D.W. Miller. Frequency Domain Structural System Identi cation by Observability Range Space Extraction. Proc. of the American Control Conference, Baltimore, Maryland, vol.1, pp 107{111, 1994.. 20] K. Liu and R.E. Skelton, \Q-Markov Covariance Equivalent Realization and its Application to Flexible Structure Identi cation," AIAA Journal of Guidance, Control and Dynamics, pp. 308{319, 1993.. 21] L. Ljung. System Identi cation: Theory for the User. Prentice-Hall, Englewood Clis, New Jersey, 1987.. 22] L. Ljung, \Some Results on Identifying Linear Systems Using Frequency Domain Data", Proceedings of the Conference on Decision and Control, pp. 3534-3538, 1993.. 23] M. Moonen, B.D. Moor, L. Vandenberghe, and J. Vandewalle, \On- and o-line identi cation of linear state-space models," International journal of Control, pp. 219{232, 1989.. 24] B.D. Moor, \The Singular Value Decomposition and Long and Short Spaces of Noisy Matrices," IEEE Transactions on Signal Processing, February 1994.. 25] L.E. Pfeer, \The RPM Toolbox: A System for Fitting Linear Models to Frequency Response Data," Proceedings of the 1993 MATLAB Conference, Cambridge, Massachusetts, October 1993.. 26] B.D. Moor, M. Moonen, P.V. Overschee, \State Space Subspace System Identi cation (S4ID)," Automatica, 1994.. 27] A. Pinkus. n-Widths in Approximation Theory. Springer-Verlag, Berlin, 1985. 20.

(27) 28] R. Pintelon, P. Guillaume, Y. Rolain, J. Schoukens, and H.V. Hamme, \Parametric Identi cation of Transfer Functions in the Frequency Domain, a Survey," Proceedings of the Conference on Decision and Control, pp. 557-566, 1993.. 29] C.K. Sanathanan and J. Koerner, \Transfer Function Synthesis as a Ratio of Two Complex Polynomials," IEEE Transactions on Automatic Control, pp. 56{58, January 1963.. 30] J. Schoukens, R. Pintelon. Identi cation of Linear Systems: a Practical Guideline to Accurate Modeling. Pergamon Press, London, U.K., 1991.. 31] T. S&oderstr&om and P. Stoica. System Identi cation. Prentice-Hall International Series in Systems and Control Engineering, 1989.. 32] G.W. Stewart and J.G. Sun. Matrix Perturbation Theory. Academic Press, Inc., 1990.. 33] P. Stoica, P. Eykho, P. Janssen, and T. S&oderstr&om. \Model structure selection by cross validation", Int. J. Control, Vol. 43, pp. 1841{1878, 1986.. 34] M. Verhaegen, \A Novel Non-Iterative MIMO State-Space Model Identi cation Technique," Proceedings of the 9th IFAC/IFORS Symposium on Identi cation and System Parameter Estimation, Budapest, Hungary, pp. 1453{1458, 1991.. 35] M. Viberg, B. Ottersten, B. Wahlberg, and L. Ljung, \A Statistical Perspective on State-Space Modeling Using Subspace Methods," Proceedings of the 30th IEEE Conference on Decision and Control, Brighton, England, pp. 1337-1342, 1991.. 36] A. H. Whit eld, \Asymptotic Behaviour of Transfer Function Synthesis Methods," International Journal of Control, Vol. 43, pp. 1413{1426, 1987.. 37] N.J. Young. An Introduction to Hilbert Space. Cambridge University Press, 1988.. A Proof of Theorem 3.2 Since G is a stable transfer function, it can be represented by the following Taylor series. G(z ) = D + C (zI ; A);1B = D +. 1 X. k=1. CAk;1Bz;k =. 1 X. k=0. gk z ;k < jzj. (A.1). where = (A) < 1 since G is a stable system. Notice that g^k can be written as. g^k =. 1 X. i=0. g(k + 2iM ) = CAk;1. X 1 i=0. A. 2iM. !. B = CAk;1(I ; A2M );1 B. (A.2). and therefore H^ qr can be factored as. H^ qr = Oq (I ; A2M );1 Cr (A.3) From the dimensions of the factors in (A.3) it is clear that H^ qr has a maximal rank n. Furthermore since < 1, (I ; A2M ) is always of rank n. Minimality of the system also implies that both Cr and 21.

(28) Oq are of rank n and hence also H^ qr if r n and q > n. In (3.15) then ^ 2 = 0 and the column range spaces of H^ qr , Oq and U^1 will be equal. A valid extended observability matrix O^q is then given by U^1 since there exists a non-singular n n matrix S such that O^q = U^1 = Oq S (A.4) ^ B ^ C ^ D^ ) which is U^1 is thus an extended observability matrix from a state-space realization (A ^ ^ similar to the original realization (A B C D). (3.16) and (3.17) will yield A and C matrices which are related to the original realization as (A.5) A^ = S ;1 AS C^ = CS as shown in 15]. Since (3.18) has a unique solution, we get B^ = S ;1 B , D^ = D: Then (A B C D) ^ B ^ C ^ D^ ) are similar. Letting q = n + 1 r = n M = n + 2 , we satisfy the condition and (A maxfq rg < M .. B Proofs in Section 5 B.1 Proof of Lemma 5.1. Without loss of generality assume G = 0. Let G^ denote the FIR system with impulse response coecients g^k k = 1 2 : : : n0 n0 X G^ (z ) = g^k z;k (B.1) k=1. where n0 = q + r ; 1. We derive after some algebra. q+r sin n20 ( k 1 2MX;1 ej ( k M ; ) Gk M ; ) 2 (B.2) G^ (ej ) = 2M k sin( M ; ) 21 k=0 where frequency response data for k > M are calculated from (3.1). Hence n0 k n0 k MX ;1 MX ;1 ( ;. ) ( ;. ; ) sin sin

(29)

(30) j jG^ (e )j 2M (B.3) sin(2kM; ) 1 + 2M sin(2kM; ; ) 1 : M 2 M 2 k=0 k=0 Notice that the right hand side of (B.3) approximates the following integral # "

(31) Z ; sin n20 t dt + Z ; sin n20t dt =

(32) Z sin n20t (2 + log n )

(33) : (B.4) 0 t 2 ;; sin 2t 2 ; sin 2t ; sin 2 The last inequality above is well known in the Fourier series literature (see, for example 7]). Since the derivatives of the integrands are bounded by n20 =2, the approximation error is less than n20

(34) =2M . Now consider the following Hankel matrix. 0 BB H^ n0n0 = B B@. g^1 g^2 .. . g^n0. g^2 g^3 ... 0. 1 g^n0 C ! 0 C ^ H A qr 1 . . . ... C CA := A2 0 : 0 22. (B.5).

(35) Since H^ qr is a sub-block of H^ n0 n0 , we have 1(H^ qr ) 1 (H^ n0n0 ): (B.6) 0 ;1 ^ ^ ^ ^ Notice that Hn0 n0 0is the Hankel matrix of G (z ) := 0G(z ) and therefore 1 (Hn0 n0 ) is equal to the Hankel norm of G^ which is bounded above by kG^ k1 by the Nehari's Theorem 37]. Thus from (B.3{B.4), and (B.6), (5.3) follows.. B.2 Proof of Lemma 5.2 Let. W^ qq := H^ qr H^ qrT Wqq := Hqr HqrT (B.7) and consider the following singular value decompositions of W^ qq and Wqq W^ qq = U^1 ^ 21 U^1T + U^2 ^ 22 U^2T Wqq = U1 21 U1T + U2 22 U2T (B.8) where 2 = 0. De ne W~ qq as (B.9) W~ qq := W^ qq ; Wqq = Hqr H~ qrT + H~ qr HqrT + H~ qr H~ qrT and let ~ ~ ! F~ := U T W~ qq U := FF~11T FF~12 (B.10) 22 12 where U := U1 U2 ]: It is known ( 32], Theorems 2.8 and 3.1) that if a Hermitian matrix is perturbed by a Hermitian matrix in the Frobenius norm then the column range spaces of U^1 and U1 + U2 P are equal for some unique matrix P satisfying 2 kF~12 kF kP kF 2 (B.11) n(Hqr ) ; kF~11 kF ; kF~22 kF provided that n2 (Hqr ) > kF~11 kF + kF~22 kF + 2 kF~12 kF : (B.12) Since the Frobenius norm is unitarily invariant, we can derive the following upper bounds p p kF~11 kF + kF~22 kF 2 kF~qq kF = 2 kW~ qq kF (B.13) p 2 kF~12 kF 2 p1 kF~qq kF = 2 kW~ qq kF : (B.14) 2 If (B.12) is replaced by the following stronger condition p (B.15) 1 > 2 2 n;2 (Hqr ) kW~ qq kF : we can bound kP kF as p kP kF 2 2 n;2 (Hqr ) kW~ qq kF : (B.16) Now the following chain of inequalities complete the proof. kW~ qq kF 2 kHqr H~ qrT kF + kH~ qr H~ qrT kF (B.17). X 1 2 1kmax g(k + 2iM ) kH~ qr kF + 1(H~ qr ) kH~ qr kF q+r;1 i=0 h i 2 kgk1 + 1(H~ qr ) kH~ qr kF : 23.

(36) B.3 Proof of Theorem 5.3. From Lemma 5.2, we can write U^1 as U^1 = (U1 + U2 P ) Q := U~1 Q. (B.18). for some invertible matrix Q and a unique matrix P . First, we show that due to the least-squares method used in the estimation of B^ and D^ , G^ is invariant to Q and ^ 1 . Suppose that A^ and C^ are calculated from (3.16-3.17). Let " ^# B E^ := D (B.19) ^ Without loss of generality, we may assume m = 1 since (3.18) is a columnwise separable leastsquares problem. Let IR and II denote respectively, projections onto reals and imaginaries. Let. . ;1. X^ k := C^ ejk=M In ; A^. . Ip k = 0 M. 2 ^ 3 X1 ] 7 66 IR. II. X 66 .^1 ] 777 X^ := 66 .. 77 G^ := 64 IR X^2 ] 75 ;II X^ 2]. 2 IR G ] 3 66 II G11] 77 66 .. 77 : 66 . 77 4 IR G2] 5 ;II G2]. (B.20). (B.21). Note that in the above, (3.1) was used. Then E^ is calculated from (3.18) as. E^ = X^ yG^:. (B.22). ^ C^ calculated for U~1 by taking ^ 1 = In : Let E^0 be the solution (B.22) for Let A^0 C^0 denote A ^ ^ A0 C0. Then it is easily shown A^ = ^ 1;1=2 Q;1 A^0 Q ^ 11=2 C^ = C^0 Q ^ 11=2 (B.23). ". Hence we get. ^ ;1=2 ;1 E^ = 1 0 Q. #. 0 E^0 : Ip. (B.24). ;1 ;1 G^ (z) = C^ zIn ; A^ B^ + D^ = C^0 zIn ; A^0 B^0 + D^ 0 (B.25) where B^0 and D^ 0 partition E^0 as in (B.19). Since q and r are xed U1 and U2 are xed. Thus A^0 and C^0 and consequently B^0 and D^ 0 continuously change by P . Hence from (B.25), G^ changes continuously by P and by Lemma 5.2, P continuously change by Hqr . Therefore G^ is a continuous function of Hqr at G.. 24.

(37) B.4 Proof of Theorem 5.4 Assume rst that e = 0 and let. Hqr := Mlim H^ qr : (B.26) !1 First, we will bound perturbation of n (Hqr ) due to the DFT. Notice that g^k can be written as g^k =. 1 X. i=0. g(k + 2iM ) k = 0 1 2M ; 1:. (B.27). ;1 Let X (z ) = Pki=0 x(i)zi be a k ; 1th order best approximation of G0 (z) := G(z;1 ) in the H1 norm. Then it is known 27] kG0 ; X k1 kG0 k1 k (A): (B.28) 0 0 Since x(i) = 0 i k and kgk1 kG k1 for any transfer function G 2 H1 , it follows that. jgk j kG0 k1 k :. (B.29). Thus from (B.27) and (B.29), we get. jg^k ; gk j . 1 X. i=1. kG0 k1k+2iM = 2M +k (1 ; 2M );1 kG0 k1. (B.30). and the perturbation of the smallest singular value of Hqr can be bounded from. jn(H^ qr ) ; n(Hqr )j 1 (H^ qr ; Hqr ) . max Xfqrg k=1. jg^k ; gk j. (B.31). where the rst inequality above follows0 from a result of Mirsky ( 32], Theorem 4.11). Now let Gq+r be a truncated impulse response of G as follows q+X r;1. Gq+r :=. k=1. Notice that ;1 (G0 ; Gq+r ) kG0 ; Gq+r k1 . gk zk :. (B.32). jgk j (1 ; );1 kG0 k1 q+r . (B.33). 1 X. k=q+r. where the rst inequality came from the Nehari's Theorem 37]. Consider the following Hankel matrix. 0 BB gg12 Sq+r := B B@ .... gq+r;1. g2 g3 .. . 0. 1 gq+r;1 C ! 0 C H A qr 1 := A 0 : .. C ... 2 . C A 0 25. (B.34).

(38) Since Gq+r is an FIR system, Sq+r represents the Hankel matrix of Gq+r . Hence n(Sq+r ) = ;n (Gq+r ) and we get " # q+X r;1 0 A 1 jn(Sq+r ) ; n(Hqr )j 1 ( A 0 ) jgk j 2. (1 ; );1. minfqrg kG0 k1 1+minfqrg):. (B.35) (B.36). 1+. An application of Adamyan-Arov-Krein theorem 1] gives 0 ;n(G ) ; ;n(Gq+r ) ;1(G0 ; Gq+r ):. (B.37). Hence from Lemma 5.1, (B.30{B.31), (B.33), and (B.35{B.37), we obtain (5.8).. C Proofs in Section 6 C.1 Proof of Theorem 6.1. Without loss of generality, we may assume that an even number of data are uniformly spaced on. 0 2]. We introduce a noise matrix 2 3 N1 0 0 66 .. 77 0 N . 77 6 2 ::: NqM = (WqM Ip) 66 . . . (C.1) .. 7 . . . . . . . 4 5 0 NM and expand it with its complex conjugate N = NqM NqM ]: (C.2) Since we are given frequency response samples Y~ (!k ) = G(ej!k ) + Nk k = 1 M directly, we let UqM = WqM . Recall equation (2.27) which is used in 18] to estimate the observability range space. In the analysis, we introduce the \sample covariance matrix" of size qp qp (C.3) G^ := M1 Y~ U ?Y~H which is bounded as M tends to in nity but has the same range space as (2.27). The equation (C.3) will then be G^ = M1 Oq X U ? X H OqT + M1 N U ?N H + M1 Oq X U ? N H + M1 N U ?X H OqT (C.4) Our aim is now to show that the range space of G^ will converge to the range space of Oq with probability one. In (C.4) we recognize one deterministic term and three stochastic terms. For the fourth term in (C.4), we notice that the elements in this qp qp matrix can be described by 1X 1 2M ? H T Sij (M ) := M N U X Oq = M Nk cij (k) ij k=1 26.

(39) for some bounded constants cij (k). All these elements are zero mean and have a fourth moment satisfying E jSij (M )j4 Mc13 + Mc22 (C.5) for some constants ci . By Chebyshev's inequality we have for each

(40) > 0 4 P (jSij (M )j >

(41) ) = P (jSij (M )j2 >

(42) 2 ) E jSij

(43) (4M )j :. Thus we obtain from (C.5-C.6). 1 X. M =1. (C.6). P (jSij (M )j >

(44) ) < 1. and it follows from Borel-Cantelli's lemma 8] that. Sij (M ) ! 0 w.p. 1 as M ! 1: which of course holds also for the third term in (C.4). Since the frequencies !k k = 1 M are uniformly distributed between 0 and 2 we have 1 UHU: (C.7) U ? = I2M ; 2M The second term in (C.4) naturally divides into two terms 1 ? H 1 1 H H H (C.8) M NU N = M NN + 2M 2 NU UN The elements of the second matrix in (C.8) is of the form. 1. M. NU H UU H UN H. . ij. 2M X 2M X = M1 2 cij (k l)Nk NlH. k=1 l=1. for some bounded constants cij (k l). The second moment of this term is then bounded as. 2 c3 c4 1 H H H E M NU UU UN M 3 + M 2 ij. for some bounded constants ci . Again appealing to Borel-Cantelli lemma, we conclude that this term also converges to zero w.p. 1. The rst term in (C.8) is simply. " i. 1 N N H = 1 h W H H qM WqM Ip diag(N1 N1 : : : N2M N2M ) M M By the assumption (6.2) and since. i " 1 h W qM WqM Ip M 27. #. !. #. H WqM Ip T WqM. H WqM Ip = 2Iqp T WqM. !. (C.9). (C.10).

(45) the expected value of NN H =M is equal to the scaled identity matrix 22 Iqp. Consider 1 NN H ; 22 I qp M which can also be written as " H # ! i. 1 h W WqM I (C.11) H 2 H 2 qM WqM Ip diag(N1 N1 ; Ip : : : N2M N2M ; Ip ) p T WqM M which now is zero mean. It is now straight forward to show that the fourth moments of the elements in the matrix (C.11) are bounded by. 4 c5 c6 1 H 2 E M NN ; 2 Iqp M 3 + M 2 ij for some constants ci and we conclude that M1 NN H ! 22 Iqp w.p. 1 as M ! 1. Finally, we arrive at. 1 O X U ? X H OT + 22 I w.p. 1 G := Mlim G^ = Mlim q qp q !1 !1 M. Assuming the signal part has the following (asymptotic) SVD. ". 1 O X U ? X H OT = h U U i 1 0 lim 1 2 q 0 0 M !1 M q where 1 is a diagonal matrix with n positive elements, we obtain. i". h. #". #". U1T U2T. #. (C.12) (C.13). #. U1T : G = U1 U2 (C.14) U2T Thus we have shown that the n left singular vectors corresponding to the n largest singular values of G^ converge w.p. 1 to a set of basis vectors for the observability range space as M ! 1. Thus A and C , up to a similarity transformation, will be consistently estimated. From (B.25), note that G^ is invariant to the similarity transformations. Thus given A^ and C^ , it follows from (B.22) that B^ and D^ are consistent. 1. + 22 I 0 0 22 I. C.2 Proof of Theorem 6.2. Let Hqr denote the Hankel matrix in (3.13) for noise free data and de ne H~ qr := H^ qr ; Hqr : (C.15) We will show that kH~ qr kF ! 0 w.p. 1 as M ! 1: (C.16) Without loss of generality, we may assume that system is single-input/single-output. Let w = 2M and = 2(l + m ; 1): The l m th element of H~ qr is then. h~ i Hqr. ;1 1 wX j k=w := 1 Sw : = N e k lm w k=0 w. 28. (C.17).

(46) Let R overbound noise covariances. We have by (C.17) E jSw j2 Rw: (C.18) It follows by Chebyshev inequality that for each

(47) > 0 we have R : (C.19) P (jSw j > w

(48) ) w

(49) 2 Taking a subsequence Sw2 , we get X X P (jSw2 j > w2

(50) ) = wR2

(51) 2 < 1: (C.20) w w Hence by Borel-Cantelli lemma, we have Sw2 ! 0 w.p. 1 as w ! 1: (C.21) w2 We have thus proved the desired result for a subsequence. Now it will be extended to the whole sequence. For w2 < k < (w + 1)2 from. Sk ; Sw2 = we have. 2 ;1 wX. i=0. h j i=k. Ni e. ; ej i=w2. ;1 i kX. +. i=w2. Ni ej =k . w2;1 2 k;1 2 h i X X E jSk ; Sw2 j2 2 E Ni ej i=k ; ej i=w2 + 2 E Ni ej =k 2 i=0. 2 sin w + 4R w 4R (1 + ) w: 4Rw k! 2. Let Then we have. Dw :=. max jS ; Sw2 j : w2 k<(w+1)2 k. i=w. 2(w+1)2 ;1 3 X E Dw E 4 jSk ; Sw2 j2 5 4R (1 + ) w2 h 2i. k w. = 2. (C.22) (C.23). (C.24) (C.25). and consequently by Chebyshev inequalty It follows as before Since for w2 < k < (w + 1)2. P (Dw > w2

(52) ) 4Rw(12

(53) +2 ) :. (C.26). Dw ! 0 w.p. 1 as w ! 1: w2. (C.27). jSw j jSw2 j + Dw . (C.28) w w2 it is clear that (C.21) and (C.27) together imply h~ i Hqr lm ! 0 w.p. 1 as M ! 1: (C.29) Thus kH~ qr kF converges to zero w.p. 1 since q and r are nite. Theorem 5.3 completes the proof since it shows that G^ : IRqr ! H1 is a continuous map at G. 29.

(54)

No results found