Linear Continuous-time System Identification and State Observer Design by Modal Analysis by Mohamed Hassan El-Shafey B.Sc. Ain Shams University, Cairo, Egypt, 1978 M.Sc. Ain Shams University, Cairo, Egypt, 1982 A Thesis Submitted in Partial Fulfilment of The Requirements for The Degree of Doctor of Philosophy in The Faculty of Graduate Studies (Department of Electrical Engineering) We accept this thesis as conforming to the required standard The University of British Columbia November, 1987 Â© Mohamed Hassan El-Shafey, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Â£~lec tcfCcj I ]~y\cj,'A eâ‚¬fjf\ <j The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date Feb. 5, /fSfr Abstract ii A new approach to the identification problem of linear continuous-time time-invariant systems from input-output measurements is presented. Both parametric and nonparametric system models are considered. The new approach is based on the use of continuous-time functions, the modal functions, defined in terms of the system output, the output derivatives and the state variables under the assumption that the order n of the observable system is known a priori. The modal functions are obtained by linear filtering operations of the system output, the output deriva tives and the state variables so that the modal functions are independent of the system instantaneous state. In this case, the modal functions are linear functions of the input exponential modes, and they contain none of the system exponential modes unlike the system general response which contains modes from both the sys tem and the input. The filters parameters, the modal parameters are estimated using linear regression techniques. The modal functions and the modal parameters of the output and its deriva tives are used to identify parametric input-output and state models of the system. The coefficients of the system characteristic polynomial are obtained by solving n algebraic equations formed from the estimates of the modal parameters. Esti mates of the parameters associated with the system zeros are obtained by solving another set of linear algebraic equation. The system frequency response and step response are estimated using the output modal function. The impulse response is obtained by filtering the estimated step response using the output first derivative modal parameters. A new method is presented to obtain the system poles as the eigenvalues of a data matrix formed from the system free response. The coefficients of the system characteristic polynomial are obtained from the data matrix through a simple re cursive equation. This method has some important advantages over the well known Prony's method. The state modal functions are used to obtain a minimum-time observer that gives the continuous-time system state as a direct function of input-output samples in n sampling intervals. Contents Abstract ii Table of Contents iiList of Tables vi Acknowledgements vii1 Introduction 1 1.1 Introduction1.2 Parametric Identification 3 1.2.1 Continuous Models from Discrete Equivalents 3 1.2.2 Modulating Functions 6 1.2.3 Multiple Integrations and Orthogonal Functions 7 1.3 Nonparametric Identification 8 1.3.1 Frequency Response 9 1.3.2 Impulse Response 10 1.4 Scope of the Thesis 1 2 Definitions and Properties of Modal Functions 13 2.1 Introduction 12.2 Mathematical Background 13 2.3 Definitions of Modal Functions and Modal Parameters 16 2.4 Properties of Modal Functions 19 iii CONTENTS iv 2.5 Summary 23 3 Estimation of The Modal Parameters 24 3.1 Introduction 23.2 Problem Formulation 4 3.2.1 Estimation from Discrete Measurements 25 3.2.2 Estimation from Continuous Measurements 26 3.3 Estimating The Output Modal Parameters 27 3.3.1 Estimates from Discrete Measurements 28 3.3.2 Estimates from Continuous Measurements 29 3.4 Estimating The Output Derivative Modal Parameters 33 3.4.1 Estimates Using Fourier Functions 34 3.4.2 Estimation Using Walsh Functions 6 3.5 Estimation Using Multi-mode Inputs 38 3.6 Summary 49 4 Parametric Identification 50 4.1 Introduction 54.2 Identification Using the Output Modal Function 50 4.3 Identification Using Modal Functions and Parameters 57 4.3.1 Estimation of a 54.3.2 Estimation of b 8 4.4 Estimate of a from the Free Response 61 4.5 Summary 75 Nonparametric Identification 73 5.1 Introduction 75.2 Estimation of the Frequency Response 74 CONTENTS v 5.3 Estimation of the Step Response 78 5.4 Estimation of the Impulse Response 82 5.5 Summary 85 6 Observer Design 7 6.1 Introduction 86.2 State Modal Functions as a State Observer 87 6.3 Estimation of I 94 6.4 Summary 6 7 Conclusion 98 7.1 Achievements7.2 Remarks on Future Research 100 Bibliography 102 List of Tables 2.1 Modal parameter vectors 23 3.1 Estimate of 6_0 using discrete data 29 3.2 Estimates of 90 using continuous data 30 3.3 Estimate of 0_o using discrete noisy data 2 3.4 Estimates of 0o using continuous noisy data 32 3.5 Estimate of 0_o using an instrumental variable 33 3.6 Estimates of Oj using Fourier functions 36 3.7 Estimates of Oj using Walsh functions 8 3.8 Estimates of 0* using Fourier functions 45 3.9 Estimates of 02 using Fourier functions3.10 True values of pj 46 3.11 Estimates of 0} using Walsh functions 43.12 Estimates of 02 using Walsh functions 6 3.13 Estimates of the modal parameters using Fourier functions 48 3.14 Estimates of the modal parameters using Walsh functions 48 3.15 True values of the modal parameters 44.1 Estimates of p0 and qâ€ž, k = 1,2,3 from discrete data 54 4.2 Estimate of the model parameters 55 4.3 Estimates of the model parameters from noisy samples 56 vi LIST OF TABLES vii 4.4 Estimates of the model parameters using instrumental variables. . . 56 4.5 Estimates of a and the system poles using Fourier functions and Walsh functions 58 4.6 Estimates of a and the system poles using Fourier and Walsh func tions 59 4.7 Estimates of b using Fourier functions and Walsh functions 61 4.8 Estimates of a and the system poles 66 4.9 Eigen values of Aa 68 4.10 Estimates of from the discrete model 71 4.11 Estimates of Sk using the new method5.1 Spectral components of u0 and y0 77 5.2 Estimated and true values of samples of the frequency response. . . 77 5.3 Samples of g(t); estimated and true values 83 5.4 Samples of h(t); estimated and true values 6 6.1 Minimum time observer parameters 95 6.2 Estimated and true values of b . 97 viii Aknowledgements I wish to thank Dr. E. V. Bohn for supervising my research work and for his fruitful suggestions. His financial support through NSERC grant is also appreciated. Chapter 1 Introduction 1.1 Introduction The analysis of a physical system starts with the formulation of a mathematical model of that system. Determining the system model may be possible in some cases by detailed analysis using the basic laws of system dynamics. Usually, this is not possible in practice because of the lack of sufficient knowledge of the system and it is therefore necessary to determine the model from observations made on the system. However, the available a priori knowledge of the system is often sufficient to determine the class of models to which a suitable model of the physical system belongs. In other words, to assume that whether a good model of the system is linear or nonlinear, finite dimensional or infinite dimensional, time-invariant or time varying, continuous-time or discrete-time, stochastic or deterministic, is often possible from the available a priori knowledge of the system. A good model should give a close approximation to the physical system according to certain chosen criterian besides being identifiable from the available measurements taken from the system input and output. The problem of identifying a linear lumped parameter time-invariant system is the simplest and most tractable. This is the main reason why linear time-invariant models have always received considerable interest in the control literature despite the fact that most systems in practical life do not fall in this category. Nonetheless, a linear time-invariant model can be made to serve as a useful approximation to a physically nonlinear system about a nominal operating point. Moreover, the parameters of a time-invariant model can be updated on-line to track possible 1 Introduction 2 variations in the system dynamics under the assumption that these variations occur relatively slowly compared to the parameter updating algorithm. Since the sixties, with the ready availability of digital computers, most of the research work in system identification was drawn towards discrete-time modeling by difference equations so that the system model matches the serial processing nature of a digital computer. However, most physical systems obey physical laws and have continuous-time outputs so that it is natural to model the system in con tinuous time by means of differential equations. The identification of a differential equation model of the system is consequently an important problem. This problem is inherently harder than the identification of a difference equation model because of the need to obtain derivatives which are usually unmeasurable. Modeling a single input-single output linear system by a differential equation of order n can in general take two forms: 1. Input-output model; Dny(t) + Za^y(t) = 5>,-P''u(*) (1-1) where u and y are the input and the output signals respectively and where the operator V is the time differentiation operator defined by dP 2. State model; (a) State equation; Dx[t) = Ax[t) + bu{t) (1.2) (b) Output equation; y{t) = cTx{t) (1.3) where x is the state vector, A is an nxn matrix1 and 6 and c are vectors of dimension n each. The superscript T means the transpose. Note that the elements of vector b in (1.2) are different from b0,..., 6â€ž_i in (1.1). The parametric system identification problem is to obtain estimates of the parameters a0,... ,an-i and 6â€ž_i of the system input-output model (1.1), or estimates of the arrays A, b and c of the 1 Sometimes called the system dynamic matrix Introduction 3 state model (1.2) and (1.3) from the available measurements made on the system. However, it is not possible to obtain all the elements of A, b and c from the input-output measurements only; some assumptions on the elements of these arrays have to be made which lead to the canonical forms which will be discussed in some detail later. Complementary to the identification problem of a finite dimensional parametric model is the identification of an infinite dimensional nonparametric model, e.g. the frequency response or the impulse response. In this case it is usually required to obtain samples of the frequency response or the impulse response. The following two sections review the important achievements in parametric and nonparametric system identification respectively. The last section in this in troductory chapter gives a overview of the thesis. 1.2 Parametric Identification A relatively recent survey of the methods of parametric identification of continuous-time systems can be found in [Young 81]. The survey classifies the parameter esti mation algorithms to output error algorithms and equation error algorithms accord ing to the cost function that is minimized. The methods reviewed in the following sections are basically of the equation error type which is adopted throughout the thesis. The term continuous systems will be used throughout the thesis to refer to continuous-time systems. 1.2.1 Continuous Models from Discrete Equivalents In order to make use of the advanced techniques of discrete model identification, the problem of identifying a continuous transfer function model2. V"-1 h V 2Equivalent to the input-output model (1.1) Introduction 4 is divided into two parts. The first part consists of identifying a discrete model of the system EA*'l) = iffe^ (1-5) from discrete measurements. The second part is the determination of H(s) from Hd{z~l) which is the inverse problem of digital filter design from the equivalent continuous filter [Oppenheim 75, Sec.5.1]. The discrete transfer function Hd{z~l) represents the ^-transform of H(s) only if the input is a sequence of impulses. Otherwise, and under the assumption that the measurements are discrete, the input has to be approximated between the sampling instants. Different approximations of the input signal, e.g. by steps or ramps, lead to different transformations, step invariant or ramp invariant, from H{s) to Hd{z~l) [Sinha,IM.K. 72]. One difficulty with this method lies in the need to factorize the denominator polynomial of H^z'1) in order to obtain the poles of of H(s) through a logarithmic operation Sk = -j,\nzk. (1.6) Other transformations which do not require polynomial factorization or logarithmic operations are given by z-x = \-sT (1.7) and _ 1 - s{T/2) z -i + *(r/2)- (L8) The latter is known as the bilinear transformation. The above transformations are based on the approximate prediction equations given respectively by r(nT) = r((n - 1)T) + Dr({n - 1)T)T (1.9) and / ^ // ^ Dr([n-1)T) +Dr(nT) r{nT) = r{{n - 1)T) + â€”^ K- -T. (1.10) where r is either the input signal u or the output signal y. As expected, the bilinear transformation usually yields better approximations than the first transformation which is based on approximating the derivative by a finite difference. Therefore, it can be concluded that all the methods of determining a continuous model from the Introduction 5 discrete equivalent model involve approximations either in the input signal or in the transformation from the continuous frequency domain to the discrete frequency domain. Another method, which assumes a zero input and consequently involves no approximation of the input signal between sampling instants, is known as Prony's method [Marple 87, Ch.ll]. Prony's method is a technique for modeling signals as a linear combination of exponentials 'W = E^' (1.11) from equally spaced samples. There are three basic steps in Prony's method. Step one determines the linear prediction model that fits the available data, i.e. determines the denominator polynomial of H<i{z~l). In step two, the obtained polynomial is factorized to obtain n roots and then obtain as in (1.6). Step three involves the solution of a set of linear equations to yield the estimates of the coefficients G\,. ..,Câ€ž. The method does not consider estimating the numer ator polynomial of H<i{zâ€”l) and therefore it is suitable only for identification of systems with impulsive type inputs [Noharski 85]. Application of Prony's method in identification of electromagnetic systems can be found in [Schaubert 79] and [Poggio 78], and in biological systems can be found in [Crittenden 83]. A state space version of Prony's method is discussed in [Crittenden 83]. The system dynamic matrix A is obtained through a logarithmic operation similar to (1.6) under the assumption that the full state vector is measurable. The method does not give the vector 6. The given reference discusses the choice of the sampling period in conjunction with the finite precision arithmetic, used to represent data samples, and the fastest and the slowest exponentials in the analyzed signal. All the above discrete models for characterizing continuous systems and signals from discrete measurements are based on the assumption that the sampling time T is fixed, and that the signal is band limited. The choice of T for a wideband signal is critical because a too small T could result in numerical instability in the estimation equations of the linear prediction model parameters and in the inverse z-Transform, while a too large T gives rise to the problem of aliasing [Astrom 83, Sec.2.5]. Introduction 6 1.2.2 Modulating Functions A classical technique, [Shinbort 57], similar to the Fourier transform and the Laplace transform, is used to transform a differential equation to a set of linear algebraic equations in the unknown parameters. Shinbort used the name method functions but the name modulating functions was used later [Loeb 65]. The differential equa tion (1.1) is multiplied by a vector f_M{t) of modulating functions and integrated over the period of observation t0 â€”â€¢ tf ff Pny(t)lM(t)dt + $>,. f Djy(t)lM(t)dt = X>; P Dju(t)lJt)dt. (1.12) Performing the integration by parts and choosing such that ^LM^o) = V'htkf) = 0 ,J = 0,..., n - 1 (1.13) yields an equation in the measured variables only and not in their derivatives (-1)" f y{t)VÂ»^dt + Zi-lYaj Py^^dt J tâ€ž j_Q Jtâ€ž = E(-l)^ P u(t)D>lM(t)dt. (1.14) Unless (1.13) is satisfied, the terminal values of u, y and their derivatives would appear in (1.14). It is seen from (1.14) that obtaining the estimation equations for the unknown parameters requires only the computation of the scalar products of the recorded signals with the vector of modulating functions and its derivatives. The use of orthogonal functions, e.g. Fourier functions, is preferred to achieve linear indepen dence of the estimation equations, a requirement which may violate the terminal condition (1.13). The elements of / can be derived from the impulse response of a set of analog filters such that the output of these filters provide the input data to an estimation algorithm which estimates the unknown parameters from (1.14). Such analog filters are known as state variable filters [Young 70]. The Poisson moment functional approach to continuous system identification [Saha 83] is a special case of the state variable method. A Poisson filter chain Introduction 7 is composed of cascaded identical filters with a first order transfer function j^, X >0. Beside not being orthogonal, the impulse responses of the cascaded filters obviously do not satisfy the terminal condition (1.14) so that the unknown initial conditions are introduced to the estimation problem. Also, passing the measured signals through many filter stages results in smoothing the signals and consequently reducing the information content of the filtered signals. Pearson [Pearson 85] used a projection of a vector of Fourier functions in order to satisfy the terminal condition, and applied the method to identification of bilinear systems using sinusoidal probing signals over sequential time intervals. 1.2.3 Multiple Integrations and Orthogonal Functions In recent years, a new technique has been established for the solution of continuous system identification problem using orthogonal functions. The main feature of this technique is that, like the modulating function technique, it reduces the problem to that of solving a system of linear algebraic equations, thus simplifying the solution. The transformation of the differential equation to an algebraic equation takes place in two steps. First, the differential equation is transformed into an integral equation by means of n integrations3. Second, the input-output signals are correlated with a vector f^{t) of orthogonal functions so that multiple integrals of the input-output signals can be obtained in terms of f0{t)- The key idea of this technique is based on the following property of a basis vector fo{t) under multiple integration f â€¢â€¢â€¢ r uÂ°)d<j â€¢ â€¢ â€¢ dtÂ«-i * pkLw u-15) JtQ J to k where P is a constant square matrix, called the operational matrix. The elements of f_o(t) 'are orthogonal on the observation interval t0 â€”> tf. Clearly, P depends on the particular choice of the basis vector fo{t). Since the method was used for optimal control synthesis using Walsh functions [Chen, C. 75], several such basis functions have been in use in solving the continuous system identification problem. Walsh functions [Rao 75] and [Bohn 82], block pulse 3 Similar to a Poisson filter chain with A=0. Introduction 8 functions [Sannuti 77] and [Chen, W. 87], Laguerre polynomials [Hwang 82], Legen-dre polynomials [Paraskev. 85], generalized orthogonal polynomials [Wang 87a] and Fourier functions [Chung 87] are examples of orthogonal functions used for linear time-invariant system identification. Because of their properties under multiplica tion, some orthogonal functions have been used in identification of time-varying systems [Tzafestas 77] and [Rangana. 87], and in identification of bilinear systems [Karanam 78] and [Wang 87b]. For example, a set of N = 2m Walsh functions of sequencies 0,..., N â€” 1 form a closed set under multiplication. The product of any two functions in the set results in a third function in the same set. This prop erty allows simple representation of the products input-output signals and also the products of input-output signals with time functions that correspond to the time varying parameters. Among these sets of orthogonal functions, Walsh functions have received much of the attention, which is apparent from the relatively large number of publica tions on the use of Walsh functions [Rao 83]. The reason for this seems to be the simplicity of the Walsh transform and of integration using the Walsh operational matrix. It has been shown, [Bohn 82], that integration using the Walsh operational matrix can be reduced to binary shifts and addition of data samples. However, nu merical errors accumulate because of multiple integrations while truncating the operational matrix. To reduce these errors, one-shot operational matrices for in tegration were introduced [Rao 81]. Another drawback of the method is that the unknown initial conditions are introduced into the estimation problem as a result of multiple integrations. These additional unknowns have to be estimated together with the system parameters which may significantly increase the dimension of the estimation problem. 1.3 Nonparametric Identification Although nonparametric system identification has lost popularity because of the current emphasis on control synthesis and design tools which hinge upon parametric system models, it has gained some ground in recent years because of the growing interest in long-range predictive control, [Richalet 78], [Rouhani 82], [De Keyser 85], Introduction 9 [Bars 85] and [Clarke 87], which mostly utilizes estimates of the system impulse or step response. 1.3.1 Frequency Response Obtaining the frequency response of the system through sinusoidal input testing is the most robust of all system identification methods. However, the long test time, usually multiples of the response settling time, required to obtain each relevant point on the frequency response curve is the reason why digital spectral analysis using the discrete Fourier transform now replaces the conventional sinusoidal input testing. To obtain the frequency response of the system, the input signal u(t) and the output signal y{t) are uniformly sampled to obtain u3(k) and y3(k). Then, the sampled signals are Fourier transformed to obtain the spectral components U(ku) and Y{kw) where OJ is the fundamental frequency over the data window. Samples of the frequency response are obtained as the ratio of the corresponding input and output spectral components or alternatively as the ratio of the input-output cross-correlation and the input auto-correlation where the astrisk denotes the complex conjugate. Unless y is the periodic response to u, the frequency response estimates obtained from (1.16) and (1.17) are subject to systematic errors, namely windowing errors and variance errors [Wellstead 81]. One way to reduce these errors is to average the input-output Fourier coefficients over a large number of independent data windows. However, sampling of the input-output signals could give rise to aliasing problems especially if the signals are wideband. Recently, [Douce 86], a method was suggested to apply a periodic waveform to the system for a time of a few periods and then remove it. It is known that if the system starts from rest, then the transient component of the response due to the input turn-on is equal and opposite to the transient component due to the Introduction 10 input turn-off. The Fourier analysis of the two transient components yields results which sum to zero, and therefore gives correct estimates of the periodic response. However, this method requires response measurements over a relatively long time interval, one equal to the input waveform period plus the response settling time. Advantageous use of this method seems possible only in the case of highly resonant systems where applying a periodic waveform to the system for a long time results in wild system response. Different methods have been devised since the late fifties [Levy 59] to obtain a parametric model of the system from its frequency response estimates by fitting those estimates to a complex transfer function. A comparison of five related linear least squares methods for transfer function synthesis from frequency response data is given in a recent publication [Whitfield 86]. 1.3.2 Impulse Response The possible damage inflicted by direct use of an impulsive test signal on the system hardware and the presence of output noise have led to the decline in its use. The currently accepted approach is to obtain the impulse response of the system by a process of deconvolution. Deconvolution is the inverse of the convolution of the system output, or the input-output cross-correlation, with the system input, or the input auto-correlation, respectively where h is the system impulse response and where Ruy and Ruu are, respectively, the input-output cross-correlation and the input auto-correlation. Theoretically, white noise is an optimum test input for estimating the system impulse response because its auto-correlation is an impulse. Thus, Ruy is itself an estimate of the impulse response. Unfortunately, white noise is difficult to realize in real life. The simplicity of generating a Pseudo Random Binary Sequence (PRBS), beside having an auto-correlation that approximates an impulse over one period of the sequence, gave a boost to its use in nonparametric system identification. y{t) = h{t) * u{t) (1.18) Ruy{r) = h{r) * ^{T) (1.19) Introduction 11 However, the auto-correlation of a PRBS is periodic, with finite width triangular pulses and nonzero value between pulses which make it just an approximation of an impulse. These three differences from the ideal impulse lead to aliasing problems, smearing and bias in the estimated impulse response [Wellstead 81]. 1.4 Scope of the Thesis The thesis considers the problem of identifying a linear continuous-time time-invariant model of a single-input single-output dynamical system from the available measurement taken from the system input and output signals. A deterministic model is used to describe the dynamics of the system. No doubt a stochastic system model that includes a noise model to account for the effect of unknown dis turbances on the system could be useful, especially if a priori knowledge of the noise model is available. However, considering a stochastic model makes the identifica tion problem more difficult especially in the continuous-time case. To account for the effect of unknown system disturbances, the method of instrumental variables, which utilizes minimum a priori knowledge of the noise model, is used. A new approach to the identification problem is considered based on the def inition of continuous-time functions, given the name modal functions, defined in terms of the system output, the system output derivatives and the system state. The modal functions are defined in chapter 2 and their properties are studied. Associated with the definitions of the modal functions are the modal parameters. Estimating the modal parameters from the available input-output measurements is considered in chapter 3 with numerical examples. The method of instrumental variables is used in simple cases. The relations between the modal functions and the modal parameters with the desired system model parameters are given in chap ter 4 in theorem form and then used to identify the model parameters. Use of the modal functions and the modal parameters in nonparametric system identification, namely frequency response, step response and impulse response, is considered in chapter 5. Chapter 6 describes a design procedure to obtain a minimum-time ob server for the continuous system state based on the use of the state modal functions and parameters. Finally chapter 7 concludes by summarizing the research results Introduction 12 as well as suggesting some points that are not covered by the thesis and that could be suitable for future research. The thesis deals mainly with developing and proving theoretical algebraic re lations between modal functions, modal parameters, input-output data, and the desired unknowns of the system models. These relations are in the standard linear-in-the-parameters estimation equations form. Consequently, even though some il lustrative numerical examples are given, the thesis does not deal extensively with the use of well established numerical techniques for solving linear estimation equa tions and the known methods for reducing parameter bias due to noise. The thesis places emphasis on showing how the modal functions approach offers new options in parameter estimations and in state observer design. Chapter 2 Definitions and Properties of Modal Functions 2.1 Introduction The goal is to identify parametric and nonparametric models of the system; the tool is the modal functions and the modal parameters which will be defined in this chapter. Section 2.2 provides the mathematical background needed in subsequent sections. The modal functions and the modal parameters are defined in section 2.3 and the properties of the modal functions are studied in section 2.4. 2.2 Mathematical Background It is well known that the general solution of the state equation (1.2) is given by [O'Reilly 83, page 3] (2.1) where <f> is the transition matrix of (1.2) given by 4>{t) = e At (2.2) The following properties of <j> will be used in subsequent analysis: <j>{h + t2) 4>{ti)<t>{h) D<f>(t) = J M{t) (2.3) 13 Definitions and Properties of Modal Functions 14 The vector x(0) in (2.1) is the system state vector at time t = 0. Similar to the initial conditions y(0), Dy(0),..., D n~ 1y(Q) required to obtain the solution of the differential equation (1.1), the elements of x(0) are necessary to obtain the solution of the equivalent state equation (1.2). The terms state and state variables are used to refer to x and the elements of x respectively. The fact that the input-output relation is invariant under a linear transforma tion of the state variables given by x(t) A b = T- Xx{t) T-iAT T~ lb c TT (2.4) implies that the parameters of A, b and c T cannot all be determined from the input-output measurements . To be able to obtain unique solutions to state model identification problems, it is necessary to find models which contain the smallest number of parameters, i.e. canonical forms. Two canonical forms that will be used in the thesis are the Jordan canonical form and the observable canonical form. (l) Jordan canonical form: The matrix A can be transformed to Jordan canonical form A* by the transforma tion A* = M~ lAM (2.5) where M is the modal matrix [Sinha.P.K. 84]. The matrix A* is block diagonal, or diagonal for short, with diagonal blocks of the form sk 1 O ,k = 1, (2.6) â€¢â€¢ 1 O sk where n* is the number of distinct eigenvalues of A and the dimension of Ak is equal to the multiplicity of sk. The vector x* = M~ lx will be refered to by the vector of system modes. Usually, the term modes is used to refer to the eigenvectors of a matrix. For the purpose of convenience only, the term modes will be used to refer to the state variables of x*. Definitions and Properties of Modal Functions 15 (2) Observable canonical form: If the pair [A,c) satisfies the observability condition [O'Reilly 83, page 9], a non-singular matrix M can be found that transforms the matrix A into phase variable form1 given as A = M~x AM = 0 0 O 1 (2.7) such that c TM 1 0 Â°] (2.8) The elements of the bottom row of A are the coefficients of the characteristic equation of A given by (2.9) A n+J2 aiAJ = 0 3=0 where O is the zero matrix. The corresponding arrays A, x, i and ~<f> are defined by relations similar to (2.4) using the substitution T = M. The state model in the phase variable form has a very close relation to the input-output model (1.1). The elements of the bottom row of A are the parameters a0,... ,aâ€ž_i of (1.1). The phase variables, i.e. the elements of x, are defined in terms of the output y and its derivatives. It can be shown by differentiating (2.1) j times that D'y(t) = c TV jx(t) = c TA'x{t) + j^^A^bD^'uit) J > 0. (2.10) Â«=i Equation (2.10) is invariant under a transformation of the state variables. Then it follows that substituting for A, b, c and x by A, b, c and x, respectively, and using (2.7) and (2.8) gives y{t) = hit) t=i V ny{t) = a Tx{t) + T,hDn'iu{t). , 0 < j < n (2.11) t=i 1Also known as companion form. Definitions and Properties of Modal Functions 16 Finally, the relation between the elements of vector b and b0,..., 6n_i of the input-output model follows by substituting into (1.1) from (2.10) and using (2.9) as n bi = Y, akh-j 3 = 0, â€¢ ,n - l,an = 1. (2.12) k=j+l 2.3 Definitions of Modal Functions and Modal Parameters Let fa be a scalar function defined as a weighted sum of the state x at the distinct time instants t, t â€” Ti,..., t â€” Tn fa{t) = a Tx(t) +J2PÂ°i c Tx{t - Ti) (2.13) t=i where a is a vector of dimension n, and pal,... ,pan are constants chosen such that fa is independent of the instantaneous value of the state x. From (2.1) it can be shown that x{t - Ti) = 4>{-Ti)x{t) + j* T' <f>{t -Ti- T)bu(r)dr. (2.14) Substituting for x at different time instants from (2.14) into (2.13) gives, after arranging terms /Â«(<) = [Â« T + !>Â«.â€¢ cTH-Ti)]x{t) +J2pai f~ T' cT<f>{t -T- T)bu{r)dT (2.15) i=i .=1 Jt The condition that makes fa independent of x(t) for all time t is or where o T is a vector of zeros. For this condition fa is given by rt-Ti /Â«(*) = EPÂ« / -Ti- T)bu{r)d7 t=i Jt Let Pa (2.16) (2.17) (2.18) such that (2.16) becomes a T + p T cT^(-ro CTH~Tn) T = O . (2.19) Definitions and Properties of Modal Functions 17 A necessary condition for the existence of a vector pa of dimension n that satisfies (2.19) is given by the following theorem [El-Shafey 87]. Theorem 2.1 : A vector pa, of dimension n, that satisfies (2.19) exists only for observable systems. Proof: The existence of pa depends on the invertibility of matrix xjj given by cTH-Ti) 4,= cTH-Tn) _ (2.20) In order to study the invertability of ip, the diagonal form of (j> is used. Substituting for <f> by M<f>*M~x and for cT by c*TM_1 gives c^Pi-T^M-1 (2.21) c*T<j>*{-Tn)M-x From matrix multiplication rules, if) can be put as the product of three matrices 0 = V D(c*)M~1 (2.22) where D(c*) is a block diagonal matrix with n* diagonal blocks; each diagonal block corresponds to a distinct eigenvalue of A and is given by (2.23) O ckl where nk is the multiplicity of the corresponding eigenvalue. The matrix V is composed of arrays Vk of dimension nkxn, k=l,... ,n* given by (eSkTl â€¢ â€¢ -eSkT") V, J-{e3kTl â€¢ â€¢ â€¢ ea*T") d"k-\ i 3tr, . .,pskTn\ L K^7^ 6 I (2.24) stacked vertically on top of each other. It is clear from (2.22) that a necessary condition for the invertibility of 0 is that D(c*) is invertible which means that Ck, k=l,... ,n* are invertible. This condition is Definitions and Properties of Modal Functions 18 exactly the necessary and sufficient condition for the observability of the system in Jordan canonical form [Sinha.P.K. 84, page 173]. The matrix V is of Vandermonde type with known determinant expansion with uniform 2i,... ,Tn. This completes the proof. It follows from theorem 2.1 that if the system is observable and 0 is invertible then there exists a vector pa of dimension n given by Â£ = -QLT^X = -aTMD~1(c*)V~1 (2.25) whose elements satisfy (2.16). Different choices of a lead to different fa and pa. Particular choices of a are of interest as for example aT = cT = c*TM"1 which leads to the following definition. Definition 2.1 : For an observable system of order n, the output modal function y0 is defined as: n Vo{t) = y[t) + X>w-y(t - (2.26) Â»=i where pol,... ,pon are the output modal parameters given according to (2.25) by Â£ = -cT0-x = -c*TD-l{c*)V~\ (2.27) The name modal time shifts will be given to Ty,..., Tn. In later chapters it will be shown how the output modal function can be used in parametric and nonparametric system identification. Another choice of aT is cTAJ, 0 < j < n. It is known that cTylJx is invariant under a linear transformation of the state x. For example, A*x â€” c*TA*3 x* = cTAx. It can easily be verified that cTAj = ij+1 ,0<j<n (2.28) where is the [j + l)th row in the unity matrix. As a consequence of (2.28), the (j + l)th state variable in x. This leads to the following definition Definition 2.2 : For an observable system of order n the state modal function xjj, 0 < j < n is Definitions and Properties of Modal Functions 19 defined as y;-(t) = ~xj+1{t) + X) PiiV{t ~Ti) , 0 < j < n (2.29) where pji,... ,pjn are the state modal parameters given according to (2.25) by f. = -cT^'V_1 = -c?TA*jD-\c*)V-1 ,0 < j < n (2.30) It is seen that the state modal function y0 is the same as the output modal func tion given by definition 2.1. The state modal parameters will play a key role in parametric system identification and in the design of state observers as will be seen in later chapters. Other modal functions of use in system analysis are given in the following definition Definition 2.3 : For an observable system of order n, the output derivative modal function is defined as n Vj{t) = Pjy{t) + YIpM* ~ TÂ«) Â»0 < 3 (2-31) i=l where pji,... ,pjn are the same state modal parameters given by (2.30). The name output derivative modal parameters will be used as well to refer to the state modal parameters pji,... ,pjn. There is a direct relation between a state modal function y, and the correspond ing output derivative modal function y; for 0 < j < n. Subtracting (2.29) from (2.31) gives yj{t)-yj{t) = t>jy{t)-xj+i{t) ,0<j<n. (2.32) From (2.11), relating output derivatives to state variable in x, and (2.32) it follows that j Vii*) -Vi(t) = Y^biD'^uit) ,0<j<n. (2.33) Equation (2.33) will prove useful in estimating the vector 6. 2.4 Properties of Modal Functions After defining different modal functions it is suitable to investigate some of their properties. Equation (2.17) is a general formula for defining different modal func-Definitions and Properties of Modal Functions 20 tions given by definitions 2.1 or 2.2. The output modal function or the state modal function is obtained by substituting for pai,... ,pan by the corresponding modal parameters given by (2.30). Using the substitution of variables o = t â€” T â€” T in the integration in (2.17) and interchanging the limits of integration yields Comparing (2.1) with (2.34), it is seen that the system general response is composed of two parts; one which depends on the system initial state, and the other which is the convolution integral of the input and the system impulse response. Conversely, a modal function is defined to be independent of the state, also it is clear that the integrals in (2.34) are different from the standard convolution integral. These are the reasons that make a modal function possesses properties which are different from the general system response. Four properties of modal functions in connection with input superposition, scal ing, time shifting and time differentiation are given below. Let the pair (fa,u) satisfy (2.34). The following properties are obvious from 1. Superposition: If the pair (Afa, Au) satisfies (2.34), then (fa + Afa, u + Aw) also satisfies (2.34). 2. Scaling: The pair (afa,au) satisfies (2.34) where a is a scaling factor. 3. Time shifting: The pair (fa(t + At),u(t + At)) satisfies (2.34) where At is a time shift. 4. Time differentiation: The pair [Vfa,Du) satisfies (2.34). It is important to note that the state vector, and consequently any linear function of the state variables, e.g. the output, do not possess any of the above properties. The following analysis shows the relation between the modal functions and both the system modes and the input modes. Let u be the free response of a continuous controller of order nc given by Â«=i (2.34) (2.34). Dxc{t) = Acxc{t) (2.35) Definitions and Properties of Modal Functions 21 fa{t) = Â£>a, <? f T' <t>iÂ°)!> Â£M-Ti - o)doxc{t). (2.37) Â£ = X>- <? HÂ°)!> c^d-Ti - o)do (2.38) 1 JO Ue(<) = <LTcxc{t). (2.36) Substituting for u in (2.34) from (2.36), noticing that ^.{t â€” Tt) = <j>c{â€”Tj)xc(i) where <j)c is the transition matrix of the controller, yields Pat C1 I 1 ' Let vector q be â€”r,-Pat CJ / .=1 JL Note that is of dimension nc as the number of modes in the input. Substituting from (2.38) into (2.37) yields fait) = qTaxcit) (2.39) The important conclusion that is drawn from the last equation is that a modal function is just a linear function of the modes of the input, and it contains none of the modes of the system. Equation (2.39) gives a general modal function for arbitrary a in terms of the input modes. It follows from (2.39) and the definitions of the output modal function and the state modal functions (2.29) that yj{t) = xj+1it)+J2pjiy{t-Ti) = q.xc{t) ,0 < j < n (2.40) t=i where q is given according to (2.38) by ?< = zZ'PH / ' cT<t>iÂ°)!> gtd-Ti -o)do ,0<j<n (2.41) i=l Similarly, an equation like (2.40) can be obtained for the output derivative modal functions. From (2.35) and (2.36), by differentiating j times, it follows that Viu{t)=gAixc{t). (2.42) Substituting from (2.42) into (2.11) yields i Djy{t) = x;-+i(0 +52bigA'-%it) ,0<j<n. (2.43) t=i Definitions and Properties of Modal Functions 22 Then from (2.40) and (2.43) it follows that n Vj{t) = Djy{t) + J2pM* - T<) = ,0 < i < n (2.44) Â«=i where q is given by T \ g + ELx<>iÂ£xr ,0<j<n q-=\i ,y = o (2-45) Equations (2.40) and (2.44) have a continuous time auto-regressive moving average (ARMA) structure. Unlike the well known discrete-time ARMA modal, which describes the system behavior only at the sampling instants, (2.40) and (2.44) are true for continuous time t. Uniform intersample spacing is a characteristic of the discrete-time ARMA model, while as seen from (2.40) and (2.44) Ti,T2 â€” 2\,..., Tn â€” Tn_i are not necessarily equal. The intermix between the state variables and the output in (2.40) gives it the form of a state observer which gives the system state as a function of input-output measurements. The following numerical example gives the modal parameter vectors from a known system dynamic matrix. Example 2.1 : Consider a third order system for which the dynamic matrix in diagonal form and the vector x*T are given by -1+J5 0 0 0 -1-J5 0 0 0-2 c*T = I 1 1 1 and Obviously, the system is observable. The poles of the system are the diagonal elements of A*. It is desired to obtain the modal parameters of the output and of its first, second and third derivatives for the modal time shifts 2i = 4ors T2 = 80TS T3 = 120T3 where Ts = 7r/420. The results of computing p., j=0,... ,3, according to (2.30) are listed in table 2.1. Definitions and Properties of Modal Functions 23 j 0 1 2 3 -0.6605E+00 +0.2999E+01 +0.5715E+01 -0.7848E+02 +0.6106E+00 -0.2564E+00 -0.1476E+02 +0.3497E+02 -0.3022E+00 -0.4330E+00 +0.5721E+01 +0.5819E+01 Table 2.1: Modal parameter vectors 2.5 Summary The main results for this chapter are the continuous-time ARMA models (2.40) and (2.44) relating the state and the output derivatives to input-output data. From a practical stand point it is seen that a choice must be made for the system order n and the modal time shifts. The choice of model order is a problem common to all parametric methods whether discrete or continuous. One approach, commonly used when estimating the model parameters, is to increase the model order until there is no significant decrease in the discrepancy between the real system and the estimated model data. A second approach is to over-parametrize and then eliminate relatively small pa rameters. The obtained continuous-time ARMA models are valid for nonuniform time shifts, and this allows greater flexibility in choosing 7^. Such flexibility is not possible in the methods that are based on the standard discrete ARMA models such as Prony's method. The flexibility in choosing the time shifts allows for determining both the fast and the slow modes in the system response. Chapter 3 Estimation of The Modal Parameters 3.1 Introduction In chapter 2 the modal functions of the system output, the output derivatives and the state variables were defined and their properties were studied. Associated with each modal function is a modal parameter vector whose dimension is equal to the order of the assumed observable model of the system. So far the modal parameters were obtained through equations (2.27) and (2.30) which assume the knowledge of the system modes. This chapter deals with the problem of estimating the modal parameters assuming a linear time-invariant observable system of order n and the availability of input-output measurements only. The estimation problem is formulated in section 3.2. Section 3.3 presents with numerical examples the estimation of the output modal parameters. Estimating the output derivative modal parameters introduces the problem of unavailable output derivatives so it is presented separately in section 3.4. The problem of increased dimension of the unknown parameter vector when multi-mode inputs are used is discussed in section 3.5 with a suggestion of a solution using a measurement prefilter to partially annihilate the input modes. 3.2 Problem Formulation From the definition of modal functions in chapter 2 we arrived at equation (2.44) which has a continuous-time ARMA structure. It is desired to obtain estimates of the parameter vectors p. and q. for j > 0 from the available input-output 24 Estimation of The Modal Parameters 25 measurements. The approach to the problem is to minimize some scalar function of the equation error which is given according to (2.44) by e,-(t) = Djy{t) + Y,Pjiy{t - Tt) -J2gjixci{t). (3.1) Â«=i where qji and the elements of q . and x^. respectively. Let Oj be the vector of unknown parameters and r(t) be the vector of measurements given respectively by and rT{t) = [y{t - TO â€¢ â€¢ â€¢ y{t - Tâ€ž) - xcl{t) xcn,:{t)}. Hence, the equation error in (3.1) is given in terms of Oj and r by ej{t)=Viy{t) + rT[t)0_]. (3.2) (3.3) (3.4) An estimate 0_j of the unknown vector Oj is sought such that a scalar loss function J, defined as some norm of the equation error e; over the observation interval t0 â€”> tf is minimized. The loss function Jj reflects the discrepancy between the continuous ARM A model (2.44) and the real system so that the estimate Oj that minimizes Jj can be used to obtain estimates of the output derivatives Vjy as Vjy{t) = -rT(t) Oj. (3.5) The most common cost function of the continuous-time error e; is based on the integral of a weighted norm in tj. However, a cost function that is formulated as the sum of squares of samples of tj at discrete time instants is possible. In the following two sections estimates of the unknown parameters are obtained by minimizing cost functions which are formulated in terms of discrete-time and continuous-time measurements, respectively. 3.2.1 Estimation from Discrete Measurements In the case of discrete measurements , let N > n + nc be the number of samples obtained over the observation interval at the discrete time instants tk, k â€” 1,..., N. Estimation of The Modal Parameters 26 The least squares estimates of the unknown parameters are obtained by minimizing the sum of squares of samples of the equation error at tk, k = 1,..., TV, that is, it is required to minimize Jj<i = Y,e2j(h) = eJdejd. (3.6) k=i In terms of the measurements , Jjd is given by Jjd = (zjd + Rd ej)T{z^d + Rd 6j) (3.7) where Zjd= [ vjy{h) â€¢â€¢â€¢ Djy{tN) }T (3.8) and Rd has rows i%Rd, k = 1,..., TV, given by llRd = rT{tk). (3.9) It is known that the estimate 6_jd of the unknown parameter vector 9j that minimizes Jjd is given by ejd=-{RTdRd)-1RTd zjd. (3.10) Obtaining estimates of the unknown parameter vector from (3.10) seems possible only for the case of j = 0. This is because although Rd is available by direct measurements , samples of the output derivatives in Zjd, j > 0, are not available. 3.2.2 Estimation from Continuous Measurements Because the system is continuous and can be sampled, at least theoretically, at any frequency, a large amount of data could result from the sampling process. Therefore, some form of data compression is needed. Correlating the measured signal with a function fk over the observation interval results in a single number. Correlating the measured signal with a vector / of N functions results in the con centration of the information contained in the collected data into N components, where each component represents a best fit to the measured signal in a least squares sense. Fourier functions and Walsh functions are two suitable candidates for this purpose because each of them forms a set of orthogonal functions. From the available continuous measurements the loss function is formulated as N J>f = Y,ehk=4f^f (3-H) k=l Estimation of The Modal Parameters 27 where t^f is *ne correlation between the equation error tj and the vector of functions / over the observation interval 1,7 = P l{t)ej{t)dt. (3.12) J to In terms of the available measurements Jjf is given by Jjf = (zjf + RdjfUjf + Rflj) (3.13) where [tf l{t)P jy{t)dt (3.H) and Rf= [*'f{t)rT{t)dt. (3.15) The estimate 0Jf- of the unknown parameter vector 03 that minimizes Jjf is given by 4,7 = -(RfRf)-lRfZs/- (3-16) Equation (3.16) seems suitable for obtaining estimates of the unknown parameter vectors for j > 0 as well as for j' = 0 unlike (3.10) in the discrete case which is suitable only for j = 0. This is because of the integration in (3.14) that gives z^. Later in this chapter it will be shown that an estimate of z_jf can be obtained either from an available estimate of zy_^f, or from an estimate of D 3~ ly obtained from (3.5) using an available estimate of 3.3 Estimating The Output Modal Parameters The output modal parameter vector po is estimated as a part of 0O. An estimate of d0 can be obtained either from discrete data, 0_od in (3.10), or from continuous data, O^f in (3.16). By putting j â€” 0 in (3.8) and (3.14) it becomes clear that zod and z0f required in (3.10) and in (3.16) respectively are obtained from the measured output. The matrices Rd and Rf are obtained from the input-output measurements. Estimation of The Modal Parameters 28 3.3.1 Estimates from Discrete Measurements This is the simplest case because Rd and Zgd are formed from input-output samples. Input-output samples should be collected at a rate which is at least equal to twice the highest frequency in the output or in the input whichever is greater. Example 3.1 : Consider a third order system with a transfer function U[S) ~ s* + 4s2 + 305 + 52 The system has three poles at Si = -2 s2 = -1 + j5 53 = -1 - j5 and a zero at 54 = â€”4. The system has the same poles as the system considered in example 2.1. It is desired to obtain the output modal parameters for a third order model of the system using discrete input-output measurements . To start with a simple case let the input u contain a single mode, a constant for example. In this case the continuous ARMA model is given according to (2.44) by 3 J/(0 + Y^Poi y[t - Ti) = q00u0{t) i=i where uQ is the constant input and where q00 is the single parameter which cor responds to the input single mode and which will be estimated together with the output modal parameters. Assuming a zero initial state of the system the response of the system to a step input of magnitude +1 followed by another step input of magnitude â€”1 is computed each sampling time Ts = 7r/420. The system poles are not far apart so it is reasonable to assume uniform modal time shifts ; Ti = T2 â€” Ti = T3 â€” T2 = 40Tâ€ž. The continuous ARMA model in this example is valid only for data windows wi : Tz < t < tt, where tt = 444TS is the transition instant when the input switches from +1 to â€”1, and w2 : tt + T3 < t < tf. Therefore, the sampling instants tk, k = 1,..., N should belong to either wi or w2. Twenty estimation equations are formed Estimation of The Modal Parameters 29 for 20 different tk\ 10 in wi and 10 in w2. In each window the spacing tk â€” tk_i is uniform and equals 20T3. Estimates of the unknown parameters given by (3.10) are listed in table 3.1. The estimated values of the three output modal parameters can be compared with the true values computed for the same modal time shifts listed in table 2.1. -0.6605E+00 Id +0.6106E+00 -0.3022E+00 +0.6480E+00 Table 3.1: Estimate of 90 using discrete data. 3.3.2 Estimates from Continuous Measurements Estimates of the unknown parameters are obtained from continuous data through (3.16) with j = 0. As seen from (3.14) and (3.15), the estimation equations are formed by correlating the available input-output measurements with a vector of functions / over the observation interval. As mentioned earlier, Fourier functions and the Walsh functions are suitable for this purpose because each forms a set of orthogonal functions. Let T0 be the length of the observation interval t0 â€”> tf and let t0 = 0 such that tj = tf â€” t0 = T0. Fourier functions defined over a time period 0 â€”> T0 occur in pairs, sinw^t and cosw^t where uk = ku>0 and ui0 = 27r/T is the fundamental frequency. A vector of N = 2m Fourier functions defined over the time period 0 â€”+ T0 is given by fF(t) = [ s'mu>0(t) â€¢â€¢â€¢ s'mmu!0t cosw0(t) â€¢â€¢â€¢ cosmw0ij. (3-17) On the other hand, the vector of Walsh functions defined over the time period 0 â€”> T0 is given by Â£w(t) = [ wal(0,t/ro) â€¢ â€¢ â€¢ wal(iV - l,t/T0) ] (3.18) Estimation of The Modal Parameters 30 where wa\(k, o) is a Walsh function of sequency k defined over the interval 0 < o < 1 [Tadokoro 78]. The correlation of a signal with a Walsh function requires no multiplication as opposed to the correlation of a signal with a Fourier function. Only addi tion and subtraction of samples are needed. A Fast Walsh Transform (FWT) [Shanks 69]reduces the computational burden even further. In this case a number N = 2m functions, where m is integer, should be used. The time interval of length T0 is divided into N equal subintervals and the samples in each subinterval are averaged to obtain a vector of N elements. This vector of averaged samples is transformed into a vector of N Walsh components using the FWT. Example 3.2 : Consider the same system as in example 3.1 with the same input-output measure ments. The first step of data compression is to concatenate the data windows wx and W2 where the continuous ARMA model in example 3.1 is valid. The second step is to correlate the concatenated input-output measurements with a vector of N orthogonal functions over the period of orthogonality in order to form matrix Rf and vector z0f. Estimates obtained by employing a vector of six Fourier functions, and estimates obtained by employing a vector of eight Walsh functions are listed in table 3.2. (using Fourier functions) A 6-oW (using Walsh functions) -0.6605E+00 +0.6106E+00 -0.3022E+00 +0.6480E+00 -0.6605E+00 +0.6106E+00 -0.3022E+00 +0.6480E+00 Table 3.2: Estimates of 0O using continuous data. Both examples 3.1 and 3.2 use a single input with a single mode in order to keep the dimension of the estimation problem low. A constant input does not provide persistent excitation of the system, a characteristic desired to avoid ill-conditioned parameter estimation equations. To provide persistent excitation of the system, Estimation of The Modal Parameters 31 two different levels of constant input were used in the previous examples. Input transition from one level to the other provides excitation to the system modes such that the response of the system after the transition carries the required information about these modes. Parameter estimation using discrete or continuous input-output measurements, described in examples 3.1 and 3.2 respectively, can be used with a Pseudo Random Binary Sequence (PRBS); a commonly used test input which provides persistent system excitation, such that the resulting parameter estimation equations are likely to be linearly independent. However, the input does not necessarily have to be of binary type; a piece-wise constant input that remains constant at least for time Tn before a new transition can be used in the same way as in examples 3.1 or 3.2. This restriction on the frequency of the input transitions results from the fact that the continuous ARMA model is valid for t > Tn with a constant input. Later in this chapter examples with multi-mode inputs are considered and a solution to the problem of increased number of unknowns is suggested. Both examples 3.1 and 3.2 assume a deterministic environment where measure ments are perfect and there is no unknown disturbances affecting the system. In practice, measurements are often contaminated with noise and the system is subject to uncertainty caused by external unknown disturbances. Therefore it is important to investigate the effect of noise on the estimates of the unknown parameters. It is known that the estimates of the parameters of a regression like (2.44) obtained by least squares using noisy measurements will be biased to a degree dependent on the ratio of the covariance of the noise and the covariance of the noise free signal. The following example shows the problem of biased estimates obtained from discrete and from continuous measurements as described in examples 3.1 and 3.2. Example 3.3 : A normally distributed noise signal with zero mean and standard deviation 0.1 is added to the output samples used in both examples 3.1 and 3.2. The noisy samples are used in the estimation equations as has been done in the above examples. The results of least squares estimation using the noisy samples are listed in tables 3.3 and 3.4. Estimation of The Modal Parameters 32 -0.5832E+00 iod +0.5881E+00 -0.2825E+00 +0.7105E+00 Table 3.3: Estimate of 0o using discrete noisy data. ioF (using Fourier functions) ioW (using Walsh functions) -0.6370E+00 +0.6151E+00 -0.3052E+00 +0.6730E+00 -0.6317E+00 +0.6326E+00 -0.3017E+00 +0.7003E+00 Table 3.4: Estimates of 0O using continuous noisy data. By comparing table 3.4 with table 3.2 it can be seen that acceptable estimates resulted when continuous measurements were used. This is because correlating the noisy signal with a Fourier or Walsh function works to average out the noise which has zero mean. Good estimates with low bias are likely to occur if the Fourier or the Walsh functions used have high correlation with the noise free output while they have no correlation with noise. Conversely, the problem of biased estimate is noticeable in the case of discrete measurements as is clear from comparing table 3.3 with table 3.1. Fortunately, the noise induced bias can be removed by incorporating an instrumental variable which is highly correlated with the true value of the signal, but not with errors in observation caused by various types of additive noise. One possibility of an instrumental variable is a shifted output [Soderstrom 81]. The amount of shift is chosen small enough so that the shifted output still has high correlation with the unshifted output, while this amount of shift is large enough so that the noise in the shifted output has no correlation with the noise in the unshifted output. The instrumental variable estimate of the unknown parameter vector is given by L = -{VTRd)-lVTzod (3.19) Estimation of The Modal Parameters 33 where V has rows ikV, k = 1,..., N, given by llV = \y{h - 2i + r) â€¢ â€¢ -y{tk -Tn + r) - xcl{tk) xme{tk)] (3.20) and T is the instrumental variable shift. The only difference between iÂ£V and t[Rd, in (3.9) that gives the least squares estimate, is in the shift r of the samples of y. Example 3.4 : It is desired to obtain an instrumental variable estimate of 0o using the same noisy input-output data as in example 3.3. Equation (3.19) is used to obtain the desired estimate with the instrumental variable shift r = 257V. Theoretically, any small amount of shift is enough in the case of white noise. The results are listed in table 3.5. It is clear that the instrumental variable method has worked successfully in reducing the estimate bias. 3.4 Estimating The Output Derivative Modal Parameters In this section the problem of estimating the unknown parameter vectors for j > 0 is considered. Obtaining least squares estimates of 6j for j > 0 from (3.16) would have been straightforward if Zjf for j > 0 were available by direct measurements . In most practical situations the number of measured variables is less than the number of state variables and therefore it is necessary to have some means of obtaining all state variables from the available measurements . Output derivatives, regarded as state variables, can be obtained by employing so called modulating functions. The use of modulating functions, as will be shown, allows the expressing of high order derivatives of the output in terms of lower order derivatives and ultimately in -0.6650E+00 +0.6346E+00 -0.3075E+00 +0.6537E+00 Table 3.5: Estimate of 0O using an instrumental variable. Estimation of The Modal Parameters 34 terms of the measured output. The following two sections show the use of Fourier functions and Walsh functions as modulating functions. 3.4.1 Estimates Using Fourier Functions Consider equation (3.14). Let the vector of modulating functions / be the vector of N = 2m Fourier functions / given by (3.17), and let the integration in (3.14) be from 0 to T0. Performing the integration by parts yields rTâ€ž The vector f_p satisfies the following relations liF = -/^P/jWP'-yO*^ (3.21) /F(0)=/F(To) = [0...0 l-.-l]3 P/F(*) = n/F.(t) where o n -ft o and J7 is a diagonal matrix given by n = o O mu0 Substituting from (3.22) and (3.23) into (3.21) gives where (3.22) (3.23) (3.24) (3.25) (3.26) cy_! = D^yiT,) - P'^yiO). (3.27) Equation (3.26) is recursive in the index j so that it can be used to obtain the spectral components of the j th derivative of the output in terms of the spectral components of the (j â€” l) th derivative and c;_i. Except for c0, cy_i is not available from the measurements . Even in the case of c0 = y(T0) â€” y(0) it is not wise to let the result of (3.26) depend heavily on the difference of two measured samples especially in a stochastic environment. However, cv,_i can be estimated with other Estimation of The Modal Parameters 35 unknowns in Oj. Hence, only an estimate c;_x can be used in (3.26) which means that ZjF obtained this way is just an estimate of the spectral components of the jth derivative of the output and not the true values. Substituting with the estimate of z_jF obtained from (3.26) into (3.13) gives JjF = ( â€”fi Z (j-l)F + RF â€”i ci-i )T(-0 Z.(J-\)F + RF -3 ci-i ) where RF = \RF , Â£F(o)}. (3.28) (3.29) Note that a hat is put on top of the variable to indicate that it is an estimate of the variable and not the true value. It is obvious that the estimate of (#J , c_,-i)r that minimizes JjF is given by o IjF [RFRF) XRF 0 (3.30) In order to obtain estimates for Oj and Cj_y for j > 0, equations (3.30) and (3.26) are used sequentially. Firstly, 0_x and c0 are obtained from (3.30) using zoF which is available from the measurements, then an estimate of z1F is obtained from (3.26) using Zgp and the estimate of c0. This estimate of z1F is substituted in (3.30) to obtain estimates for 02 and ci, and so forth. Example 3.5 : Consider the same system as in example 3.1. Let the system this time be excited by a sine wave input u(t) = sin4Â£. A sine wave has two imaginary modes thus the dimension of q., j = 0,..., n, is two and the continuous ARMA model is given by ^M*) + Y^Pi'Vi* ~ T>) = J2aiixci{t) = 9jisin4f + qj2cos4t. i=l i=l This model is valid for all time t > Ts unlike the case of a piece-wise constant input. The vector q. is dependent on the choice of the controller state variables xci â€”3 and xc2. A phase variable form is assumed such that r . iT xc = sui4t cos4i 0 4 -4 0 T [1 0] The system output is sampled each Ts â€” n/420. It is desired to obtain estimates Estimation of The Modal Parameters 36 j 0 1 2 3 -0.6617E+00 +0.2974E+01 +0.5757E+01 -0.7815E+02 IjF +0.6111E+00 -0.2529E+00 -0.1481E+02 +0.3491E+02 A ijF -0.3029E+00 -0.4526E+00 +0.5737E+01 +0.6237E+01 â€”JF +0.9608E+00 +0.5126E+01 +0.5580E+01 -0.1107E+03 +0.1448E+01 +0.5592E+01 -0.8856E+01 -0.1557E+03 Table 3.6: Estimates of 0j using Fourier functions. of 0j for j = 0,1,2 and 3 for the same modal time shifts as in example 3.1 using (3.26) and (3.30) sequentially. A vector of six Fourier functions is used with the fundamental frequency OJ0 = 2n/384Ts. The results are listed in table 3.6. 3.4.2 Estimation Using Walsh Functions As mentioned earlier, one advantage of using Walsh functions is the simplicity of computing the Walsh Transform which requires only sample addition or subtrac tion. A Walsh function is a piece-wise constant function that takes only binary values; +1 or â€”1. A transition from one value to the other may occur only at the distinct time instants kA, k = 1,...,7V â€” 1, where A = T0/N where N is the dimension of f_w and T0 is the time period over which the elements of / are defined. Consider equation (3.14). Let the vector of modulating functions be the vector of Walsh functions / given by (3.18) and let the integration in (3.14) be from 0 to T0. The integration over the time interval 0 â€”* T can be divided into inte grations over subintervals of length A such that / is constant over each of these subintervals ZjW = T,U(k*~) L , ^0* = EZHr(*A-)[P'-1y(*A)-P'-1y((fc-l)A)] k=l J(k-1)A k=1 (3.31) where /^(/cA-) is the limit of Â£w{t) as t approaches kA from the left. Define z^_l Estimation of The Modal Parameters 37 and z.j-i as 4-i (3.32) P^yfO) â€¢â€¢â€¢ D^yUN-ljA) Then substituting from (3.32) into (3.31) gives ZJW = *U/-i -Sj-i) (3-33) where the mth row of \T/ is composed of samples of Walsh function of sequency (m â€” 1) sampled at the discrete time instants kA~, 1 < k < N. The matrix \P is the Walsh matrix used to compute the discrete Walsh transform of a sampled signal [Tadokoro 78]. Samples of the output derivative of order j â€” 1 required in (3.32) are not available by measurements , so it is required to obtain estimates of these samples by some means. The required samples can be obtained from (3.5) provided that there is an estimate of 0j_1 already available (3.34) where R+d = and Rd = ' LT(A) rT{NA) '(0) rr((JV - 1)A) (3.35) and r(t) is given by (3.3). Substituting for z-_x and zj_1 in (3.33) by their estimates from (3.34) gives zjW = *{R; - itf^-i- (3.36) Finally, substituting for zjF in (3.16) by its estimate ZjW from (3.36) gives hw = -{RwRw)~1Rw^{Rd ~ Rd)ft{j-i)w- (3.37) Equation (3.37) shows that it is possible to obtain estimates for Oj, j > 0, sequen tially from the available measurements since Rw, Rd and R% are formed from the input-output measurements . Estimation of The Modal Parameters 38 3 0 1 2 3 (Ljw Ejw -0.6605E+00 +0.6113E+00 -0.3022E+00 +0.3002E+01 -0.2543E+00 -0.4325E+00 +0.5715E+01 -0.1477E+02 +0.5718E+01 -0.7850E+02 +0.3490E+02 +0.5870E+01 +0.8717E+00 +0.1370E+01 +0.4750E+01 +0.5265E+01 +0.6102E+01 -0.8434E+01 -0.1009E+03 -0.1470E+03 Table 3.7: Estimates of 6j using Walsh functions. Example 3.6 : Consider the same system as in example 3.5 with the same input-output measure ments . It is desired to obtain estimates of Â£3 for j = 0, 1, 2, and 3 for the same modal time shifts considered in the mentioned example using equation (3.37). A vector of eight Walsh functions defined on a period of length T0 = 384T,, is used. The results are listed in table 3.7. In the case of noisy measurements, the instru mental variable method can be used. A continuous-time instrumental variable is needed. A sophisticated method of generating this instrumental variable using the estimated model of the system is given in [Young 70] and [Young 79]. 3.5 Estimation Using Multi-mode Inputs A good test input signal should contain enough modes to excite all the modes of the system in order to obtain well conditioned parameter estimation equations. Such an input is referred to as a persistently exciting input. Examples of persis tently exciting inputs are periodic inputs with many frequency components, e.g. square waves and pseudo-random binary sequences. On the other hand, from the computational point of view, it is desirable to keep the number of unknowns to be estimated from a single equation as low as possible. It is seen from (2.44) that there is one unknown to be estimated for each input mode. This could cause a serious problem if the input contains a large number of modes as for example the case of a square wave input which contains an infinite number of modes. Estimation of The Modal Parameters 39 To solve this problem, a prefilter that partially, or totally, annihilates the input modes will be used to filter the input-output measurements before being used in forming the estimation equations [Trivett 81]. The filtered input and output signals have a smaller number of modes than the original signals. Let u be the response of a continuous controller as given by (2.35) and (2.36). Let the matrix Ac be diagonal so that it can be partitioned into square blocks along the main diagonal. Let Ac be partitioned into two diagonal submatrices Aci and AC2, and let xc, cc and <f>c be partitioned accordingly (3.38) (3.39) (3.40) Acl 0 -Lei 0 Ac2 ZLc2 u{t) = [ cjx c T Â±c2 2Â£ci{t) Xc2{t) xc2(t + T) MT) o o <M0 *cl(*) The controller is divided into two subsystems Cy and C2 defined by the pairs (-Aci,ccl) and {Ac2,cc2) respectively such that the controller output is the sum of two components Uy and u2 generated respectively by Cy and C2\ u(t) = Uy{t) + u2(t) = cJi ^(t) + cj2 x^it). (3.41) The input component uy contains modes of Cy only and the input component u2 contains modes of C2 only. Let ncl and nc2 be the dimension of Cx and C2 respectively such that ncl-j-nc2 = nc. Both pairs (J4cI,Cc1) and (Ac2, Cj.^ are assumed to be observable hence, we can obtain two vectors 1^ and l_2 of dimension ncy and nc2 respectively whose elements are the output modal parameters of subsystems Cy and C2 respectively such that + Yllki Â£ck{t ~ Tki) = o ,fc = l,2 (3.42) where Tki, â€¢ â€¢ â€¢, 1^, k = 1,2 , are suitable modal time shifts. The vectors of filter parameters ly and /2 are given according to (2.27) by Ik â€” ~Â£ck4 ,ck Â» ^ â€” 1> 2 (3.43) where ipck is given by (2.20) with the substitution cclc for c and <f>Ck for <f>. The right-hand side of (3.43) is zero because Uy and u2 are assumed to be the free responses Estimation of The Modal Parameters 40 of Ci and C2. Equation (3.42) represents a linear filtering operation, hence the linear filters t\ and Â£2 can be defined as as =2 ,fc = l,2. (3.44) The filter Zi annihilates all the modes in while t2 annihilates all modes in Xj,2 , so if the input u is filtered through Ci, the resulting filtered input will contain none of the modes of xcl. Conversly, if u is filtered through Â£,2, the resulting filtered input will contain none of the modes of X&. Filtering u through Zi and Â£2 results in annihilating all the modes of the input Liu{t) = cj1Â£lxcl{t) +cf2Â£1xC2(i) =cJ2Jl1xc2 MM*) = giMiZdit) + Â£2Â£2Â£i3U2{t) = Q. (3.45) Example 3.7 : Consider the input u(t) = e 2t + e 'sin5< + e 'cos5Â£. This input can be regarded as the free response of a third order continuous con troller with one real pole at â€”2 and two complex conjugate poles at â€” 1 Â± j5. The partitioned controller is given by xcl(t) = [e_2t] x^it) = [e~fsin5* e_tcos5i]T cci = [l] ^ = [1 l]r Ael = [-2] Ac2 = [-1 - 1] Mr) = [<T2T] </>C2{T) = e~T cos 5r sin 5r sin 5r cos 5r Let t\ and Â£2 be two filters which filter out the real mode and the complex modes respectively. Let the modal time shifts for Â£1 and for Â£2 be given by Tn = TT/10 r21 = TT/10 T22 = 7r/5 Estimation of The Modal Parameters 41 For these modal time shifts the corresponding ip matrices are given according to (2.20) by I'd = e~T/5 A2 = g-ir/5 Hence the vectors of output modal parameters l\ and I2 are given according to (2.27) by -71-/5 U= 0 e2*/5 To filter out the real mode, Ci is used Â£iu(t) = u{t) - e*/5u(t - TT/10) = (1 + e37r/10)e-f sin 5t + (1 - e3^10)^ cos 5t. To filter out the complex modes, Z2 is used C2u{t) = u(*) + e2,r/5u(i + TT/5) = (1 + t-2Â«l*)z 2t Obviously if Â£1 and Â£2 are applied the resulting filtered output is zero. Example 3.8 : Consider a periodic input with period T0 such that u(t) = u(t â€” T0). The fun damental frequency is u0 = 2ir/T0. In general a periodic u could have an infinite number of imaginary modes Â±jkw0, k > 0, and consequently the controller that generates u could then have order infinity. Consider filtering out a single frequency ku0. The controller can be partitioned into two parts; one part that generates the frequency ku0 0 kojn Ad -ku)0 0 Cd = [1 11 <I>C\{T) cos ku/0T sin/cwpf - sin ku)0T cos ku0T Estimation of The Modal Parameters 42 and the other part that generates all other frequencies. To perform the filtering operation, a filter Zk with the modal time shifts Tx = T0/4 T2 = T0/2 is used. For these modal time shifts the corresponding ifik is given by cos kir/2 â€” sin fc7r/2 sin /CTT/2 + cos /c7r/2 cos kir cos kn For any odd integer k0, x^kâ€ž is given by where â€”m m -1 -1 m + 1 if (k0 â€” 1) is divisible by 4 â€” 1 if (k0 â€” 1) is not divisible by 4. It follows from (2.27) that for odd integers k0 the vector of output modal param eters for filters Zko that filters out frequency k0u0 is given by [0 1] regardless of the value of k0. This means that a filter Zko with the given modal time shifts and modal parameters, annihilates all the odd frequency components in u. If u0 is a periodic signal with odd symmetry, e.g. a square wave, then Â£ko completely annihilates the modes of u Ckou0{t) = u0{t) + u0{t - T/2) = 0. It is clear from the above example that although the input may have an infinite number of modes, it was possible by a simple filtering operation to annihilate all the modes of the input. Consider the continuous ARMA model (2.44). Let vector q. be partitioned according to the partitioning of x,. in (3.38) such that â€”J V jy{t) - T.) = qfxcl(t) +q* Txc2(t). i=l Applying Ci and C2 , given by (3.44), to (3.46) respectively yields V jÂ£iy(t) + JZpJiHiy{t - Ti) = qf clSe2{t) (3.46) Estimation of The Modal Parameters 43 1=1 In (3.47) the filtered input-output measurements pairs (Â£iy , Liu) and (Â£2y , Â£2^) are used instead of the direct measurements pair (y , u) such that the number of unknowns in the two cases of filtered measurements is n + nc\ and n + nc2 respectively, while the number of unknowns if the direct measurements were used is n + nci + nc2. Example 3.9 : It is desired to obtain estimates of Oj for j = 0,1,2 and 3 of the system considered before in examples 3.5 and 3.6 from the system response to the input u(t) â€” cos6i â€” sin3i. The response of the system to this input starting from zero initial state is sampled each time interval Ts = T0/384 where T0 is the period of the input component of frequency w3 = 3. The modal time shifts TX,T2 and T3 are taken such that Tx = T2 - Tx = T3 - T2 = 40TS. The input has four imaginary modes at Â±j6 and Â±j3 which means that the dimension of each of q , j > 0 , is four. To reduce this number the technique of the measurements prefilter Â£ is used. Two filters Ci and Â£2 are used to eliminate modes of frequencies 6 and 3, respectively as in (3.47) such that in each case the number of modes in the filtered input is reduced from four to two. The choice of the modal time shifts for each filter completely defines the corresponding filter parameters. The following values for the modal time shifts are assumed; Tn = T0/8 Tl2 = T0/4 T21 = T0/4 T22 = T0/2 The prefilter parameters which result from the above choice of the modal time shifts are the same for both filters (1 = [0 If L2 = [0 If. Estimation of The Modal Parameters 44 Hence, the filtered inputs Â£iu and Â£2u. are given by Â£iu(t) = (cos6*-sin3Â£) + (cos6(i - T0/4) - sin 3(t - T0/4)) = (cos 6* - sin3t) + (cos(6t - TT) - s'm(3t - TT/2)) = â€” sin 3t + cos 3t and Â£2u(t) = (cos6i -sin3t) + (cos6(Â£ - T0/2) - sin3(i - T0/2)) = (cos 6* - sin3t) + (cos(6t - 2TT) - sin(3* - 7r)) = 2cos6i. As expected, Â£iÂ« has a single frequency w3 = 3 and Â£2u has a single frequency u>6 = 6. The corresponding filtered outputs are Â£iy(0 = y(0 + y(< - T0/4) t2y{t) = y{t) + y{t~T0/2). Estimates of the unknown parameters are obtained as before by minimizing some function of the equation error defined in this case according to (3.47) by ci(t) - D jLiy{t) + J2PjitMt - Ti) - gf Â£ixc2(f) i=i n e2{t) = D''L2y{t) +J2PjiC2y{t - Tt) - q)T C2xcl{t). i=i The state vectors xcl and xc2 are chosen such that Â£iÂ£c2{t) = [sint cosÂ£]r and Â£2Xci{t) = [sin4Â£ cos4Â£]r. Either Fourier functions or Walsh functions can be used as described in subsections 3.3.1 and 3.3.2 with the exception that the filtered measurements pairs (Â£iy, Z\u) and (Â£2y, C2u) are used instead of the direct measurements pair (y, u). Two vectors of estimated parameters 0^ = [p*T <?*T]T and = [p2^ 3**\ Tâ€¢> fÂ°r eacn j, result Estimation of The Modal Parameters 45 3 0 1 2 3 -0.1390E+01 +0.1188E+01 -0.4190E+00 +0.7137E+00 +0.2994E+01 -0.1701E+01 +0.2271E+02 -0.2640E+02 +0.6527E+01 -0.3997E+02 -0.4291E+02 +0.4671E+02 -0.3362E+00 +0.4154E+00 -0.1629E+01 +0.2898E+01 +0.4215E+01 +0.6785E+01 +0.7236E+02 -0.4658E+02 Table 3.8: Estimates of 0j using Fourier functions. 3 0 1 2 3 -2 EjF -0.1387E+01 +0.1124E+01 -0.4168E+00 +0.7382E+00 +0.2968E+01 -0.1684E+01 +0.2252E+02 -0.2622E+02 +0.6424E+01 -0.4038E+02 -0.4224E+02 +0.4634E+02 ~2 -0.9154E+00 -0.8869E-01 -0.5276E+01 +0.1354E+01 +0.3682E+01 +0.2380E+02 +0.2279E+03 +0.5318E+02 Table 3.9: Estimates of 9j using Fourier functions. from the filtered measurements pairs (Ciy, Liu) and (Â£22/, Â£2^) respectively. The results of estimating 0* and 0? using Fourier functions are listed in tables 3.8 and 3.9. Notice that the results obtained in this example for the modal parameters are different from those obtained in the previous examples because the modal time shifts considered in this example are different from those considered in the previous examples. The true values of the modal parameters are listed in table 3.10. Those values are obtained from equation (2.30). The results of estimating the same vectors using Walsh functions are listed in tables 3.11 and 3.12. Example 3.10 : Consider a fourth order system with the transfer function H( \ 400s + 400 ^ ; s4 + 6s3 + 115.25s* + 2215 + 338 Estimation of The Modal Parameters 46 3 0 1 2 3 -0.1389E+01 +0.7304E+00 +0.2263E+02 -0.4022E+02 Pj +0.1126E+01 +0.2975E+01 -0.2631E+02 -0.4258E+02 -0.4178E+00 -0.1688E+01 +0.6471E+01 +0.4649E+02 Table 3.10: True values of pj. 3 0 1 2 3 -0.1388E+01 +0.1126E+01 -0.4179E+00 +0.7306E+00 +0.2975E+01 -0.1688E+01 +0.2263E+02 -0.2632E+02 +0.6474E+01 -0.4024E+02 -0.4256E+02 +0.4649E+02 -0.3434E+00 +0.4151E+00 -0.1683E+01 +0.2893E+01 +0.4137E+01 +0.6802E+01 +0.7320E+02 -0.4654E+02 Table 3.11: Estimates of 0* using Walsh functions. 3 0 1 2 3 t S-jW EJW -0.1388E+01 +0.1126E+01 +0.4179E+01 +0.7314E+00 +0.2973E+01 -0.1687E+01 +0.2262E+02 -0.2630E+02 +0.6467E+01 -0.4023E+02 -0.4254E+02 +0.4647E+02 -0.9128E+00 -0.8730E-01 -0.5325E+01 +0.1364E+01 +0.2981E+01 +0.2371E+02 +0.2261E+03 +0.5283E+02 Table 3.12: Estimates of 6? using Walsh functions. Estimation of The Modal Parameters 47 The system has two pairs of complex conjugate poles at â€”lÂ±jl.5 and â€” 2Â±jl0, and a zero at â€”1. A fourth order model is assumed for the system and it is desired to obtain the vectors of modal parameters p. for j =0, 1, 2, 3, and 4 from the input-output measurements for suitable modal time shifts. The system is excited by a square wave of period TQ = TT such that the funda mental frequency OJ0 = 2. Starting from the following initial state y(0) = 0.2 Py(0) = 0.0 P2y(0) = 0.0 D3y(0) = 0.0 the response of the system to the square wave input is sampled such that 420 sam ples are taken every time period T0, so that the sampling time Tg equals 7r/420. The values assigned to the modal time shifts are Ti = 20TS T2 = 40TS Ts = nor3 T4 = 160T3 The differences between the modal time shifts are not uniform in this example. The square wave input has odd symmetry similar to the case considered in example 3.8. Hence, a filter Â£ with modal time shifts Ti = T0/4 and T2 = T0/2, and modal parameters l\ = 0 and l2 = 1 annihilates all the input modes Â£u(t) = u(t) + u{t - T/2) = 0. The filtered output Â£y{t) = y{t)+y{t-T/2) represents a free response of the system. In this case the dimension of vector q. in the continuous ARMA model is zero so that for each j it is required to estimate only a vector p. of dimension four. Similar to the examples considered before estimates of the unknown parameters can be obtained using Fourier functions or Walsh functions. Estimates p.F obtained using Fourier functions, are listed in table 3.13 and estimates p.w using Walsh functions are listed in table 3.14. Estimation of The Modal Parameters 48 LF hw ElF EZF hp -0.1042E+01 +0.2992E+01 +0.2658E+02 -0.4063E+03 -0.9384E+03 -0.2220E-01 -0.2784E+01 -0.1545E+02 +0.4268E+03 -0.1537E+03 +0.2179E+00 +0.1236E+01 -0.8259E+01 -0.7817E+02 +0.1074E+04 -0.5770E-01 -0.9301E-01 +0.5245E+01 +0.7556E+01 -0.6100E+03 Table 3.13: Estimates of the modal parameters using Fourier functions. EnW hw hw hw EAW -0.1042E+01 +0.2990E+01 +0.2760E+02 -0.4126E+03 -0.1016E+04 -0.2200E-01 -0.2808E+01 -0.1635E+02 +0.4350E+03 -0.9821E+02 +0.2179E+00 +0.1262E+01 -0.8266E+01 -0.8111E+02 +0.1088E+04 -0.5769E-01 -0.1059E+00 +0.5330E+01 +0.8519E+01 -0.6231E+03 Table 3.14: Estimates of the modal parameters using Walsh functions. p Â£, Ei E, Â£4 -0.1048E+01 +0.2988E+01 +0.2759E+02 -0.4118E+03 -0.1017E+04 -0.2200E-01 -0.2806E+01 -0.1635E+02 +0.4343E+03 -0.9438E+02 +0.2179E+00 +0.1262E+01 -0.8255E+01 -0.8105E+02 +0.1085E+04 -0.5768E-01 -0.1059E+00 +0.5326E+01 +0.8543E+01 -0.6219E+03 Table 3.15: True values of the modal parameters. Estimation of The Modal Parameters 49 The number of Fourier functions used is eight functions with the fundamental frequency UJ0 = 2, and the number of Walsh functions used is eight functions defined over the time period 0 â€”> n. The true values of the modal parameters, computed using (2.30), are included in table 3.15 for the purpose of comparison with the estimated values. 3.6 Summary Chapter 3 is devoted to the problem of modal parameter estimation. Estimation equations, which are linear in the parameters, are formed from discrete data as well as from continuous data. Modulating functions, e.g. Fourier and Walsh, are used in the estimation process. One advantage of estimating the parameters of the continuous-time ARMA model (2.44) over directly attempting to estimate the parameters of the input-output model (1.1) is that only one stage of integration is required. This results in avoiding the introducing of the unknown initial condi tions to the estimation problem. Recursive estimation equations, (3.26),(3.30) and (3.37), are used to obtain the vectors of the unknown parameters. An instrumental variable is used in a simple case to reduce the bias that results from noisy measurements. Signal averaging is another simple way of reducing the effect of noise. In example 3.1, for example, the input could be considered as a segment of a PRBS so that data could be averaged over different segments of the sequence. In general, multi-mode inputs are required to realize persistent excitation and it is shown how modal filters can be used to reduce the order of the parameter estimation equations. Chapter 4 Parametric Identification 4.1 Introduction The motive behind the definition of the modal functions and the modal parameters in chapter 3 is to use them in system identification. Obtaining the modal parame ters of the output and the output derivatives from input-output measurements was presented in chapter 3 with numerical examples. This chapter will concentrate on obtaining parametric models of the system through the modal functions and the modal parameters. Identification of an input-output model is considered first. A method to obtain estimates of vectors a and 6 of the input-output model through the output modal function is presented in section 4.2. An important rela tion between the vectors of modal parameters and the model characteristic equation is established in section 4.3 and then used to estimate the coefficients of the char acteristic equation. Also in the same section the modal functions of the output and the output derivatives are used to obtain an estimate of vector b. In section 4.4 a simple method for estimating the system dynamic matrix from the available mea surements is described. The coefficients of the system characteristic equation are obtained from the estimated dynamic matrix through a simple recursive relation, and the system poles are obtained as the eigenvalues of that matrix. 4.2 Identification Using the Output Modal Function The Output modal function y0 is obtained from the system output by a linear filtering operation by adding weighted samples separated by the modal time shifts 50 Parametric Identification 51 Ti,..., Tn. Similarly, a filtered input uâ€ž can be denned as n u0{t) =u{t) + J>0,u(t - T{). (4.1) t=i It follows from the assumption that the model is linear and time-invariant that the filtered input-output pair (u0,y0) satisfies the same differential equation as the original input-output pair (u, y) : Pny0(t) + T,a>Pjyo(t) = Ebjpiu^t). (4.2) j'=0 j=0 Equation (4.2) means that the output modal function y0 represents a particular solution to the differential equation model when the input is the filtered input u0 given by (4.1). The particular solution y0 contains modes from the input only and contains none of the system modes while the general solution contains modes from both the system and the input. This distinct characteristic of y0 makes the output modal function very useful for system identification. Once y0 is obtained from the general system response it becomes straight for ward to obtain its derivatives analytically since it is composed of modes of the input which are known. It follows from (2.44) for j = 0, and from (2.35) that the jth derivative of y0 is given by Djy0(t)=Â£Aixc(t). (4.3) Similarly, u0 and its derivatives can be obtained in terms of the input modes as n u0{t) = cj(/ + J2Po*H-Ti))xc(t) (4.4) Â«=i where <f)c is the state transition matrix of the controller. Let Â»T = cT{I + itpÂ°iM-Ti)) (4.5) t=i so that (4.4) becomes u0{t) = v?2u{t)- (4.6) Then from (2.35) it follows that the jth derivative of u0 is given by Dju0{t)=krAixe{t). (4.7) Parametric Identification 52 Substituting for y0, u0 and their derivatives from (4.3) and (4.7) respectively into (4.2) yields Â£{K + E *iM)2U(t) = Â»T E MkW- (4-8) ;'=0 j=0 Since (4.8) is valid for all t it follows that Â£K + E *iÂ£M = E biUTA{. (4.9) J=0 j'=0 Estimates of the unknown model parameters can be obtained by minimizing the norm of the equation error vector defined as & = Afi0 + R*Q- (4-10) where a b (4.11) and Rg is an ncx2n matrix given by Ro = koAh0---A"~lTlo -ATcv...-ArlTÂ»\. (4.12) The dimension of 6 is in general 2n, hence in order to be able to obtain an estimate of 0, the number of rows of R$ should be at least 2n, i.e. the order of the controller that generates u should be at least 2n : nc > 2n. (4.13) The estimate of 0 that minimizes the norm of ee when nc > 2n is given by 0 = -(RjRe)-lRjAfqg. (4.14) Example 4.1 : Consider the third order system of example 3.1. It is desired to obtain an estimate of the unknown model parameter vector 9 from (4.14). The dimension of 0 in this case is six, therefore according to inequality (4.13) the input should contain at least six modes. An input u given by u{t) = cos 6t â€” sin 4t â€” sin 2t Parametric Identification 53 contains three frequency components and therefore it suits the purpose of esti mating 0 from (4.14). The output of the system starting from zero initial state is computed every sampling time Ts â€” 7r/420. The number of the output modal parameters that must be estimated in order to obtain the output modal function from the system output is three which cor responds to the order of the system. The dimension of qg is equal to the number of input modes which is six in this example. The elements of qg are unknown and must be estimated with the three output modal parameters, hence, the total num ber of unknown parameters is nine. In order to reduce this number of unknowns the measurements prefilter described is section 3.5 will be employed. Three different prefilters Â£1? Â£2 and Â£3 are used, each one of them annihilates two frequency components. Let Â£l5 Â£2 and Â£3 annihilate frequency pairs (6,2), (6,4) and (4,2) respectively, and let the modal time shifts of the three filters be given respectively by 1057; 2107; 357; 70TS 1207; 1707; 5ora 1007; i5or3 2007;. The corresponding filter parameters are given respectively by [1 0] [-0.4146/5 + 00 0.1413/5 + 00 0.3972/5 + 00 0.7438/5 + 00] [-0.1616/5 + 01 0.221915 + 01 -0.1656/5 + 01 0.100075 + 01] These parameters are computed using (2.27). The reason why tx is a second order filter and not a fourth order can be understood from the result of example 3.8. It has been shown in that example that a filter of modal time shifts (T/4,T/2), where T is the period of the fundamental frequency, and modal parameters [0 1], annihilates all the odd multiples of the fundamental frequency. Hence, Â£,x annihilates both the fundamental frequency w2 = 2 and w6 = 6 = 3w2. Let the filtered input-output pairs (u^y1), (u2,y2) and (u3,y3) be the result of filtering the input-output measurements pair (u,y) through Â£,x, Â£2 and Â£3 respec tively. Each filtered input-output pair (uk,yk), k = 1,..., 3 is used to obtain an estimate of of a vector 0* of dimension five composed of the output modal param-Parametric Identification 54 eter vector p of dimension three and a vector qk of dimension two. Discrete data i-0 â€” o samples are used to obtain the required estimate of 0k as in example 3.1. The time difference tk+x â€” tkl is taken equal to lOTl,. Fifteen estimation equations are formed from each of the three filtered input-output pairs and least squares estimation is A k use to obtain an estimate 0O in each case. Three sets of results, each corresponds to one filtered input-output pair are listed in table 4.1. A; 1 2 3 -0.6586E+00 -0.6597E+00 -0.6572E+00 L +0.6087E+00 +0.6120E+00 +0.6092E+00 -0.3016E+00 -0.3034E+00 -0.3030E+00 hk -0.6326E+00 -0.5275E+00 -0.1374E+01 +0.1374E+00 +0.1229E+01 -0.2457E+01 Table 4.1: Estimates of p0 and qk, k = 1,2,3 from discrete data. The listed estimates of the output modal parameters are obtained for the following modal time shifts : 2i = 40T3 Ti = 80TS ^ = 120T3 Note that the estimates of the output modal parameter vector in each case are the same, which is expected, while the estimates of vectors g*,fc=l,2 and 3, are dependent on A:. Each vector qk, A;=l,2 and 3, gives an output modal function y* that corresponds to a filtered input uk as yk0{t) =qfxck{t) , k = 1,2,3 where ^ is the state vector of the subcontroller that generates the frequency component uk given by x<.k{i) = [sinwjfei coswfci]T , k = 1,2,3. ^see example 3.1 Parametric Identification 55 The filtered input uk obtained as in (4.6) is Â«k0{t)=ZkT2Uk{t) where u k is given according to (4.5) as *f = &(I + T,PÂ«M-TJ) i=l where <f>ck is the state transition matrix of the corresponding subcontroller. Since the filtered input-output pairs satisfy the system input-output model, it follows from (4.9) that + E*M = E^kTA{k , k = 1,2,3. j=o j=0 For each k, k=l,2 and 3, the last vector equation has dimension two. Therefore, six equations are available in the six unknown model parameters. Solving these equations gives the required estimates of the unknown parameters. The result is listed in table 4.2. +0.5205E+02 a +0.3000E+02 +0.4001E+01 0 +0.5204E+02 b +0.1280E+02 -0.4759E-01 Table 4.2: Estimate of the model parameters. Example 4.2 : Consider the same estimation problem as in example 4.1. The conditions of the estimation problem are the same as in example 4.1 with the exception that noisy output samples are used instead of the noise free samples used in that example. The noise added to the output samples is normally distributed with zero mean and has standard deviation 0.01. The bias resulting from using noisy data is evident in the parameter estimates listed in table 4.3. Parametric Identification 56 +0.5671E+02 A a +0.3945E+02 +0.4170E+01 9 +0.5596E+02 b +0.1266E+02 +0.5583E-02 Table 4.3: Estimates of the model parameters from noisy samples. The noise added to the output samples affects the estimates of the output modal parameters and the elements of vector qg which are used to obtain the required estimates of the unknown model parameters. In example 3.4 an instrumental variable was used to reduce the bias in the parameter estimates. A shifted output was used as the instrumental variable. Similarly, the same instrumental variable can be used to obtain unbiased estimates of v and qk which can be used to obtain better estimates of the unknown model â€”O â€”0 parameters. The results listed in table 4.4 shows that the instrumental variable method has worked successfully in reducing the bias in the parameter estimates. +0.5253E+02 a +0.3005E+02 +0.4015E+01 9 +0.5211E+02 b +0.1282E+02 -0.4365E-01 Table 4.4: Estimates of the model parameters using instrumental variables. Parametric Identification 57 4.3 Identification Using Modal Functions and Parameters In this section the modal functions and the modal parameters of the output and the output derivatives are used to obtain estimates of the parameters of an input-output model of the system, i.e vectors a and 6. 4.3.1 Estimation of a Theorem 4.1 : The modal parameter vectors of the output and the output derivatives pg,... ,pn satisfy the system characteristic equation : Â£ + Ea^ = oT. (4.15) 3=0 Proof : It follows from (2.30) that Â£ + Â£ ^Â£ = -cT{An + E j=0 j=0 (4.16) Since the matrix A satisfies its characteristic equation then the sum in brackets on the right-hand side of (4.16) is zero which proves the theorem. Theorem 4.1 is an important result for parametric system identification. The coefficients of the characteristic equation can be obtained by solving n algebraic equations formed from the modal parameter vectors p ,... ,p as a = -P2 P 1 M where T p !-o T Â£.i (4.17) (4.18) Example 4.3 : Consider the third order system in example 3.1. It is desired to obtain an estimate of a using equation (4.17). The modal parameter vectors required in (4.17) has been obtained in examples 3.5 and 3.6 from continuous data using Fourier functions Parametric Identification 58 and Walsh functions respectively. The estimates of the modal parameter vectors resulting from each example are used to obtain an estimate of a. The results are listed in table 4.5. It is of interest to obtain the poles of the system from the estimate of a by finding the roots of the characteristic equation. The roots of the characteristic equation resulting from the estimates of a obtained by Fourier functions and by Walsh functions are listed in the same table. Using Fourier functions Using Walsh functions +0.5027E+02 +0.5229E+02 a +0.2986E+02 +0.3003E+02 +0.3923E+01 +0.4010E+01 -0.1932E+01 -0.2011E+01 -0.9953E+00 +J0.5003E+01 -0.9995E+00 +J0.5001E+01 53 -0.9953E+00 -J0.5003E+01 -0.9995E+00 -J0.5001E+01 Table 4.5: Estimates of a and the system poles using Fourier functions and Walsh functions. Example 4.4 : Consider the fourth order system of example 3.10. It is desired to obtain estimates of a and the system poles from the estimates of the modal parameter vectors obtained in that example. Two sets of estimates of the modal parameter vectors were obtained using Fourier functions and using Walsh functions. The results of estimating a from (4.17) using the results of example 3.10 are listed in table 4.6. The roots of the resulting characteristic equation in both cases are listed in the same table. 4.3.2 Estimation of b Theorem 4.2 : The modal functions of the output and the output derivatives j/i,... ,j/n satisfy nâ€”1 nâ€”1 Vn{t) + Â£ ajyj{t) = ]T biVju{t). (4.19) j=0 j=0 Parametric Identification 59 Using Fourier functions Using Walsh functions +0.3379E+03 +0.3383E+03 A a +0.2211E+03 +0.2212E+03 +0.1153E+03 +0.1154E+03 +0.5995E+01 +0.6006E+01 -0.1000E+01 +J0.1499E+01 -0.1000E+01 +j0.1500E+01 s2 -0.1000E+01 -J0.1499E+01 -0.1000E+01 +J0.1500E+01 Ss -0.1997E+01 +.7O.IOOOE+02 -0.2003E+01 +J0.1000E+02 s4 -0.1997E+01 -j0.1000E+02 -0.2003E+01 -jO.lOOOE+02 Table 4.6: Estimates of a and the system poles using Fourier and Walsh functions Proof: Equation (4.19) looks quite similar to the input-output model (1.1). Both equations have the same right-hand side. The only difference between both equations is that the output and the output derivatives in (1.1) are replaced by the corresponding modal functions. Let y{t) = \y{t-Tl)...y{t-Tn))T (4.20) then it follows from the definition of the modal functions of the output and the output derivatives that y,-(<) = t>jy{t)+P?v{t) ,o<j<n. (4.21) Hence, fcW + EwW = 0nÂ»(*) + Â£Â»(') + l>y(0'"y(') +#Â»(*))â€¢ i=o jzzo Rearranging terms yields yn(t) + J:ajyj(t) = Dny{t) + Â£ Â«yP''y(0 + \Â£ + 3=0 j=0 j=0 (4.22) (4.23) According to theorem 4.1 the term between rectangular brackets in (4.23) is zero, thus n-1 n-1 yn{t) + E <w(*) = t>nv{t) + + X ay^'yW- (4-24) Parametric Identification 60 Equation (4.19) follows directly from (4.24) and the system input-output model (1.1) hence concluding the proof of the theorem. Equation (4.19) can be used to obtain an estimate of b from the modal functions 2/i,..., j/n and an estimate of a. Substituting for yl5..., yn and for Dn~1u in (4.19) by the equivalent in terms of the input modes from (2.44) and (2.42), respectively, yields iÂ£ + EÂ«>S,-= E^Aisu{t). (4.25) 3=0 3=0 Equation (4.25) is true for all t, hence it follows that 3=0 3=0 An estimate of b can be obtained by minimizing the norm of the equation error vector which is defined as Â§i=gl- Rbb (4.27) where Rb = \ccATecc...A^~1)Tcc\ (4.28) and n-1 2t=2Â» + T,liai- (4-29) y=o Unless the number of zeros of the system is known a priori, the dimension of b is usually assumed to be n. Therefore, the dimenssion of Rb should be at least n by n. The number of rows of Rb is nc as clear from (2.34) thus the order of the controller generating u should at least be n : nc>n (4.30) The estimate of 6 that minimizes the norm of the equation error vector when the strict inequality holds, i.e. nc > n, is given by b=(RbTRb)-1RbTqt . (4.31) Example 4.5 : Consider the third order system of example 3.1. It is desired to obtain an estimate Parametric Identification 61 of b using (4.31). The system has one zero so that there is only two nonzero elements in b. However, the dimension of b will be taken equal to three to illustrate a case where the number of zeros is not known a priori. Estimation of vectors q , j=0,... ,3, were obtained in example 3.9 using Fourier functions and using Walsh functions. In each case two sets of results were obtained q1. and q\, j=l,... ,3, each corresponds to a prefiltered input-output measurements . Each of the vectors q\ and q2., jâ€”0,... ,3, has dimension two which is not adequate to estimate b which has dimension three. However, each pair of vectors q1. and q2 can be augmented to obtain a larger vector q. of dimension four which is adequate to obtain the required estimate of b. Estimates of q\ and q\, j=l,... ,3, obtained in example 3.9 and the estimate of a obtained in example 4.3 are used to form qt which is needed in (4.31) in order to obtain the required estimate of 6. The results are listed in table 4.7. Two sets of results are listed, each corresponds to using either Fourier functions or Walsh functions in obtaining the estimates of p. and q., jâ€”0,... ,3. Using Fourier functions Using Walsh functions +0.5191E+02 +0.5181E+02 6 +0.1291E+02 +0.1266E+02 -0.3879E-01 -0.6689E-01 Table 4.7: Estimates of 6 using Fourier functions and Walsh functions. 4.4 Estimate of a from the Free Response If the input to the system is zero, the vector of unknown parameters 03 in (3.2) reduces to the modal parameter vector p^. In this case, equation (3.37) that gives the estimates of 0_3 using Walsh functions reduces to where W = -{RlRwylRl-*{R-d-R+d). (4.33) Parametric Identification 62 From (4.32) it is obvious that a modal parameter vector p. can be given in terms of the output modal parameter vector po as hw = WJlw (4-34) Substituting from (4.34) into (4.15) gives (Wn + E^Wj)poW. (4.35) i=o It can easily be shown by multiplying (4.35) by Wk, fc=0,..., n â€” 1 and using (4.34) that [Wn + Â£ ajW')PT = O (4.36) j'=o where P is the estimate of P in (4.18). Since P is invertible, it follows from (4.36) that W satisfies the system characteristic equation. This means that W, which is obtained from the system free response, can be used to estimate the system characteristic equation. The following analysis aims at establishing the relation between data matrices, similar to W obtained using Walsh functions, in a more general case. Let u in (1.2) be zero such that Xj and j// are the free system state and output respectively given as Dxf{t) = Axf(t) (4.37) yj{t)=cTxf{t). (4.38It follows from (2.1) that for any time shift r, the shifted output is given by Vfif + T) = cT(f){T)xf. (4.39) Integrating (4.37) from a time instant t0 to another time instant tf yields xf(tf) -xf{t0) = A / xf(t)dt. (4.40) It is obvious that many data equations (4.40) can be obtained for different pairs of time instants (t0i,t/i), (^2,^/2), â€¢â€¢â€¢â€¢ Writing n data equations (4.40) gives the following [xf{tfl)-xf{tol)...xf{tfn)-xf{ton)\ = A[fn xf{t)dt... ft/n xf{t)dt}. (4.41) Parametric Identification 63 Let and EÂ» = [xf(tfl) - xf{tol).. .xf{tfn) - x,(*â€žâ€ž)] Tâ€ž = [/ xf{t)dt... xf{t)dt]. *tol â€¢'ton Then A can be obtained as A â€” *- n (4.42) (4.43) (4.44) It would be straightforward to obtain A from (4.44) if all the state variables were measurable. The output is seldom the full state so that (4.44) has no practical use. However, the columns of Hâ€ž and Tn can be transformed from the state space to the output space through an invertible matrix tp as follows : [ Xf(tfl) - Xf{t0i) . . .Xf{tfn) - Xf(ton) ] y/(t/1-r1)-y/(fol-T1) ... y/^/n-ro-y/^-ro _ yf{tfl - Tn) - yf{tol - Tn) ... yf{tfn - Tn) - yf{ton - Tn) _ (4.45) Zn = ipTn = cTH-Tn) [ xf(t)dt... & xf(t)dt} f&yAt-Tjdt ... ///; yf(t - Tjdt J& y,{t - Tn)dt ... ///; yf(t - Tn)dt J (4.46) A transformation of A can be obtained from the measured output through (4.45) and (4.46) based on the following theorem: Theorem 4.3 : For an observable system, the matrix A given by A = YnZ-x (4.47) where Yn and Zn are given by (4-45) and (4-46) respectively, has the same eigen values as A. Parametric Identification 64 Proof : Substituting for Yjy and Zn from (4.45) and (4.46) respectively into (4.47) yields A = VSnT_V_1- (4-48) It follows from (4.48) and (4.44) that A = 0A^_1. (4.49) Equation (4.49) shows that A is a transformation of A which maintains the eigen values of A unchanged if is invertible which completes the proof. Although theorem 4.3 does not give the matrix A directly in terms of the available measurements, it is possible to obtain a transformation of A, matrix A, from which it is possible to obtain the eigenvalues of A and its characteristic equation. The eigenvalues of A are obtained by solving for the eigenvalues of A, and the coefficients of the characteristic equation of A are obtained from the simple recursive relation [Sinha.P.K. 84, page 33] = I ~Tk ,k = l (A tin) an-k \ -i(an_Jk_1r1 + an_tr2 + --- + oâ€ž_1rfc_1 + rfc) ,k = 2,...,n 1 J where Tk is the trace of A*. As seen from (4.50), only n matrix multiplications are needed to obtain the coefficients of the characteristic equation of the matrices A and A which makes it attractive from the computational point of view especially because A is obtained from the available measurements through simple averaging and difference operations plus a single matrix inversion of an nxn matrix. The number of columns in 5n and Tn is not restricted to the order of the system n. It is possible to form matrices 5JV and Tjv with a number of columns N > n as follows : SAT = [3â€ž xf{tf{n+1)) - xf(to{n+1)). ..xf(tfN) - xf[toN)\ (4.51) and xf{t)dt... xf{t)dt] (4.52) -o(n+l) JtoN It follows from (4.41) that a least squares estimate of A obtained from Sjv and TAT is given by A = ENTTN{TNTTN)-1. (4.53) Parametric Identification 65 Consequently an estimate of A is given by A = 0A^_1. (4.54) It can easily be shown that A is given in terms of the available measurements by l = YNZTN{ZNZTN)~l (4.55) where YN and Z^ are obtained from equations similar to (4.45) and (4.46) respec tively with N replacing n. Example 4.6 : Consider the fourth order system in example 3.10. It is required to obtain the coefficients of the characteristic equation of the system and the system poles from the free response. The free response of the system is computed for the following initial state : y(0) = 0.2 Dy{0) = 10. P2y(Â°) = 0. P3y(o) = o. Samples of the system response are computed each time interval Ts = 7r/420. First, the matrices Yjv and Zjy, N = 15 > n = 4, are formed from the free response. The values of the time shifts Ti,..., T4 are taken as follows : Tx = 40T, T2 = 80Ta T3 = 100T3 T4 = 120Ta The time instant pairs (t0k,tfk), k=l,... ,15, are taken such that tfk â€” tok = 20T3 and i0(fc+i) â€” t0k = 20TS for the specified range of A;. An estimate of A is obtained as in (4.55). The coefficients of the characteristic equation are obtained recursively from (4.50), and the poles of the system are obtained as the eigenvalues of the estimate of A. The results are listed in table 4.8. Estimation of of the system characteristic equation using theorem 4.3 is based on the assumption that y/ is the unforced response of the system. However, theorem Parametric Identification 66 a System poles +0.3388E+03 +0.2216E+03 +0.1156E+03 +0.6014E+01 -0.2007E+01 +J0.1001E+02 -0.2007E+01 -JO.IOOIE+02 -0.1000E+01 +J0.1500E+01 -0.1000E+01 -J0.1500E+01 Table 4.8: Estimates of a and the system poles. 4.3 can be extended to include systems driven by the free response of continuous controllers give by (2.35) and (2.36). The system and the controller are augmented to obtain a higher order system given by V ' x(t) ' ' A bg' 0 Ac . suit) . and y(t) = [cT <?} x(t) au[t) (4.56) (4.57) The order of the augmented system is n + nc, and the output y is considered as the free response of the augmented system as clear from (4.56) which does not include a forcing term. The augmented system dynamic matrix is given by A bg O Ac (4.58) The n + nc eigen values of Aa are those of A and Ac as can be seen from the diagonal form of A and Ac. If matrix Ma is formed from the modal matrices M and Mc of A and Ac respectively as M O O M, (4.59) Since M and Mc are invertible it follows that Ma also is invertible and Ma 1 is given by M;1 = M- O O M. -i A matrix A0 given by Aa = MaAaM~l has the same eigenvalues as Aa. From (4.58) to (4.61) it follows that A* MbgM;1 O Al (4.60) (4.61) (4.62) Parametric Identification 67 where A* and A*c are the diagonal matrices of A and Ac respectively obtained as A* = MAM-1 and A*c = McAcM~l. It is obvious that Aa is upper diagonal, and since the eigenvalues of a triangular matrix are its diagonal elements [Strang 80, page 187], it follows from (4.62) that the eigenvalues of Aa are the eigenvalues of A* and A*. The characteristic equation of the augmented system and its eigenvalues are obtained as in example 4.6 assuming a system of order n + nc and a free system response. The eigenvalues of A which correspond to the poles of the actual system can be distinguished from the eigenvalues of Ac which are known. The charac teristic polynomial of the augmented system is the product of the characteristic polynomials of the actual system and the characteristic polynomial of the con troller. So if P(s), Pc{s) and Pa{s) are the characteristic polynomials of the actual system, the controller and the augmented system respectively then The characteristic polynomial of the actual system can be obtained from (4.63) either by equating the coefficients of the like powers of 5 on both sides of the equa tion or by long division, under the assumption that the characteristic polynomial of the controller Pc(s) is known. Example 4.7 : The third order system of example 3.1 is excited by a sinusoidal input of fre quency 4 The order of the controller is nc = 2 and has two imaginary poles at Â±j4. The output of the system is computed for zero initial state with the sampling time TS â€” 7r/420. An augmented system of order n + nc = 3 + 2 = 5 is assumed. The matrices YN and ZN, N = 20 > 5, are formed from the system response for the following time shifts : Pa(s) = P(s)Pc(s) (4.63) u(t) = sin 4t . Xi = T2 = TS = 110TS 4or3 8ors Parametric Identification 68 T4 = HOT/, T5 = 170Ta The time instant pairs (t0k,tfk), k=l,... ,20, are taken such that tfk â€” t0k = 20T3 and t0(k+i) â€” t0k â€” 20TS for the specified range of k. An estimate of A0 is obtained as in (4.55). The coefficients of the characteristic equation of the augmented system are obtained from (4.50) using the estimate of Aa. The estimated characteristic polynomial is Pa(s) = s5 + 3.992s4 + 45.93s3 + 115.5s2 + 478.4s + 827.1. The characteristic polynomial of the controller is Pc(s) = s2 + 16 . Substituting in (4.63) and equating the like powers of s on both sides gives the following equations in the unknowns a0, ax, and a2 16a0 = 827.1 16a! = 474.4 16a2 + a0 = 115.5 a2 = 3.992 These five equations can in general be solved by least squares to obtain the three unknowns. However a more direct solution can be obtained from the first, the second and the last equations. The results are listed in table 4.9. A a K Eigenvalues of A0 +0.5169E+02 +0.2965E+02 +0.3992E+01 +0.3146E-02 +J0.3995E+01 +0.3146E-02 -J0.3995E+01 -0.1003E+01 +j0.5000E+01 -0.1003E+01 -J0.5000E+01 -0.1993E+01 +J0.0000E+00 Table 4.9: Eigen values of Aa. The poles of the augmented system are obtained by solving for the eigenvalues of the estimate of Aa. The results are listed in table 4.9. It is obvious that the first Parametric Identification 69 two eigenvalues correspond to the controller two imaginary poles while the other three correspond to the actual system poles. The new method for obtaining the poles of the system as the eigenvalues of an output data matrix, as described above, has important advantages over the known methods based on discrete prediction models of the system like Prony's method. The new method uses continuous measurements so that possible aliasing, associated with sampling continuous signals, is not a problem. The flexibility of choosing the modal time shift allows for modeling the signal low frequency components as well as the high frequency components without aliasing. The poles of the continuous system are obtained as the eigenvalues of a data matrix, while if a discrete model is used, the continuous system poles are obtained from the roots of the characteristic equation of the estimated discrete model by a logarithmic operation. The following example illustrates a case of a wideband signal where Prony's method fails to determine the exponential modes of the signal accurately, while the new method yields accurate results. Example 4.8 : Consider a signal y(t) = e~f sin3f + cos 100*. Let u>i and w/, be 3 and 100 respectively, so that the period of the low frequency component and the high frequency component are TJ = 27r/3 and Th = 27r/100 respectively. Samples of y are available every sampling interval Ts = Ti/666 Â« T/i/20. It is required to determine the complex numbers sk, k â€” 1,... ,4, so that y can be modeled by y(t) = tcÂ»eSkt-It is clear that the complex exponents are Sl = -1 + j3 s2 = -1- j3 53 = +.7100 s4 = -jlOO. First, estimating the complex exponents is done indirectly through the discrete Parametric Identification 70 prediction model y(mT) + ^M((m-fc)T)-0. jfe=i The parameters of the discrete model are estimated from samples of y using least squares p= -{RTRy1RTy where R is an nxN matrix whose mth row is given by ilR= [y(m + 2) ... y(m - 1) ] and y = [ y(4) ... y{N + 3)]T . The poles of the discrete model zk, k â€” 1,..., 4, are obtained as the roots of 4 Finally each sk is obtained from zk by the logarithmic operation (1.6). The choice of T in the discrete model is a problem because y has two widely separated frequency components. Choosing T comparable with TJ results in aliasing as seen from the results listed in table 4.10, which are obtained for T = 33TS and TV = 10. On the other hand, choosing a small T according to the sampling theory, T < Th/2, results in losing the accuracy of modeling the low frequency component as seen from the results listed in the same table , which are obtained for T = 4TS Â« 2ft/5 and the same N. It is seen from the results in the table that the accuracy of the estimate of w( is deteriorated by the small T. The condition of the estimation equations for the discrete model parameters can be measured by testing the condition number, [Strang 80, page 282], of the matrix (RTR). The higher the condition number is, the more ill-conditioned the matrix is. The condition numbers of the two cases of T = 33TS and T = 4TS are 0.6768E+01 and 0.6816E+04, respectively. The condition number in the case of small T is three orders of magnitude larger than the condition number in the case of large T. This can be interpreted as due to the high correlation between the columns of R in the case of small T. Parametric Identification 71 Using T = 33TS Using T = 4TS -0.1000E+01 +J0.3000E+01 -0.1000E+01 -J0.3000E+01 -0.5744E-06 +j0.2109E+02 -0.5744E-06 -j0.2109E+02 -0.1187E+01 +J0.2261E+01 -0.1187E+01 -J0.2261E+01 -0.5260E-03 +J0.1000E+03 -0.5260E-03 -JO.IOOOE+03 Table 4.10: Estimates of sk from the discrete model. Applying the new method and using the modal time shifts T = %T where T = 33Ta, and using time instant pairs [tok,tfk), k=l,... ,10, such that tfk â€” tok = to(k+i) â€” tok = 33Ta, give the results listed in table 4.11. It is obvious that there is no probem of aliasing despite the large modal time shifts. Using T = 33TS -0.1000E+01 +J0.3000E+01 -0.1000E+01 -J0.3000E+01 +0.3586E-03 +j0.1008E+03 +0.3586E-03 -j0.1008E+03 Table 4.11: Estimates of sk using the new method. 4.5 Summary Chapter 4 deals with the use of modal functions and parameters for determining model parameter estimation equations that do not depend on the system initial conditions or the use of repeated integrations to eliminate derivatives. The use of the output modal function, which is composed of known exponential modes, in place of the general system output gives an algebraic vector equation (4.9) that relates the exponential modes of the output modal function to those of the filtered input through the unknown model parameters and which can be solved for these parameters. Theorem 4.1 shows that the modal parameter vectors satisfy the system charac teristic equation. Using this result yields (4.17) which proves that the a-parameters can be identified independently of the 6-parameters. When the a-parameters are Parametric Identification 72 known, theorem 4.2 can be used to identify the 6-parameters as in (4.26). Theorem 4.3 proves that the system eigen values can be identified from an output data matrix. By use of an augmented state vector this approach can be extended to the case of nonzero inputs. The output data matrix allows the identi fication of the coefficients of the system characteristic polynomial recursively from (4.50). It is shown that the new approach yields accurate results with wideband signals whereas Prony's method fails. Chapter 5 Nonparametric Identification 5.1 Introduction The identification of nonparametric system models is considered in this chapter. Samples of infinite dimensional models, namely the frequency response, the step response and the impulse response, are obtained. Obtaining the system frequency response by sinusoidal input testing usually requires long experimental time to let the system response reach the periodic steady state, and therefore it is not suitable for real-time applications. On the other hand, obtaining the frequency response from finite data records has its own drawbacks. Data windowing results in a smearing of the signal spectrum in a manner dependent upon the precise shape of the window, e.g. rectangular, Hanning or Hamming [Marple 87, Sec. 5.3]. It is necessary to average the spectrum estimates over a large number of data windows in order to reduce the variance induced in these estimates by the response transient component. Finite data frames also induce bias in the spectrum estimates which is not decreased by increasing the number of frames. The bias is reduced by increasing the duration of the data frame, i.e. longer test time. Estimating the impulse response by deconvolution using the input auto correlation and the input-output cross-correlation functions involves windowing, bias and variance errors similar to the case of estimating the frequency response [Wellstead 81]. An advantage of direct nonparametric modeling over parametric modeling of systems is that no a priori knowledge of the system order is required. However, the indirect nonparametric methods described in this chapter do not have this advantage because they are based on the definitions of the modal functions 73 Nonparametric Identification 74 and the modal parameters which require the knowledge of the system order. Section 5.2 presents a method of estimating the frequency response of the sys tem from the output modal function for a periodic waveform input. The output modal function is utilized again in section 5.3 to estimate the step response from the system response to a general piecewise constant input. The system impulse response is obtained as the derivative of the estimated step response using the output first derivative modal parameters in section 5.4. 5.2 Estimation of the Frequency Response The output modal function y0 played a key role in estimating the unknown model parameters in section 4.2. As mentioned earlier the output modal function repre sents a particular solution to the model differential equation when the input is the filtered input u0 given by (4.1). This particular solution has the property that it contains modes of the input only and the transient component of the system output due to system modes is not present in y0. Consequently, if the input is periodic, i.e. a summation of complex conjugate imaginary modes, the output modal function is periodic with the same spectral components as the input. The periodic filtered input-output pair (u0,y0) can be used to estimate samples of the system frequency response. Let u be a sinusoidal input with frequency u>0 such that the controller matrix Ac in phase variable form is given by and the controller state vector is given by x,.(t) = [s'mu>0t costj0t]T (5.2) The filtered input u0 and the filtered output y0 are given by (4.6) and (2.26), respectively. Choosing the controller state equation in phase variable form makes it possible to obtain vectors v_ and go, in (4.6) and (2.44), by correlating u0 and y0, respectively with the vector [sino;0i cosw0]r over the period of the input. (5.1) Nonparametric Identification 75 It follows from equation (4.9), by multiplying both sides by the the inverse of U? + Â£-=o1MD,that nâ€”1 nâ€”1 3=0 3=0 (5.3) The diagonal form of AC in (5.1) is given by A: = JU0 0 0 -jw0 (5.4) The matrix Ac is related to its diagonal form A* through the modal matrix M as Ar = MAIM where and M M'1 = -Substituting from (5.5) into (5.3) gives â€”J J 1 1 J 1 -J 1 J'=0 (5.5) (5.6) (5.7) (5.8) J=0 The frequency response of a system is given as a complex function of the frequency u> by #(jcu) Hence, it follows from (5.9) and (5.8) that qT = vTM H{ju0) 0 0 H{-ju0) M -l (5.9) (5.10) Substituting for M and M 1 in (5.10) gives T T J ^ â€” # (jw0) + F(-ju;0) -j(ff(jw0) - H(-JOJ0)) J(H{JUJ0) - H{-ju0)) H{JUJ0) + H(-ju0) (5.11) Let Re(H{](jj)) and Im{H{ju>)) be the real part and the imaginary part of H(JOJ), respectively. It follows from (5.11) that Â£ = Â»T Re(H(jcj0)) Im{H{ju0)) -Im(H{juj0)) Re(H(ju0)) (5.12) Nonparametric Identification 76 From matrix and vector multiplication rules, it can easily be verified that (5.12) can be written as â€”0 where V\ and u2 are the elements of vector u_. A sample of the system frequency response evaluated at frequency UJ0 is obtained from (5.13) as ' Re{H{3u0)) Im(H(ju0)) } = Â£ Equation (5.14) gives one sample of the required frequency response. In order to obtain samples at different frequencies, a periodic input with multiple frequencies should be used. For each harmonic frequency uk = koj0, two vectors uk and qk, each of dimension two, can be obtained by correlating u0 and y0, respectively with a vector [sinwji coswkt]T over a time period of the fundamental T0 = 27r/w0. Each resulting vector pair i/* and q^ can be used in (5.14) to obtain one sample of the frequency response at frequency u>k. Example 5.1 : From the same input-output measurements of the system considered in exam ple 3.1 it is required to obtain samples of the frequency response of the system. Starting from rest, the system is excited by a square wave with a fundamental fre quency u0=l. Only odd harmonics of u>0 are present in the spectrum of the input, so it is possible to estimate only the frequency response samples of the system at UJ0 and odd multiples of u0. The output modal parameters were obtained in example 3.1. these parameters are used to filter the input-output measurements (u,y) to obtain {u0,y0). Pairs of vectors uk and qk, for odd k, are obtained by correlating u0 and y0, respectively, with a vector [sin/cf cos kt\T over the time period of the fundamental which is equal to 2TT in this case. Vector pairs obtained for odd A; from fc = ltoA; = ll are listed in table 5.1. These vectors are used in (5.14) to obtain the required estimates of the system frequency response samples at the corresponding frequencies. The results are listed in table 5.2 with the true values of the corresponding samples as computed from (5.9). Re(H{juj0)) Im{H(ju0)) -u2 vx (5.13) Uy U2 -u2 ux (5.14) Nonparametric Identification 77 k uk ak %i 9*2 1 +0.1436E+01 +0.2359E+01 +0.1974E+01 +0.1751E+01 3 -0.9008E+00 +0.3559E+00 -0.4778E+00 +0.8423E+00 5 +0.1395E+00 -0.3717E+00 -0.6021E+00 -0.1003E+00 7 +0.4124E+00 +0.6293E+00 -0.1089E+00 -0.3875E+00 9 -0.6618E+00 +0.7798E+00 +0.1795E+00 -0.1675E+00 11 -0.8865E+00 -0.2632E+00 +0.1213E+00 +0.4375E-01 Table 5.1: Spectral components of u0 and y0. k Re{H{juk)) Im(H(juk)) Re{H{juk)) Im(H(juk)) 1 +0.9135E+00 -0.2811E+00 +0.9135E+00 -0.2811E+00 3 +0.7784E+00 -0.6276E+00 +0.7785E+00 -0.6277E+00 5 -0.2964E+00 -0.1509E+01 -0.2974E+00 -0.1509E+01 7 -0.5101E+00 -0.1612E+00 -0.5099E+00 -0.1610E+00 9 -0.2384E+00 -0.2782E-01 -0.2383E+00 -0.2795E-01 11 -0.1392E+00 -0.8015E-02 -0.1393E+00 -0.8181E-02 Table 5.2: Estimated and true values of samples of the frequency response. One advantage of the method described above to obtain the frequency response of the system is that it uses input-output measurements over a relatively short time interval which makes the method suitable for unstable systems as well as for stable systems. The periodic response of the system is obtained by filtering u and y to get u0 and y0 through the modal filter which annihilates the modes contributed by system which represents the transient response of the system. Practically, to obtain a periodic response the experiment time should be long enough to let the transient component of the response settle down which is not possible for an un stable system. If on the other hand digital spectral estimates of the auto and cross correlation functions of u and y from finite data records are used to obtain the required estimates of the system frequency response, the obtained estimates will suffer from the inherent errors to the digital spectral estimation method, namely windowing errors and aliasing [Wellstead 81]. Nonparametric Identification 78 5.3 Estimation of the Step Response A simple input that can be applied to a system for the purpose of identification is a step input. If the system is initially at rest when a step input is applied, then the output of the system is the step response. If this is not the case and the system starts from a nonzero state, which is common in practice, the output of the system is the sum of the step response and the response due to the nonzero initial state. A more general input is a piece-wise constant input given by k MO = J2 v{m)A(t - mT) ,kT<t<(k + l)T (5.15) tn=0 where , , f u(mT) - u((m - 1)T) ,m>0 ,r . v m) = \ ) ' u > > ' 5.16) v ' I u\m) , m = 0 K ' and A(t) is the unit step function. The impulse response h of the system is given by MO = ^H^h- (5.17) Substituting from (5.17) into (2.1) using the piece-wise constant input given by (5.15) yields 2/(0 = cT<j>{t)x(0) + h(t-T)Y^, v(m)A(T - mT)dT. (5.18) JÂ° m=0 The upper limit of the summation in (5.18) is A; although r does not necessarily satisfy kT < T < (k + 1)T as required in (5.15). However, it is possible to use A: as the upper limit since the step function A(t) is zero for t < 0. Interchanging the order of integration and summation in (5.18) and using the property that A(t) is zero for t < 0 yields 2/(0 = cT<f>{t)x{0) + E /' ~ T)dTV(m)- (5-19) m=0 JmT The step response g is the integration of the impulse response h: g{t) = j h(r)dT (5.20) consequently (5.19) can be written in terms of g as k y{t) = cT(j){t)x{0) + Es(*- mT)v(m). (5.21) m=0 Nonparametric Identification 79 The fact that g(0) is zero1 is used to obtain (5.21) from (5.19). It is clear from (5.21) that y is composed of the response due to the nonzero initial state and the convolution summation of the input steps and the system step response. The determination of the step response from (5.21) by deconvolution is com plicated by the existence of the initial state response term. The effect of this term is usually accounted for by considering N terms in the convolution summation with iV large enough so that the effect of the nonzero initial state is negligible for t > NT. Therefore, for a stable system with an impulse response that vanishes for t > NT, the output is given approximately by k y(t) Â« g(t - {k - N)T)u(k - N) + E g{t-mT)v{m). (5.22) m=k-N+\ An alternative exact approach to the deconvolution problem is considered below using the output modal function. The output modal function as defined by equation (2.26) is the sum of weighted output samples. Substituting for the output samples in (2.26) from (5.21) and taking the modal time shifts T,=tT, iâ€”,... ,n, yields k n kâ€”i Vo{t) = Â£ <7(* ~ mT)v{m) + Â£p0,- Â£ d{t - mT - iT)v{m). (5.23) m=0 t=l m=0 The response component due to the nonzero initial state is annihilated as a result the output filtering operation, so that it is possible to obtain the step response by deconvolution using y0 as given below. Each of the summation from m=0 to m=k-i, i=0,... ,n, can be divided into two summations; one from m=0 to m=k-n and the other from m=k-n+l to m=k-i such that (5.23) becomes, after arranging terms, kâ€”n n Vo{t) = E \g{t - mT) + Â£p0,<7(* -mT- iT)]v{m) m=0 j=l k nâ€”1 kâ€”i + E 9{t-mT)v{m) + J2Poi E g{t - mT - iT)v{m\h.2A) mâ€”kâ€”n+l t=l mâ€”kâ€”n+l The term in rectangular brackets in (5.24) is the output modal function of the system which corresponds to a step input. It has been shown before that the follows from equation (5.20) Nonparametric Identification 80 modal function which corresponds to a constant input is a constant, hence for t > nT the term in the rectangular brackets in (5.24) is constant and can be taken outside the summation n kâ€”n y0{t) = \g{t-{k-n)T) + Y,Poi9{t-{k-n)T-iT)]J2vn t'=l m=0 k nâ€”l kâ€”i + Â£ g{t-mT)v{m) + Y,Poi Â£ g{t-mT -iT)v{m). m=kâ€”n-f-1 Â«=1 m=kâ€”n+1 Substituting t = kT + s ,0<s<T (5.25) and augmenting the last two summations over i and changing the summation ar gument m to l=k-rn-i, i=0,... ,n-l, gives y0{kT + s) = [g{s + nT) + Y,Poig{s + {n-i)T)\u{k-n) x=i + J^Poi^gis + lT^ik-i-l). In (5.26) p00 = 1. Equation (5.26) can be put in vector form as y0{kT + s) = [v{k) v(k - 1) ... v(k - n + 1) u(k - n) ] 1 O Pol 1 L Pon â€¢â€¢â€¢ Pol 1 g(s + T) g{s + nT) . Let the row vector uT be (5.26) (5.27) uT(k) = [ u(k) u{k-l) ... u(fc-n)] (5.28) such that v(k) v(k-l) ... u(fc-n + l) u(A; - n) ] = uT(/e)A (5.29) where A = 1 -1 1 0 -1 O o 1 0 -1 1 (5.30) Nonparametric Identification 81 Let g(s) be a vector of dimension n+1 given by g(s)= g{s) g{s + T) ... g{s + nT) lT (5.31) A least squares estimate of g(s) can be obtained from N > n + 1 equations like (5.27) for different k as where g(s) = {U U)~ U y_g{s) yo(s) = [y0{s) y0{S + T) ... y0{s + {N - 1)T) \ 1 O L = and the kth row of U is given by Poi 1 Pon â€¢â€¢â€¢ Pol 1 lTkU = u(k) . (5.32) (5.33) (5.34) (5.35) Equation (5.32) gives estimates of samples of g(t) over the time interval 0 â€”> (n + l)T only. Estimates of samples of g(t) beyond this time interval can be obtained by extrapolation as given below. For t > nT > 0 the step function is constant and consequently the output modal function of the step response is a constant given by g{t) + T,Poig{t-iT) =q00 ,t>nT>0. (5.36) t=i The extrapolation relation of g(t) for t > s + (n + k)T, k > 1 follows from (5.36) as n g(s + {n + k)T) = q00-J2Poig{s + {n + k-i)T) i=i n = g{s + {n + k - l)T) + J2Poi9{s + {n + k - i - 1)T) t=i n - X>oi0(s + (^ + â„¢ ~ 0R) ,fc>l- (5.37) i=l Example 5.2 : It is required to obtain estimates of samples of the step response of the third Nonparametric Identification 82 order system in example 3.1 using (5.32) and (5.37). The system is excited by a piece-wise constant input +1, +2, -1 and -2 such that each constant input level lasts for 80rs where Ts = 7r/320. Starting from rest, the output of the system is computed every sampling time Tâ€ž. The time interval T is taken equal to 20TS. The output modal parameters are obtained as in example 3.1 from discrete data, and (5.32) is used to obtain estimates of g(t) over the time interval 0 â€”> 80TS. Finally, (5.37) is used to obtain estimates of samples of g(t) for t > 80TS. Estimates of samples of g(t) obtained every 10TS are listed in table 5.3. The true values of the corresponding samples are listed beside the estimated values. From the listed results it is seen that the estimated values are exactly the true values represented by four decimal digits. 5.4 Estimation of the Impulse Response Theoretically, the impulse response is the system response to a unit impulse input. With some exceptions, such as biological and biomedical experiments where impul sive test signals are often the most practical forms, using an impulsive test signal to estimate the system impulse response is not practical. The impulse response of a system is usually estimated by a process of deconvolution knowing the input auto-correlation and the input-output cross-correlation functions. In this section, a method is considered for obtaining the impulse response as the derivative of the step response using the output derivative modal function. The system impulse response is related to the step response by The output derivative modal function that corresponds to a step input is obtained from (2.44) by setting the input single mode to unity as h{t) = Vg{t). (5.38) n Dg{t) + T,Pu9{t-iT)=ql0. (5.39) Â«'=i From (5.38) and (5.39) it follows that n h{t) + J2Pug{t-iT) = ql0. (5.40) Nonparametric Identification k Estimated values True values 0 +0.6047E-06 +0.8327E-16 1 +0.2291E+00 +0.2291E+00 2 +0.2291E+00 +0.2291E+00 3 +0.1086E+01 +0.1086E+01 4 +0.1154E+01 +0.1154E+01 5 +0.1004E+01 +0.1004E+01 6 +0.8457E+00 +0.8457E+00 7 +0.8207E+00 +0.8207E+00 8 +0.9160E+00 +0.9160E+00 9 +0.1028E+01 +0.1028E+01 10 +0.1071E+01 +0.1071E+01 11 +0.1039E+01 +0.1039E+01 12 +0.9821E+00 +0.9821E+00 13 +0.9542E+00 +0.9542E+00 14 +0.9680E+00 +0.9680E+00 15 +0.1000E+01 +0.1000E+01 16 +0.1021E+01 +0.1021E+01 17 +0.1018E+01 +0.1018E+01 18 +0.1002E+01 +0.1002E+01 19 +0.9804E+00 +0.9894E+00 20 +0.9887E+00 +0.9887E+00 21 +0.9967E+00 +0.9967E+00 22 +0.1005E+01 +0.1005E+01 23 +0.1006E+01 +0.1006E+01 24 +0.1003E+01 +0.1003E+01 25 +0.9981E+00 +0.9981E+00 26 +0.9965E+00 +0.9965E+00 27 +0.9980E+00 +0.9980E+00 28 +0.1001E+01 +0.1001E+01 29 +0.1002E+01 +0.1002E+01 30 +0.1001E+01 +0.1001E+01 31 +0.9999E+00 +0.9999E+00 Table 5.3: Samples of g(t); estimated and true val Nonparametric Identification 84 Once the samples of g are estimated they can be used in estimating samples of h using (5.40) provided that there are available estimates of the output derivative modal parameter vector pt and the scalar qlo. It is seen that (5.40) can be used for estimating samples of h only where t > nT. It is possible, however, to obtain samples of h for t > 0 if negative modal time shifts are used. First, the estimates of p1 and ql0 are obtained for negative modal time shifts and then these estimates are used in (5.40) to obtain h(t) for t > 0 in terms of g(t) for t > â€”T. Estimation of the output derivative modal parameters has been considered in chapter 3. The estimation process can be made easier by making use of the step response, already estimated in example 5.2, instead of the general system response. Multiplying (5.40) by a vector of modulating functions f(t), e.g. Fourier or Walsh functions, then integrating over the period of orthogonality T0, making use of the fact that g(0) = 0, gives - / " g{t)Df{t)dt + g{T0)f{T0) + Y.Pii / 9{t-iT)f{t)dt = qlo / Â° f{t)dt. (5.41) Jo â€” â€” i=l Jo ~ Jo ~ Least squares estimates of p , qlo and g{T0) can be obtained from (5.41) using a vector / of dimension N > n + 2. It should be noted that employing a vector of N Fourier functions as given by (3.17) is not suitable for estimating the unknown parameter q\0 because the integral of / over the period of orthogonality is the zero vector. However, it is possible to estimate all other unknowns in (4.44) so that it becomes possible to estimate qio by correlating (5.39) with a unit function f0(t) = 1 over the time interval t â€”> T0 as follows: qio = ^r(ff(T0) + X>Â« / Â° g{t - iT)dt). (5.42) Example 5.3 : It is required to obtain estimates of the impulse response samples of the system considered in example 5.2 from the estimates of the step response samples obtained in that example. Assuming the following modal time shifts Tx = -20TS r2 = -4ors T3 = -60TS and using a vector of six Fourier functions with a fundamental time period T = Nonparametric Identification 85 192Ta in (5.41) result in the following estimates of the output derivative modal parameters: pu=-0.1678E+01 pi2= 0.8164E+01 pi3= 0.1126E+02 The estimate of qi0, as obtained from (5.42) is qlo= 0.1774E+02. Samples of the impulse response are estimated from (5.40) using these obtained estimates. The results are listed in table 5.4. The true values of the corresponding samples are also listed beside the estimated values. The listed samples are obtained every 10T3. 5.5 Summary Methods of identifying nonparametric models are described in chapter 5. Equation (5.14) gives the system frequency response at frequency u0 in terms of vectors qo and v obtained from the input-output measurements. Such a simple direct relationship cannot be obtained when conventional frequency response methods, such as digital spectral analysis with the Discrete Fourier Transform, are used in the general case with output transients due to the system poles. The output transient component is not present when the output modal function is used. The step response is obtained from the output modal function in the case of a piece-wise constant input through (5.32) and (5.36). The impulse response is obtained by a simple filtering operation of the estimated step response, (5.40). Nonparametric Identification k Estimated values True values 0 +0.3080E-03 -0.9992E-15 1 +0.2135E+01 +0.2132E+01 2 +0.2464E+01 +0.2461E+01 3 +0.1156E+01 +0.1155E+01 4 -0.3732E+00 -0.3724E+00 5 -0.9625E+00 -0.9616E+00 6 -0.5194E+00 -0.5195E+00 7 +0.2483E+00 +0.2472E+00 8 +0.6266E+00 +0.6253E+00 9 +0.4378E+00 +0.4369E+00 10 -0.6800E-03 -0.1019E-02 11 -0.2808E+00 -0.2809E+00 12 -0.2454E+00 -0.2457E+00 13 -0.2713E-01 -0.2778E-01 14 +0.1460E+00 +0.1452E+00 15 +0.1551E+00 +0.1544E+00 16 +0.4564E-01 +0.4508E-01 17 -0.6123E-01 -0.6168E-01 18 -0.8530E-01 -0.8580E-01 19 -0.3548E-01 -0.3604E-01 20 +0.2612E-01 +0.2547E-01 21 +0.4851E-01 +0.4788E-01 22 +0.2734E-01 +0.2672E-01 23 -0.7232E-02 -0.7797E-02 24 -0.2450E-01 -0.2507E-01 25 -0.1698E-01 -0.1757E-01 Table 5.4: Samples of h(t); estimated and true values. Chapter 6 Observer Design 6.1 Introduction The design of a state observer that gives the internal state of the system in terms of the input-output measurements in case of piece-wise constant inputs is considered in this chapter. It is assumed that there is no available a priori knowledge of the internal structure of the system so that the choice of the state variables is free. The most natural choice of the state variables in this case is the phase variables, which have a close relation with the system output and its derivatives. The state modal functions defined in chapter 2 will form the basis for the design procedure. In section 6.2, a design procedure for a continuous-time and a discrete-time observers based on the definition of the state modal functions is given. The material given in section 6.3 completes the identification of a state model of the system by introducing a method for estimating the vector 6. 6.2 State Modal Functions as a State Observer Each of the state modal functions yj(t), j = 0,... ,n â€” 1 as given by definition 2.2 is a weighted sum of the state variable Xj+i at time t and n earlier samples of the output y. On the other hand, the same state modal function can be derived from equation (2.17), which gives a modal function fa for arbitrary vector a, by stting a = cTAj such that aTx(t) = cTAjx(t) = cTAjx(t) â€” ij+lx(t) is the state variable These specific choices of a have led to the definition of the state modal parameter vectors p.,j = l,...,n â€” l, given by (2.30). The state modal 87 Observer Design 88 parameter vectors are identical to their respective output and output derivative modal parameter vectors which can be estimated using the methods described in chapter 3. In the following analysis wherever the index j appears, it is assumed that j can take the values from 0 to n-1. Let u be a piece-wise constant input as given by (5.15). Let the modal time shifts be integral multiples of the input sampling interval T such that Ti â€” iT, i = 1,..., n. It follows from (2.17) by substituting pai,... ,pan for the state modal parameters pji,... ,pjn that the state modal function y; is given at the time instant t = kT + s by n r(k-i)T+s y3{kT + s) = Y^Pji / cT(/>({k - i)T + s - r)6u(r)aY. (6.1) Let r(k-i)T+a Zi(kT + s)= cT<j>((k - i)T + s - r)bu(r)dT (6.2) JkT+a such that y,-(fcr + s) = ^p^T + s). (6.3) Â»=i The analysis given below aims at expressing 2, in terms of the input steps in order to obtain a relation between the state modal function y3 and the piece-wise constant input u. Substituting for u in (6.2) as a sum of steps yields f(k-i)T+3 Zi(kT + s) = / cT<f){{k - i)T + s - T)bdru(k-i) JkT+s * r(k-i+m)T + Â£/ cT<j){{k-i)T + s-T)bdTv(k-i + m). (6.4) m=l JkT+<> Using the substitution o â€” (k â€” i)T + 5 â€” r in the integrals gives rÂ° Zi{kT + s) = / cT<j>(o-)bdo-u(k -1) * rs-mT + Â£ / cT4>(o-)bdo-v{k-i + m). (6.5) m=l J-iT Each integral over the time intervals t = â€”iT â€”> s â€” mT in (6.5) can be expressed as sum of two integrals over the time intervals t = â€” iT â€”> 0 and t = 0 â€”> s â€” mT Observer Design 89 such that (6.5) becomes Zi(kT + s) = / cT(j)(a)bdou(k - i) + V / J <f>{o)bdov{k - i + m) J~iT m=lJ~iT * rs-mT + Yl cT^{a)bdov{k-i + m). (6.6) m=l^Â° Since from the definition of v i Â£ v(k - i + m) = u(/c) - u(fc - i) (6.7) m=l it follows from (6.6) by adding the first two terms on the right-hand side of the equation that /0 * rs-mT cT(j){a)bdau{k) + V] / cJ'4>{o)bdov{k -i+m). (6.8) -iT m=iJo Substituting back for z, in (6.3) yields Â« rO n i rs-mT Vj{kT + s) = YlPji / cT<t>{o-)bd(Ju(k) + J^Pi* / cT<j)(a)bdav(k -i + m). i=l ~iT t'=l m=l 0 (6.9) Introducing rm{s) = / cT<j>{c)bdo- (6.10) Js-mT gives n n i yj(kT + s) = Ep,-,-rt-(0)u(fc) - J^Pji Â£ rm{s)v{k - i + m). (6.11) t=l m=l Let yj0 be the state modal function that corresponds to a constant input u(fc), then Vjo{kT + s), k > n, can be obtained from (6.11) by setting all the terms v(k â€” i + m) to zero as follows: Vj{kT + s) = X>;.r,(0)u(fc)- (6-12) t=i On the other hand, yJ0 is given by (2.40) as Vjo = Â§,â€¢<,Â«(*;) (6.13) where qJ0 is obtained from (2.41) by setting = <pc(t) = 1 as Qjo = J2Pji cT(f>(a)bd(T. (6-14) .=i y-'T Observer Design Then it follows from (6.12) and (6.13) that n i=l Substituting from (6.15) into (6.11) yields n i y,(kT + s) = qjou{k) - Â£p>,- Â£ rm(s)v(k - i + m). t=l m=l Equation (6.16) can be written in vector form as Vj{kT + s) = [u(k) -v{k) -v{k-l) ...-v{k-n + l)]Tj where and 90 (6.15) r(s) Â»â€¢â€ž(Â«)] ' 1 0 ... 0 0 pji â€¢â€¢â€¢ Pjn . 0 Pjn o be shown that \u{k) -v(k). . â€” v(k â€” n + l)] = (6.16) (6.17) (6.18) (6.19) {k)AT (6.20) where uT(/c) and A are given by (5.28) and (5.30) respectively. Therefore, it follows from (6.17) and (6.20) that yj(kT + s) â€” uT(k)A Tj Qjo r(s) (6.21) Finally, the required relation between the state modal function ]jj and the piece-wise constant input u follows from (6.21) as MkT + s) = uT{k)q.{s) where q.(s) = ATr;- Qjo r{s) (6.22) (6.23) Substituting definition 2.2 for y3- in (6.22), and moving all terms with the exception of Xj+1 to the right-hand side of the equation yields ~Xj+1{kT + s) = -J2Pjiy{{k - i)T + a) + g.T(*)Â£(*) (6.24) Observer Design 91 As seen from equation (6.24) it is necessary to obtain estimates of the vectors p. and q. in order to use the equation as an observer of the state variable Xj+i. Estimating the modal parameters from the system response to a piecewise-constant input has been considered in chapter 31. It has been shown that a time interval greater than or equal to nT should be allowed between any two successive input transitions. This restriction resulted from the fact that the continuous ARMA model used to estimate the modal parameters is valid only over the time interval of a constant input. An alternative method based on the estimated step response can remove this restriction on the frequency of input transitions. The continuous ARMA model in case of a step input is given by D*g{t) + itPjig{t - iT) = qjo ,j > 0. (6.25) t=i The methods of obtaining the output derivative modal parameters p}\,... ,pjn for j > 0 using Fourier functions or Walsh functions which are described in section 3.4 are applicable to the case of the ARMA model (6.25). The only difference is that the estimated step response samples will replace the measured system response samples. The estimates of q^.{s) can be obtained from (6.23). It is seen from (6.23) that under the assumption that there is an available Tj, i.e. there is an available estimate of p*, it is possible to obtain q^(s) from the scalar q]0 and the vector r(s) which is independent of the index j. Actually, it is sufficient to obtain an estimate of r(s) in order to obtain ^.(-s) from (6.23) because the scalar qj0 can be obtained from (6.15) as 9,-o = pTr(0). (6.26) The required estimate of r(s) can be obtained from (6.22) by setting j=0 as <loo r{s) (ATT0)-\(s) (6.27) A least squares estimate of the vector qo(s) required in (6.27) can be obtained from N > n + 1 equations like (6.22) formed from the input-output measurements for different k as Us) = {UTU)-Wy(s) (6.28) isee examples 3.1 and 3.2 2 see (6.19) Observer Design 92 where Vo(s) and U are given by (5.33) and (5.35) respectively. Comparing (6.28) with (5.32) gives the relation between the estimate of Q00{s) and the vector of step response samples g(s) as io{s) = ALg{s) (6.29) where L is the lower triangular matrix given by (5.34). Equation (6.29) gives the estimate of qo{s) as a linear transformation of the estimates of the step response samples. The estimate of q^s) is then obtained from the estimates of r(s) and gJ0 using (6.23) where q\j0 is obtained from (6.26). By estimating q^s) equation (6.24) can be used to obtain the state variable x;+1 in terms of the available input-output measurements . Equation (6.24) represent a single row in the vector equation x{kT + s) == Py{{k - 1)T + s) + Q{s)u{k) (6.30) where y((k â€” 1)T + s) is given by (4.20), the matrix P is given by (4.18) and the matrix Q{s) is given by 2Â» Q(s) (6.31) The dimension of P and Q(s) are nxn and nx(n + 1) respectively. The state observer (6.30) gives the system state x at time t = kT + s as a direct function of n output samples y((k â€” l)T + s),... ,y((k â€” n)T + s) and n+l input samples u(k),..., u(k â€” n + l). To obtain a state observer of the system state at the discrete sampling instants only, the time variable 5 is set to zero such that (6.30) becomes x{kT) = Py{{k - l)T) + Q{0)u{k). (6.32) The system state at the beginning of the kth sampling interval x(kT) is independent of the input u{k) applied to the system during that sampling interval, which means that the first column of Q(0) in (6.32) is zero. This can be verified by testing the first element of the rows of Q(0). According to (6.31), the (j + l)th row of Q(0) is j^.(0). The first element in q}(0) is given by (6.23) as if2,(o)=i?A'r,- Qjo 1(0) = qj0 - pjY(o) = 0. (6.33) Observer Design 93 The last equation follows from (6.26). Let the n dimensional vector u(k â€” 1) be formed from the elements of the n + 1 dimensional vector u(k) by deleting the first element, u(k), as follows: u{k-l) = {oI}u{k). (6.34) Similarly, let the nxn matrix Q be formed from the elements of the nx(n + 1) matrix 0(0) by deleting the first zero column such that (6.35) Thus the discrete-time observer follows from (6.32) as x[k) = Py{k - 1) + Qu{k - 1). (6.36) The time interval T has been ommited from the discrete-time observer since it is implicitly understood. It is seen from (6.36) that it takes only n sampling intervals to form the vectors y(k â€” 1) and u(k â€” 1) in order to be able to reconstruct the state vector x(k). This kind of observer is called a dead beat observer [Astrom 83], or a minimum time observer [O'Reilly 83]. It should be noted that although the observer (6.36) is discrete, it gives the state x of a continuous-time model of the system since x is defined in terms of the output and its derivatives which are meaningless in a discrete model case. To summarize, the design of the state observer in case of a piece-wise constant input takes the following steps: â€¢ Estimation of the output modal parameter vector pg. â€¢ Estimation of the vector of step response samples g(s). â€¢ Estimation of the state modal parameter vectors p , ... ,Pnl from the esti mated step response. â€¢ Obtaining ?0(s),. â€¢ â€¢, 5n_x("S) from the estimated step response and the modal parameter vectors using equations (6.29), (6.27), (6.26) and (6.23). â€¢ Forming the matrix P from rows PQ , â€¢ â€¢ â€¢ >P^_j-Q = 0(o) Observer Design 94 â€¢ If a continuous observer is sought, the matrix Q(s) is formed from rows qg(s),... ,qnl(s), while if a discrete observer is sought, the matrix Q is formed as in (6.35). Example 6.1 : Consider a third order system with transfer function rr(^ - 6.5s + 52 U [ ' ~ s* + 4s2 + 30s + 52 The system is identical to the system considered in example 3.1 with the exception that the system in example 3.1 has a zero at s = â€”4 while the system in this example has a zero at s = â€”2. Starting from rest, the system is excited with the piece-wise constant input +1, +2, -1, and -2. Each constant level lasts for a du ration 80rs where Ts = 0.01. It is desired to obtain a minimum time observer for this system to give the system state at the sampling instants x(kT), and midway between successive sampling instants x(kT + T/2) where T is the input sampling time taken equal to 20TS. Following the design steps listed above results in ob taining three matrices P, 0(0) and Q(T/2) listed in table 6.1. It is seen from the results that the first column of Q(0) is zero as expected. 6.3 Estimation of 6 The identification of the parameters of a system model has been considered in chapter 4. Estimates of the coefficients of the characteristic polynomial have been obtained in sections 4.3 and 4.4. By estimating these coefficients it is possible to obtain an estimate of the system dynamic matrix in phase variable form A. Hence, it remains to estimate the vector 6 in order to complete the identification of the state model of the system. The nonzero elements in h result in a difference between the state variables x;+i, j = 1,..., n â€” 1, and the corresponding output derivative P3y. The same difference results between the state modal functions and the corresponding output deriva tive modal functions. The output derivative modal functions and the state modal functions that result from a step input are constants gJ0 and qj0, j = 1,... ,n â€” 1 respectively. Estimates of these constants are obtained in section 6.2. The relation Observer Design 95 Matrix P -0.1555E+01 +0.1264E+01 -0.4494E+00 -0.3142E+00 +0.4343E+01 -0.2163E+01 +0.2621E+02 -0.2842E+02 +0.6032E+01 Matrix Q(0) +0.4985E-05 -0.9472E-06 +0.2734E-06 +0.1455E+00 +0.1475E+00 +0.1444E+01 +0.5780E+00 -0.2339E+00 -0.2850E+01 -0.3413E-01 -0.1643E+00 +0.4581E+00 Matrix Q{T/2) +0.3561E-01 +0.7309E+00 +0.1156E+01 +0.2257E+00 +0.9755E-02 +0.1185E+01 +0.1017E-03 -0.3505E+01 -0.4398E+00 -0.1211E-01 -0.5830E-01 +0.1626E+00 Table 6.1: Minimum time observer parameters between gJ0 and q\j0 follows from (2.33) by setting the input u to unity and all its derivatives to zero ; <ljo-qjo = bj J = l,...,n- 1. (6.37) Hence the estimates of bj, j = 1,... ,n â€” 1, can be obtained from (6.37) as the difference between the corresponding output derivative and state modal functions which are obtained from the estimate of the step response. The n th element in 6 cannot be obtained from (6.37) because, although the output derivative modal functions are defined for j > 0, the state modal functions are defined for j up to n â€” 1 only. The output derivative modal function of order n is given by definition as n qno = Dng{t) + J2Pm 9{t ~ iT). (6.38) Â«=i The relation between the output derivative of orders 1,... ,n and the system state variables is given by (2.11). Substituting in (2.11) for y by g, for u by unity and for all the derivatives of u by zero, it follows that the n th derivative of g is P ng{t) = a Tx{t)+ ~bj. (6.39) Observer Design 96 From definition 2.2 of the state modal function it followos that the state modal function that results from a step input is given by n Â£j+i(0 + YlPaSi1 - iT) = ~QiÂ° (6-40) t'=l Substituting for the state variables in (6.39) from (6.40) it follows that nâ€”1 n Dn9{t) = -E afao -Y,Pji9{* ~ &)) + k- (6.41) j=o Â«=i Substituting back for Vng in (6.38) and rearranging terms yields nâ€”1 n nâ€”1 Qno = ~ Â£ a,Qjo + Â£[Pni + Â£ o.jPji]g{t - iT) + bn. (6.42) j'=0 i=l j-0 However, from theorem 4.1, it follows that the term in rectangular brackets in (6.42) is zero, and therefore bn can be obtained as n-1 k = qno + Â£ ajqjo- (6.43) 3=0 The constant output derivative modal function qno can be obtained exactly the same way as other output derivatives modal functions qj0, j < n. The parameters a,j, j = 0,..., n â€” 1 are obtained from the modal parameters as in section 4.3. Example 6.2 : It is required to obtain the vector b for the system in example 6.1 using (6.37) and (6.43) . The estimates of p., qj0 and gJ0, ,7=0,..., nâ€”1, already obtained in example 6.1 will be used. The estimates of pn and qno are obtained from the estimated step response. Using these estimates in (6.37) and (6.43) gives the estimate of b listed in table 6.2. The true values of the elements of b are listed with the estimated values for comparison purposes. 6.4 Summary The main result of chapter 6 is the minimum-time state observer in continuous-time, (6.30), and in discrete-time, (6.36). Although (6.36) is discrete in nature, it gives the continuous-time system state x which is defined in terms of the derivatives Observer Design 97 ^ I b +0.7399E-02 0.0 +0.6453E+01 +0.6500E+01 +0.2624E+02 +0.2600E+02 Table 6.2: Estimated and true values of 6 of the continuous output. The minimum-time observer gives the system state in terms of the input-output samples in n sampling intervals. A method is given to determine the vector 6, of the state model, from the modal functions that correspond to the step response, and the estimates of the a-parameters obtained before in chapter 4. Chapter 7 Conclusion 7.1 Achievements A new method for continuous system identification is presented. The new method uses modal functions which are defined in terms of the system output, the system output derivatives and the system state. A necessary condition for the existence of the modal parameters, associated with the modal functions, is established in theorem 2.1. The properties of the modal functions are studied in section 2.4. It is shown that these functions are independent of the system poles and the system instantaneous state and that they are linear functions of the input modes. Methods of estimating the modal parameters from the available input-output measurements are presented in chapter 3. Orthogonal functions, Fourier functions and Walsh functions, are used to estimate modal parameter vectors. It is shown that only one stage of integration is required to determine one modal parameter vector at a time resulting in only one unknown added to the estimation problem in the case of Fourier functions1, and none in the case of Walsh functions2. This is an advantage over the known methods utilizing modulating functions which usually result in n additional unknowns [Saha 83], or trade the orthogonality of the modulating functions for a lower number of unknowns [Pearson 85]. It is also an advantage over the methods which utilize orthogonal functions and operational matrices for integration. These methods require the estimation of unknown initial conditions, and they suffer from accumulated errors due to multiple integrations 1see equation (3.30) 2see equation (3.37) 98 Conclusion 99 with truncated operational matrices. The output modal function is obtained from the system output by a simple filtering operation3. After showing that the output modal function is a particu lar response of the system which is composed of modes from the corresponding filtered input, the filtered input-output pair is used to identify parametric and nonparametric models of the system. A parametric model is identified in section 4.2 through a method that utilizes the output modal function. It is shown that the unavailable derivatives which used to be a problem in continuous system identification can be obtained from the output modal function analytically4. With a periodic waveform input, the output modal function, which in this case is periodic, is used in section 5.2 to obtain samples of the frequency response of the system. The transient component of the system response, which causes the systematic errors in the frequency response estimation methods using digital spectral analysis, [Wellstead 81], is annihilated by the use of the output modal filter. Estimates of samples of the step response of the system are obtained using the output modal function in section 5.3. Samples of the impulse response of the system are obtained in section 5.4 as the derivatives of the step response by a simple filtering operation that uses the output first derivative modal parameters5. Modal functions and modal parameters of the system output and the system output derivatives are used in chapter 4 to identify a parametric input-output model of the system. A nice property of using the modal parameters in parametric system identification is decoupling of the estimation equations with a smaller num ber of unknowns to be estimated from each equation. By estimating n+1 modal parameter vectors, each of dimension n, the vector a of the coefficients of the char acteristic polynomial of the system is obtained according to theorem 4.1 as the solution of a set of n algebraic equations formed from the estimates of the modal vectors. The vector b is obtained according to theorem 4.2 by solving another set of linear algebraic equations formed from the modal functions of the output and 3see Definition 2.1 4see equation (4.3) 5see equation (5.40) Conclusion 100 the output derivatives and the estimate of a. With the existing methods utilizing modulating functions or multiple integrations, there will be n parameters associ ated with the system poles, there can be n parameters associated with the system zeros plus n unknown initial conditions all contained in one equation. It can be difficult to realize well-conditioned estimation equations even for low n. Theorem 4.3 gives a proof that a data matrix formed from the free response of the system possesses the same eigenvalues as the system dynamic matrix A. From this data matrix the coefficients of the system characteristic polynomial are obtained recursively through a simple relation6 which is well known in matrix algebra but seems not to have been used before in system identification. The system poles are obtained by solving for the eigenvalues of a data matrix. Theorem 4.3 has been extended to include systems driven by continuous controllers. In section 6.3, a method is presented to estimate the vector b thus completing the identification of a continuous state model of the system in phase variable form. Based on the definition of the state modal functions, a procedure to design a minimum-time observer in case of a piece-wise constant input is given in section 6.2. The observer gives the continuous system state as a direct function of input-output samples in a number of sampling intervals equal to the order of the system. 7.2 Remarks on Future Research The choice of the modal time shifts affects the linear independence of the columns of xjj. This affects the accuracy of the estimates of the modal parameter vectors and consequently the accuracy of the estimates of the model unknown parameters derived from these vectors. Therefore, by understanding the exact impact of the modal time shifts on the estimation equations, they can be chosen to optimize the conditioning of the estimation problem. The effect of noise on the performance of the minimum-time state observer when used in a feedback control loop should be investigated. A Kalman filter could be used to obtain best estimates of the output samples that are used to obtain the system state. 6see equation (4.50) Conclusion 101 Since the modal time shifts T,, i=l,... ,n in the case of a piece-wise constant input is not restricted to Ti=iT, it is possible to obtain a parametric discrete model of the system with T^/c.T, k, > i. Such a model can be useful in designing predictive controllers yielding better responses. This thesis has dealt with single-input single-output systems. Extensions of the thesis results to multi-input multi-output systems does not seem to pose any theoretical difficulty. Modal functions can be determined for each output, as is the case for single-input single-output systems. The problem appears to be one of ordering output modes and accounting for free subsystem modes and coupling modes between subsystems. Bibliography [Astrom 83] K.J.Astrom and B. Wittenmark, 1983 . "Computer controlled sys tems: Theory and design". Prentice-Hall, NJ. [Bars 85] R. Bars and L. Keviczky, 1985 . "Self-tuning regulators based on nonparametric identification"'. Proc. IFAC Identification and sys tem parameter estimation, York, UK, 815-821. [Bohn 82] E.V. Bohn, 1982 . "Measurements of continuous-time linear system parameters via Walsh functions". IEEE Trans. Ind. Electron. IE-29, 88-46. [Chen, C. 75] C.F. Chen and C.H. Hsiao, 1975 . "Walsh series analysis in optimal control". Int.J. Control 21, 881-897. [Chen, W. 87] W.L-. Chen and C.-Y. Chung, 1987 . "New orthogonal operational matrix in block-pulse series analysis". Int.J.Systems Sci., 403-408. [Chung 87] H.-Y. Chung, 1987 . "Systen identification via Fourier series". Int.J.Systems Sci., 1191-1194-[Clarke 87] D.W. Clarke, C. Muhtadi and P.S. Tuffs, 1987 . "Generalized pre dictive control - Part I. The basic algorithm". Automatica 28, 187-148. [Crittenden 83] J.L. Crittenden, R.J. Mulholland and E.F. Martinez, 1983 . "Sam pling and data analysis of Prony's method applied to model iden tification". Int.J.Systems Sci. 14, 571-584-102 BIBLIOGRAPHY 103 [De Keyser 85] R. De Keyser, G. Van de Velde and F. Dumortier, 1985 . "A com parative study of self-adaptive long-range predictive control meth ods". Proc. IFAC Identification and system parameter estimation, York, UK, 1317-1822. [Douce 86] J.L. Douce, 1986 . "A note on frequency response measurements". Proc. IEE 188, 189-190. [El-Shafey 87] M.H. El-Shafey and E.V. Bohn, 1987 . "The use of modal functions for continuous system identification and state observer design". Int.J. Control 45, 1728-1736. [Hwang 82] C. Hwang and Y.P. Shih, 1982 . "Parameter identification via La-guerre polynomials". Int.J.Systems Sci. 18, 209-217. [Karanam 78] V.R. Karanam and P.A. Frick, 1978 . "Bilinear system identifi cation by Walsh functions". IEEE Trans. Aut. Control. AC-28, 709-718. [Levy 59] E.C. Levey, 1959 . "Complex curve fitting". IRE Trans. Aut. Con trol AC-4, 37-44. [Loeb 65] J.M. Loeb and G.M Cahen, 1965 . "More about process identifi cation". IEEE Trans. Aut. Control AC-IS, 859-361. [Marple 87] S. Lawrence Marple Jr., 1987 . "Digital spectral analysis: with applications". Prentice-Hall, NJ. [Noharski 85] Z. Noharski, L. Bogdan and J. Studzinski, 1985 . "Estimation of system structure and parameters from noisy sampled impulse re sponse". Proc. IFAC Identification and system parameter estima tion, York, UK, 1747-1758. [Oppenheim 75] A.V. Oppenheim and R.W. Schafer, 1975 . "Digital signal process ing". Prentice-Hall, NJ. [O'Reilly 83] J. O'Reilly, 1983 . "Observers for linear systems". Academic Press, London. BIBLIOGRAPHY 104 [Paraskev. 85] [Pearson 85] [Poggio 78] [Rangana. 87] [Rao 83] [Rao 751 [Rao 811 [Richalet 78] [Rouhani 82] [Saha 83] P.N. Paraskevopoulos, 1985 . "Legendre series approach to identi fication and analysis of linear systems". IEEE Trans. Aut. Control AC-SO, 585-589. A.E. Pearson and F.C. Lee, 1985 . "Efficient parameter identi fication for a class of bilinear differential systems". Proc. IFAC Identification and system parameter estimation, York, UK, 161-165. A.J. Poggio, M.L.van Blaricum, E.K. Miller and R. Mittra, 1978 . "Evaluation of a processing technique for transient data". IEEE Trans. Antennas and Propag. AP-26, 156-171. A.N. Ranganathan and N.S. Rajamani, 1987 . "Identification of lumped linear time-varying systems via Laguerre polynomials". Int.J'.Systems Sci. 18, 673-680. G.P. Rao, 1983 . "Piecewise continuous orthogonal functions and their application to systems and control ". Lecture Notes in Control and Information Sciences, Springer Verlag. G.P. Rao and L. Sivakumar, 1975 . "System identification via Walsh functions". Proc. IEE 122, 1160-1161. G.P. Rao and K.R. Palanisamy, 1981 . "Improved algorithm for pa rameter identification in continuous systems via Walsh functions". Proc. IEE 130, 9-16. J. Richalet, A. Rault, J.L. Testud and J. Papon, 1978 . "Model predictive heuristic control: Applications to industrial processes". Automatica 14, R. Rouhani and R.K. Mehra, 1982 . "Model Algorithmic Control (MAC); Basic theoretical properties". Automatica 18, 401-414-D.C. Saha and G.P. Rao, 1983 . "Identification of continuous dy namical systems: The Poisson moment functional approach". Lec ture Notes in Control and Information Siences, Springer Verlag. BIBLIOGRAPHY 105 [Sannuti 77] P. Sannuti, 1977 . "Analysis and synthesis of dynamical systems via block-pulse functions". Proc. IEE 124, 569-571. [Schaubert 79] D.H. Schaubert, 1979 . "Application of Prony's method to time-domain refiectometer data and equivalent circuit synthesis". IEEE Trans. Antennas and Propag. AP-27, 180-184-[Shanks 69] J. Shanks, 1969 . "Computation of the Fast Walsh-Fourier Trans form". IEEE Trans. Computers C-18, 457-459. [Shinbort 57] M. Shinbort, 1957 . "On the analysis of linear and nonlinear sys tems". Trans. ASME 79, 547-552. [Sinha.N.K. 74] N.K. Sinha, A. Sen and R.H.Y. To, 1974 . "Time domain approx imation using digital methods". Int.J.Systems Sci. 5, 278-382. [Sinha.N.K. 72] N.K. Sinha, 1972 . "Estimation of transfer function of continuous system from sampled data". Proc IEE 119, 612-614-[Sinha.P.K. 84] P.K. Sinha, 1984 . "Multivariable control: An introduction". Mar cel Dekker, NY. [Soderstrom 81] T. Soderstrom and P.Stoica, 1981 . "Comparison of some Instru mental Variable methods-Consistency and accuracy aspects". Au-tomatica 17, 101-115. [Strang 80] G. Strang, 1980 . "Linear algebra and its applications". Academic Press, London. [Tadokoro 78] Y. Tadokoro and T. Higuchi, 1978 . "Discrete Forier Transform computation via the Walsh Transform". IEEE Trans. Acoust. Speach Signal Process. ASSP-26, 236-240. [Trivett 81] D.H. Trivett and A.Z. Robinson, 1981 . "Modified Prony method approach to echo-reduction measurements". J. Acoust. Soc. Am. 70, 1166-1175. BIBLIOGRAPHY 106 [Tzafestas 77] S.G. Tzafestas and N. Chrysochoides, 1977 . "Time-varying reac tivity reconstruction via Walsh functions". IEEE Trans Aut. Con trol AC-22, 886-888. [Wang 87a] M.-L. Wang, R.-Y. Chang and S.-Y. Yang, 1987 . "New approach for parameter identification via generalized orthogonal polynomi als". Int.J. Control 46, 569-579. [Wang 87b] M.-L. Wang, R.-Y. Chang and S.-Y. Yang, 1987 . "Analysis and parameter estimation of bilinear systems via generalized orthogo nal polynomials". Int.J. Control 46, 719-729. [Wellstead 81] P.E. Wellstead, 1981 . "Non-parametric methods for system iden tification". Automatica 17, 55-69. [Whitfield 86] A.H. Whitfield, 1986 . "Transfer function synthesis using frequency response data". Int.J. Control 43, 1413-1426. [Young 70] P.C. Young, 1970 . "An Instrumental Variable method for real time identification of a noisy process". Automatica 6, 271-287. [Young 81] P.C. Young, 1981 . "Parameter estimation for continuous-time models-A survey". Automatica 17, 23-39. [Young 79] P.C. Young and A.J. Jakeman, 1979 . "Refined Instrumental vari able methods of recursive time-series analysis. Part I: single-input single-output systems". Int.J. Control 29, 1-30.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Linear continuous-time system identification and state...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Linear continuous-time system identification and state observer design by modal analysis El-Shafey, Mohamed Hassan 1987-12-31
pdf
Page Metadata
Item Metadata
Title | Linear continuous-time system identification and state observer design by modal analysis |
Creator |
El-Shafey, Mohamed Hassan |
Publisher | University of British Columbia |
Date | 1987 |
Date Issued | 2010-09-24T03:35:30Z |
Description | A new approach to the identification problem of linear continuous-time time-invariant systems from input-output measurements is presented. Both parametric and nonparametric system models are considered. The new approach is based on the use of continuous-time functions, the modal functions, defined in terms of the system output, the output derivatives and the state variables under the assumption that the order n of the observable system is known a priori. The modal functions are obtained by linear filtering operations of the system output, the output derivatives and the state variables so that the modal functions are independent of the system instantaneous state. In this case, the modal functions are linear functions of the input exponential modes, and they contain none of the system exponential modes unlike the system general response which contains modes from both the system and the input. The filters parameters, the modal parameters are estimated using linear regression techniques. The modal functions and the modal parameters of the output and its derivatives are used to identify parametric input-output and state models of the system. The coefficients of the system characteristic polynomial are obtained by solving n algebraic equations formed from the estimates of the modal parameters. Estimates of the parameters associated with the system zeros are obtained by solving another set of linear algebraic equation. The system frequency response and step response are estimated using the output modal function. The impulse response is obtained by filtering the estimated step response using the output first derivative modal parameters. A new method is presented to obtain the system poles as the eigenvalues of a data matrix formed from the system free response. The coefficients of the system characteristic polynomial are obtained from the data matrix through a simple recursive equation. This method has some important advantages over the well known Prony's method. The state modal functions are used to obtain a minimum-time observer that gives the continuous-time system state as a direct function of input-output samples in n sampling intervals. |
Subject |
Linear time invariant systems Discrete-time systems |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2010-09-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0097969 |
URI | http://hdl.handle.net/2429/28666 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- UBC_1988_A1 E47.pdf [ 6.07MB ]
- Metadata
- JSON: 1.0097969.json
- JSON-LD: 1.0097969+ld.json
- RDF/XML (Pretty): 1.0097969.xml
- RDF/JSON: 1.0097969+rdf.json
- Turtle: 1.0097969+rdf-turtle.txt
- N-Triples: 1.0097969+rdf-ntriples.txt
- Original Record: 1.0097969 +original-record.json
- Full Text
- 1.0097969.txt
- Citation
- 1.0097969.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 19 | 1 |
China | 18 | 12 |
Russia | 8 | 0 |
Republic of Lithuania | 6 | 0 |
Japan | 3 | 0 |
Iran | 2 | 0 |
Luxembourg | 2 | 0 |
Macedonia | 2 | 0 |
France | 1 | 0 |
Egypt | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 18 | 2 |
Ashburn | 10 | 0 |
Shenzhen | 10 | 11 |
Beijing | 7 | 0 |
San Francisco | 3 | 0 |
Tokyo | 3 | 0 |
Saint Petersburg | 2 | 0 |
Strassen | 2 | 0 |
Shtip | 2 | 0 |
Washington | 2 | 0 |
Livingston | 1 | 0 |
Chengdu | 1 | 1 |
Wilmington | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0097969/manifest