UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Linear systems identification and optimization with application to adaptive control Hanafy, Adel Abdel Raouf 1972-03-21

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
831-UBC_1972_A1 H35.pdf [ 5.83MB ]
Metadata
JSON: 831-1.0101279.json
JSON-LD: 831-1.0101279-ld.json
RDF/XML (Pretty): 831-1.0101279-rdf.xml
RDF/JSON: 831-1.0101279-rdf.json
Turtle: 831-1.0101279-turtle.txt
N-Triples: 831-1.0101279-rdf-ntriples.txt
Original Record: 831-1.0101279-source.json
Full Text
831-1.0101279-fulltext.txt
Citation
831-1.0101279.ris

Full Text

LINEAR SYSTEMS IDENTIFICATION AND OPTIMIZATION WITH APPLICATION TO ADAPTIVE CONTROL ADEL ABDEL-RAOUF HANAFY B.Sc, University of Cairo, Cairo, Egypt, 1965 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of Electrical Engineering We accept this thesis as conforming to the required standards THE UNIVERSITY OF BRITISH COLUMBIA July, 1972 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 shal1 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 representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Depa rtment The University of British Columbia Vancouver 8, Canada Date 26-1-7 Z ABSTRACT This thesis is concerned with the problem of identifying and controlling linear continuous systems. Algorithms which are feasible for real-time digital computations are developed for the design of both the controller and the identifier. The generalized equation error is shown to be applicable to a mean-square method of rapid digital estimation of linear system para meters. Due to the imposed structure of the estimator the manipulation of high order matrices is avoided. Examples illustrate the effective ness of the estimator for a variety of cases dealing with quantization noise as well as measurements noise. In some cases, this digital identi fier requires the computation of the generalized inverse of a matrix. A simple algorithm for computing the generalized inverse of a matrix is developed. This algorithm eliminates the need for Gram-Schmidt ortho-gonalization and its associated (in the interest of accuracy) reortho-gonalization as new vectors are introduced. A two-stage estimator is developed for estimating time-invariant and time-varying parameters in linear systems. During the second stage, the time-varying parameters are considered as unknown control inputs to a linear subsystem of known dynamics. For the first stage, the digital identifier is shown to be effective in identifying the time-invariant parameters. Numerous examples illustrate the effectiveness of the method. To design a feedback controller a method of successive approxi mations for solving the two point boundary value problem for optimum constant gain matrices is developed. The method is shown to be compta-tionally equivalent to a deflected gradient method. Convergence can ii always be achieved by choice of a scalar step-size parameter. An on line approximate method is developed which appears suitable for systems whose parameters must be identified. The two point boundary value problem is replaced by an algebraic problem, the solution of which gives a sub-optimal constant gains. The problem of trajectory sensitivity reduction by augmented state feedback is shown to be well-posed if constrained structure of the gain matrices is taken. The simple structure of constant gain matrices is considered. Based on the developed identifier and controller, one possible strategy for optimal adaptive control which is particularly attractive from an engineering point of view is studied. The strategy is to identify system parameters at the end of an observation interval and then to use the parameters to derive an "optimal" control for a sub sequent control interval. Two examples are considered in order to illustrate the effectiveness of this strategy. iii TABLE OF CONTENTS Page ABSTRACT ii TABLE OF CONTENTS iv LIST OF TABLES vLIST OF FIGURES vii ACKNOWLEDGEMENT ix 1. INTRODUCTION 1 2. RAPID DIGITAL IDENTIFICATION OF LINEAR SYSTEMS 2.1 Introduction 8 2.2 Problem Formulation 9 2.3 Recursive and Block Data Processing 14 2.4 Measurement Requirements 17 2.5 Applications 18 2.6 Comparison with a Limited Memory Filter 27 2.7 Some Statistical Properties of the Digital Identifier... 30 2.7.1 Correlation of Noise Samples 32.7.2 The Error Covariance Matrix 6 2.7.3 Estimation of the Confidence Interval Length 37 2.8 Discussion 42 3. RECURSIVE COMPUTATION OF THE GENERALIZED INVERSE OF A MATRIX. 3.1 Introduction 44 3.2 Properties of the Pseudo-Inverse 43.3 Recursive Evaluation of R^ 46 3.4 Recursive Evaluation of A1 9 3.5 Examples 50 4. TWO-STAGE ESTIMATION OF TIME-INVARIANT AND TIME-VARYING PARAMETERS IN LINEAR SYSTEMS 4.1 Introduction.. 56 4.2 Problem Formulation 7 4.3 Estimation of Av. Algorithm 1 59 4.4 Estimation of a Periodic \ Algorithm II 62 4.5 Estimation of a Periodic Av. Algorithm III 5 4.6 Two-Stage Estimation of Af and Av 67 4.7 Examples 64.7.1 Example 1 70 4.7.2 Example 2 1 Page 4.7.3 Example 3 71 4.7.4 Example 4 2 4.7.5 Example 5 3 4.7.6 Example 65. OPTIMUM CONSTANT FEEDBACK GAINS FOR LINEAR SYSTEMS AND THE PROBLEM OF TRAJECTORY SENSITIVITY REDUCTION 5.1 Introduction 76 5.2 Problem Formulation 8 5.3 Unified Treatment of Different Formulations 82 5.4 Gradient Descent and Successive Approximations 4 5.5 Optimum Constant Gain: Trajectory Sensitivity Reduction. 88 5.6 On-Line Computation of Gains 90 5.7 Examples 93 6. OPTIMAL ADAPTIVE CONTROL 6.1 Introduction 102 6.2 Optimum Control During the First and Subsequent Optimization Intervals 106.3 Examples 104 6.3.2 Example 1 105 6.3.3 Example 2 7 6.3.4 Discussion . 110 7. CONCLUSIONS 116 APPENDIX I 9 APPENDIX II .' . 122 APPENDIX III.- 125 APPENDIX IV 7 APPENDIX V 130 REFERENCES • 132 v LIST OF TABLES Pagi Table (4.1) 75 Table (4.2)Table (4.3)Table (4.4) 75 LIST OF FIGURES Figure Page 1.1 The Proposed General Optimal Adaptive System 2 2.1 Example 1. Estimation of a^ with Quantization Noise 21 2.2 Example 1. Estimation of with Quantization Noise 21 2.3 Example 1. Estimation of a^ with Quantization and Unbiased Measurement Noise 22 2.4 Example 1. Estimation of a^ with Quantization and Unbiased Measurement Noise2.5 Example 1. Estimation of a-^ with Quantization and Biased Measurement Noise 23 2.6 Example 1. Estimation of a^ with Quantization and Biased Measurement Noise2.7 Example 2. Estimation of a^ with Quantization noise 24 2.8 Example 2. Estimation of a^ with Quantization Noise 24 2.9 Example 2... Estimation of. b^ with. Quantization. Noise.. 25. 2.10 Example 2. Estimation of a^ with Quantization and Unbiased Measurement Noise 25 2.11 Example 2. Estimation of a^ with Quantization and Unbiased Measurement Noise 6 2.12 Example 2. Estimation of b^ with Quantization and Unbiased Measurement Noise 22.13 Example 2. Estimation of a., with Quantization and Unbiased Measurement Noise Using a Standard Limited Memory Filter... 28 2.14 Example 2. Estimation of a~ with Quantization and'Unbiased' Measurement Noise Using a Standard Limited Memory Filter... 28 2.15 Example 2. Estimation of b1 with Quantization and Unbiased Measurement Noise Using a Standard Limited Memory Filter... 29 4.1 The Recursive Observation-Estimation Intervals 64 5.1 Example 1. Control as a Function of Time For: (a) Lee's Method. (b) Constant Gain and Augmented States. (c) Constant Gain and Non-Augmented States 98 5.2 Example 1. Trajectory Senstivity Function z^ For: (a) Lee's Method. (b) Constant Gain and Augmented States, (c) Constant Gain and Non-Augmented States 99 Figure Page 5.3 Example 1. Trajectory Sensitivity Function z.^ for: (a) Lee's Method. (b) Constant Gain and Augmented States. (c) Constant Gain and Non-Augmented States 100 5.4 Example 2. Block Diagram of a Sychronous Generator Voltage Control System 101 6.1-a The Initialization Interval and the Operation Period 103 6.1-b The Successive Observation Intervals 106.1-c The Successive Optimization Intervals (X, = r) 103 6.1-d The Successive Optimization Intervals (I = 2r) 103 6.2 Example 1. Identification of Parameter a^ with Quantization., Noise 112 6.3 Example 1. Identification of Parameter a^ with Quantization Noise6.4 Example 1. The Optimum Piecewise Constant Feedback Gain K^... 113 6.5 Example 1. The Optimum Piecewise Constant Feedback Gain K^... 113 6.6 Example 2. Identification of Parameter a with Quantization Noise 114 6.7 Example 2. Identification of Parameter a^ with Quantization Noise6.8 Example 2. Identification of Parameter a^ with Quantization Noise 114 6.9 Example 2. The Optimum Piecewise Constant Feedback Gain K^... 115 6.10 Example 2. The Optimum Piecewise Constant Feedback Gain l^... 6.11 Example 2. The- Optimum Piecewise Constant Feedback Gain K,.-..- 115 viii ACKNOWLEDGEMENT I am indebted to my supervisor, Dr. E. V. Bohn, for his continuous guidance and assistance throughout the period of this research. Thanks are due to Dr. V. G. Haussmann for his valuable comments on the draft. Reading the draft by Dr. M. S. Davies and Dr. G. F. Schrack is duely appreciated. 1 also wish to thank Miss Linda Morris for her patience in typing the manuscript. The financial support by the National Research Council of Canada, grant 67-3134 (1968-1971), the Scholarship by the B. C. Hydro Company (1970), and the U. B. C. Graduate Scholarship (1971-1972) are appreciated. ix 1 1. INTRODUCTION An optimal adaptive control is a control that will cause a plant to operate in an "optimum" manner without a Priori knowledge about the constraints on either the magnitude or the rate of variation of the plant uncertainties. These uncertainties may be in the plant parameters, initial conditions, and/or noisy measurements. One possible mathematical design approach for an adaptive control is to use statistical methods. Plant uncertainties are then described by probability density functions. Such statistical methods are generally not computationally feasible even when a priori knowledge of the statistics of the plant uncertainties exist. In most practical situations, simplifying assumptions must be introduced and these detract considerably from the usefulness of a mathematically precise statistical approach. An alternative approach is the optimally sensitive control [1] where no a priori knowledge of statistical parameters is required. However, the range of plant parameter and/or initial condition variations is restricted. This thesis is concerned with the development of an optimally sensitive feedback control which does not require detailed statistical data for its realization. However, in order to increase the range of optimal control action, parameter estimation is used. Algorithms which are feasible for real-time digital computations are developed for the design of both the controller and the identifier. Figure (1.1) illustrates the general structure of the optimal adaptive control system which has been developed in this thesis. The -n LEAST Sq. I.C ESTIMATOR A-NOISE 11 FEEDBACK GAINS Y- -1 SENSITIVITY L_ ^ M?Z?£Z. ^ PLANT DIGITAL IDENTIFIER ADJUSTABLE MODEL ~—7> 11 C=3 /7/.7F/P r 11 'L -SCOMPOSITE T ^JDENTIFI^ _p OUTPUT FIG. (I.I) THE PROPOSED GENERAL OPTIMAL ADAPTIVE SYSTEM, 3 blocks drawn in dashed lines represent elements which may not be required in many practical applications. This optimal adaptive control system can generally be divided into two main parts; the estimator and the controller. A detailed discussion of the components which comprise these two parts together with the thesis layout is considered In the following two sections. 1.1 The Estimator The estimator is used to obtain a "best" estimate for the unknown plant parameters and initial conditions. It is also used to obtain a "best" estimate for the noisy and/or the unmeasured plant state variables. These estimates are needed in order to construct a reasonably accurate deterministic linear model of the plant. A deterministic linear model is necessary if use is to be made of the powerful tools that have been developed from optimal control theory. Linear continuous models are considered in this thesis. The choice of linear models is made because of their relative simplicity from the instrumentation and computation point of view. This simplicity is largely due to the high degree of development of the theory of linear systems. Even when the plant is nonlinear, however, the use of quasilinearization techniques can result in a linearized model of the plant. Over a fixed period of time-, such, a linear... model is. of ten. adequate. to. des cribe the dynamic behaviour of a nonlinear plant. The estimator is composed of four major elements (see Figure (1.1)). These four elements are: a- Initial Conditions Estimator This estimator is a standard linear least-squares estimator which minimizes 4 a quadratic error function based on a system of linear algebraic equations. It is necessary to obtain a good estimate of the plant initial conditions since the values of the controller constant feedback gains depend on initial conditions. This estimator is used whenever the plant states are conta minated by noise and/or some of these states are unmeasurable. A linear least-squares estimator has been chosen due to its simplicity. Also, no apriori knowledge of noise statistics is needed. b - Filter The filter gives a continuous estimation of the plant states. If some of the states are unmeasurable then a Luenberger observer is a good choice to reconstruct the unavailable states. On the other hand a Kalman filter can be used for best estimation of the noisy states. c - Digital Identifier A continuous-time parameter identification method that has been extensively investigated is the parameter-tracking technique [2,3,4], This is based on a steepest descent parameter adjustment. These parameter-tracking identifiers have been largely developed for use with analog computers. Practical on-line identification, however, requires the use of digital computers and the discrete version of parameter-tracking introduces several difficulties. The optimum gain or optimum step-size depends on the input signals and is difficult to determine on-line [5]. In Chapter 2 a digital identifier is developed which gives rapid and accurate estimation of the unknown, time-invariant parameters of the plant. The identifier uses data obtained from augmented input controls and output states which results in parallel data processing. 5 A discretized version of the continuous augmented model is used. Digital simulation results have shown that the identifier developed in this thesis has very good rates of convergence. This digital identifier requires the computation of a matrix inverse which may not exist in some cases. In such cases a generalized-inverse of a matrix is required. A number 7 of methods have been proposed for finding the generalized inverse of a matrix [6, 7, 8], These methods generally consist of some combination of a Gram-Schmidt orthogonalization process, the expansion of products of matrices, and the solution of sets of linear algebraic equations. In Chapter 3 an algorithm is developed which replaces these operations and computes- the generalized inverse from simple recursive formulas. Inherent in the method is a simple check on the linear inde pendency of the column vectors of the matrix. d - Composite Identifier A two-stage estimator is developed in Chapter 4. In the first stage the digital identifier is used to obtain an estimate for the time-invariant parameters. In the second stage the unknown time-varying parameters are considered as control inputs and a linear state regulator quadratic cost function dynamic optimization problem is formulated. The solution of the associated two-point boundary-value problem for the optimum control" results in an estimate for the time-varying parameters. This two-stage identification algorithm is suitable for the cases where only a few of the unknown parameters of the plant are time-varying. The use of the complicated techniques developed for the identification of time-varying parameters is not justified for such cases. The usual approach taken to identify time-varying parameters is to consider them as unknown augmented states described by differential equations of known form and then use some form of state estimation [9]. Alternatively, integral transforms have been used [10] to extend the capability of parameter-tracking system to track time-varying parameters. The implementation difficulties for such techniques make their use questionable if not completely unpractical. 1.2 The Controller In Chapter 5 the design of a controller with optimum constant feedback gains is discussed. A method of successive approximations is developed and it is shown to be equivalent to a deflected gradient des cent approach with a variable weighting matrix. This property is utilized to develop an algorithm for off-line design studies which results in rapid convergence. An approximate method is developed to replace the nonlinear two-point boundary-value problem by an algebraic optimization problem. This method appears to be suitable for on-line computation of suboptimum constant feedback gains (see e in Figure (1.1)). Trajectory sensitivity function- minimization is an important requirement... The in troduction of the sensitivity functions in the feedback structure of the adaptive controller makes it possible to up-date the constant feedback gains less frequently and makes the plant controller less sensitive to modelling errors. The sensitivity functions are generated by the sensi tivity model shown by f in Figure (1.1). The parameters of the sensitivity model are up-dated each time a new estimate of the unknown plant parameters is made. In Chapter 6 a recursive observation-optimization strategy is described. During a finite observation Interval, measurements of the state are taken and used to determine estimates of the plant parameters. The linear model parameters are adjusted according to these new estimates. The linear model is used to compute a control strategy which is applied in a subsequent control interval. This scheme seems particularly attractive for adaptive systems. In summary, the objective of this thesis is to investigate practical methods for designing an optimal adaptive control system based on a reasonably large domain of plant uncertainties. Emphasis is placed on techniques which can be implemented on real-time digital process control computers. To reduce computational requirements, the control strategy should not require frequent up-dating. This objective has been accomplished by developing a rapid digital identifier for continuous-time systems and by developing simple algorithms for the computation of optimally sensitive feedback gains. 8 2. RAPID DIGITAL IDENTIFICATION OF LINEAR SYSTEMS 2.1 Introduction Digital control of adaptive systems is feasible only if parameter identification and evaluation of control laws can be per formed in real time. To reduce the complexity of the computations, simplifying assumptions must be made. One feasible method is to identify the system as an approximate linear system within a given time interval and then use a linear optimal control [5,11], For such a control policy to be effective in an adaptive system the identification should be rapid. A continuous-time parameter identification method that has been extensively investigated is the parameter tracking technique [2,3,4]. This is based on a steepest descent parameter adjustment of a quadratic error function or a Liapunov function. The error can be taken to be a response or a generalized equation error. The generalized error approach is attractive since it results in an exact continuous-time steepest descent adjustment. Furthermore, by use of state variable filters, it allows considerable flexibility in the generation of additional input-output signals and additional- equation errors. These signals play an important:- role- in the estimator that will be discussed in this Chapter. They are considered to be the time response of augmented input controls and output states. These states can be chosen so that the quadratic error function has a unique minimum in parameter space which then allows very rapid identi fication. Practical on-line identification, however, requires the use of digital computers and the discrete version of parameter tracking introduces several difficulties. 9 The optimum gain or step-size depends on the input signals and is difficult to determine on-line [5]. Furthermore, a steepest descent search for a minimum does not generally result in rapid identi fication. A popular method for on-line identification by digital computers is the mean-square technique. In the linear system case all state variables must be measured and the estimated parameters are biased. Furthermore, if rapid identification is attempted, difficulties can arise with ill-conditioned matrices. A digital estimation method is developed here which can overcome these difficulties by use of augmented input controls and output states. Due to the imposed structure of the estimator the manipulation of high-order matrices, which is generally the case in the mean-square approach, is avoided. 2.2 Problem Formulation Consider a completely observable and controllable continuous-time linear system x = Ax + Bu + w (2.1) where x is an n-state vector, u is an m-control vector and w is additive plant noise. Let T be the sampling interval and let tQ, t^, t^, be a set of sampling instants. For the moment, to ease the notational burden, let z^(k) = z^(t^) be the sampled measurement at t^ where z^t) = x(t) + v(t) (2.2) and where v(t) is additive measurement noise. (The estimator to be developed does not require that all states be measured as implied by (2.2). The general measurement requirements are discussed in Section (2.4) Let Zi+l(t) - fo\^d'> U£+l(t) = i\™*> Wfc) = l\^d^ <2'3> The vectors z^, u^, (£ = 2), are defined as augmented states and aug mented controls, respectively. It is seen from (2.1), (2.2) and (2.3) that z„ = Az„ + Bu„ + w. % I I I (2.4) A A • where u^ = u, w^ = v - Av + w. Integrating (2.4) over the interval tk = x = tk+1 and letting Z£(k) = z£(tk>, U£(k) - u£(tk), yields Z£(k+1) = AZ^(k) + BU£(k) + N£(k) (2.5) where A = I + TA, B = TB and Vk) " [ k+1t5(z£<T) " z£<tk)) + B(uJl(T) " u£(tk}) + VT)]dT (2'6) k The discrete set of equations forms the basis for all further discussion. To ease the notational burden all n(n+m) elements of A and B are consi dered unknown. The problem is to estimate A and B from a set of measured Z^(j)» U-j^(j) and computed (from (2.3)) augmented states and controls. Let Vk)- -z4(k) A c = b' n (2...Z) where prime denotes transpose and where a^ and b^ are the j-th n and m row vectors of A and B respectively. These rows are ordered into an 2 (n + nm) vector c as shown. Equation (2.5) can then be written as z£(k+i) = e£(k)c + N4(k) (2.8) where 0 Yj(k) 0 •Yi(k) (Yl(k))Bn (2.9) (Due to the frequent occurrence of matrices of block-diagonal form it is convenient to introduce a special notation. The first subscript B is used to indicate that a block diagonal matrix is to be formed from the vector or matrix within the bracket. The second subscript n indicates 2 the number of blocks. Note that 9 (k) is an nx(n + nm) matrix). JO Consider the sampling instants k-N+1, k, so that N blocks of data are available and let Z £ ^N Z(k) Z (k-N+1) J Tc + N(k-l) N(k-N) (2.10) where N e(k-i) K • , 6(k) = , Z(k) e(k-N) n+m mm " zx(k) Zn+m<k> (2.11) The matrix 9(k) is n(n+m) x n(n+m) and the vector Z(k) is n(n+m) x 1. The matrix 6^ is of dimension nN(n+m) x n(n+m) as can be seen from (2.11). It is seen from (2.7), (2.9) and (2.11) that the measured and computed augmented states and controls are used to form the matrix 8(k) and the vector Z(k) at the k-th sampling instant. These quantitites are evaluated for N sampling instants and are used to form the matrix 9^ and the vector Z^. The problem is to determine an estimate of the parameter vector c given 9^ and Z^, where these quantities are related by (2.10). The simplest approach is to define the optimum estimate to be the estimate which minimizes the quadratic form (Z^ - 8^0)'(Z^ - Q^c). The solution is the minimum mean-square (also called least-squares) estimate [6] c = PM6*Z.T N N N where h1 = < w (2.12) (2.13) Since the basic structure of all linear mean-square estimators is the same it is of interest to compare the estimator (2.12) with pre viously proposed estimators [5, 12-16]. Mayne has used the Kalman filter approach. This requires that (2.10) be interpreted as a measurement equation with 6^ as a deterministic measurement matrix. This is the case only if there is no measurement noise and if all states are measured. A further and very significant difference is in the definition of 0^ and N(k). These investigators take 0 00 91(k-l) 01(k-n-m) , N(k) Nx(k-1) N (k-n-m) (2.14) (References [2, 13-15], furthermore, treat only the scalar measurement case. Their results do not directly apply to the above case). Since (see (2.5)) E[Z1(j+l)N^(j)] = E[N1(j)NjL(j)] •* 0 (2.15) it follows that the estimate is biased. This has been pointed out pre viously [13]. In order to eliminate the bias, sufficient time must elapse between samples. This can make rapid identification difficult. Even when the bias is small enough to give acceptable identification, further difficulties can arise with data in the form [2.14]. If a small sampling interval is chosen to achieve rapid identification there may not be a significant difference between adjacent rows of 6(k) resulting in an ill-conditioned matrix. These difficulties tend to be avoided if data in the form (2.10) is used. Some insight into the estimation bias associated with (2.12) can be obtained by considering the noise samples N^(j), N^(k), (p = 1, ..., n + m; j ^ k) , as independent. It follows from this assumption and (2.5) that E[Z0(k)N'(k)] = 0, (Z, p = 1, .... n + m) (2.16) Consequently 6.. at sampling time t, is independent of N (k) and the estimate _L ic p (2.12) is unbiased when N = 1. From (2.6) it is seen that quantization noise is an important component of N^(k) making an analytic investigation into the dependency of the noise samples intractable. Some general observations can, however, be made. Successive integrations performed on w (T) not only reduce the effective amplitude of the noise but also reduces the correlation between noise samples N^(k), (£ ~1, ..., n + m). The truth of this statement for a certain class of noise samples is shown in Section (2.7). Consequently, for a given k, the quantization noise terms of (2.6) for different I tend to be uncorrelated. Con sequently, little if any difficulty is to be expected with estimation bias. All examples tried verified this expectation. The essential difference between the two estimators can be seen from (2.15) and (2.16). If £ = 1, sequential data processing must be used. If T is chosen large to overcome ill-conditioning then quantization noise can become excessive. The combination of heavy noise and estimation bias (2.15) can result in failure of the conventional estimator to give acceptable estimation. These difficulties tend to be avoided if augmented states and augmented controls are chosen (£ = n + m) so that parallel data processing can be used. Heuristic engineering justifications can also be given for the use of augmented states and controls in linear system parameter identifi cation. If combined parameter estimation and optimal control is attempted signal-to-noise problems can arise [12], Optimal control attempts to keep state deviations from a zero reference state small which decreases the signal-to-noise ratio. Integrators are low-pass filters which improve signal-to-noise ratio. The data used is then less noisy. Such a situa tion has been reported by Bakke and he suggests using a hybrid identifi cation method [16]. Basically his hybrid method is to use data in the form (2.14) and then replace Y^, 6^ by X^y 82* Two different estimates are taken and a decision is then, made on the best estimate. However, a full set of augmented states is not considered, nor is the data combined in a mean-square sense. 2.3 Recursive and Block Data Processing The estimator (2.12) operates in a block data processing mode. A new estimate is made after every N additional data blocks. Other processing modes can be used such as the recursive method suggested by Kishi where old data blocks are deleted and new data blocks introduced. Let, 9N-R ZN-R -N > -XT (2.17) - -R ZN-R . z = ' N 9 z R R It follows from (2.13) that PN' = PN-R + 9R9R (2.18) Suppose that all quantities in (2.12) are known and it is desired to determine ^ from these quantities (This amount to determining the best estimate when R data blocks are deleted). Making use of (2.18) it follows from (2.12) that 9N-RZN-R - ^ " W -N-R (2"19) By partitioning (2.12) in the manner indicated by (2.17), it is seen that °N = WR + PN6N-RZN-R (2'20) substituting (2.19) into (2.20) yields CVR" (I " VRV^N " VRV (2'21> With the aid of the identities <X " PN6R6R)"lpNeR = PN6R (I " Wi^ <X ~ VRV"1 = 1 + (I ~ PNeR0R)_lpNeReR (2'22) (2.21) can be reduced to the standard from for recursive estimation CVR = *N + PN9R (I " eRPN9R)_1 (VN " ZR> (2'23) applying the matrix inversion lemma to (2.18) yields PN-R = PN + Vi (I " WR^VN (2'24) Equation (2.23) gives the best estimate if the last R data blocks are deleted and (2.24) gives the associated P-matrix. If R new data blocks are sampled, a new Z and 6D are available and new values of c^ and P^ can be computed. The equations giving these values are obtainable from (2.23) and (2.24) by interchanging N with N - R, replacing 6_ by -8_ and leaving 0' unchanged. With these changes (2.23) K K is. and (2.24) gives the Kalman filter algorithm (since the measurement matrix is a random variable the Kalman filter is not optimal). The estimators proposed in References [5, 12-15] are based on a scalar measurement. The bracketed quantity in (2.24) is then a scalar and no matrix inversion is required. However, many scalar samples are required to estimate the parameters within a given confidence interval. This can result in slow identification. The matrix to be inverted in (2.24) is of dimension Rn(n + m) x Rn(n + m). This is a consequence of the use of a vector measurement of augmented states. In effect, rapid identification is attempted by the parallel processing of augmented state samples. The price paid for parallel processing is in the high order matrix to be inverted in the Kalman filter algorithm. An alternative recursive approach is therefore desirable. It follows from (2.9), and (2.14) that n+m k-1 SL=1 J=k-N where n+m k-1 QN = £ £ Vj)-Yi (j) (2'26) 1=1 j=k-N consequently PN = %\n (2'27) Equation (2.27) gives the interesting result that only an (n + m) x (n + m) matrix need be inverted to find V^. (In general, to determine 2 2 2 n + nm parameters, an (n + nm) x (n + nm) matrix has to be inverted). The structure of (2.27) allows a recursive algorithm to be developed which gives i-11 terms of QN"''. Consider (2.26) and define Vn = Y„(k - N), S, = Q_T - V.V' S. = S. - V. V'., I I ' 1 N 11 j J-1 J j' S = Q 1 (2.28) n+m N-1 applying the matrix inversion lemma yields s-1 = (Ql) - v^p-1 - QJ1 + Q"1 v1ci1 + v^r^o"1 (2.29) Note that the quantities in brackets on the right hand side of (2.29) are scalars. The set of equations (2.29) define a recursive algorithm for finding given Repeating the algorithm R times allows P^^ to be expressed in terms of (see (2.27)). This procedure replaces (2.24) in the data deletion step. Substituting (2.18) into the first identity in (2.22) yields PN-R6R = W1 " WP"' <2'30> which allows (2.23) to be expressed in the form ^N-R " ^N + PN-R6R(VN " V . (2'31) The data deletion is performed using (2.29) and (2.31). No matrix in version is required. Insertion of new data is performed using (2.29) in reverse and by changing (2.31) to the appropriate form. 2.4 Measurement Requirements 1 Due to the use of state variable filters not all states need be measured as implied by (2.2). Most systems of the form (2.1) can be decomposed into subsystems which have a phase variable structure. Suppose, for example, that a subsystem is represented by x3 = a3x3 + a2x2 + alX]L + ... x2 = X3 xx = x2 (2.32) To identify (2.32) it is sufficient to generate the complete states of one the augmented models. The first augmented model that can be generated completely is the third model (Jl = 3) . This is seen by taking Z3 = It is seen from (2.33) that only x^ need be measured. Any additional augmented model (l > 3) can be generated using (2.3) and (2.33). C xi dTi dx2 ^ ^ x„ dr, dx0 o o 2 1 2 ft x„ dr, dx0 o o 3 1 2 ^ fZ x dr1 dx„ o o 1 1 2 •^t x1 dx, oil (2.33) 2.5 Applications Example 1. Consider the linear system x.^ = a-|X2 + a2x1 + u, x^(o) X2 Xl , x2(°) 1.0 0.0 (2.34) where a^ = -2.0, a^ = -3.0, tQ = 0, t^ = 16., and take I = 4 (see The fact that only two elements need be identified results in some simplification. The full augmented system is of the form xl "X2 xl" "Ul" x3 X4 x3 "al" U2 = • + X5 X6 x5 U3 _x8 X7 _U4 _ (2.35) where x^ = x^, xg = x5» x8 = x7' Equation (2.35) is discretized and expressed in the form of (2.8) by appropriate definition of Z^(k), ^(k) an& c' = (ap » The estimation of c is performed by a block processing mode with R = 1. (if no old data is deleted, the estimation algorithm is essentially the Kalman filter algorithm). The following sample intervals and control are used: (1) T = 0.1, (2) T = 0.2, (3) T = 0.4, (4) T = 1.0 with the corresponding values of N: (1) N = 21, (2) N = 11, (3) N = 6, t (4) N = 3. The control u(t) = 0.1 + sin 3t + sin lOt is taken. Figures (2.1) and (2.2) give the results with zero additive measurement noise (n^ = o = n2) . However, discretization of (2.35) introduces quantization errors which can be considered as additive correlated quantization noise. It is seen from the figures that identification is rapid even for large quantization noise. Figure (2.3) and (2.4) illustrate the results for state dependent noise of the form x,(k) n.(k) = a. • (-i ) ng., (i = 1, 2) (2.36) m where ng^ is a pseudo-random computer generated noise sequence which has a maximum amplitude of n . The choice a. = a„ = 0.1, E(ng.) =0. is m 12 l made. Figures (2.5) and (2.6) illustrate the results for the choice = 0.1, a2 = 0.25, E (ng^) = 0.1, E (ng2) =0.14. It is seen from these results that the estimator is insensitive to correlated quanti zation noise, state dependent measurement noise and to bias in the measure ment noise. Example 2. Consider the linear system (See Reference [3]) Xi = x2 , x 3^(0) = 1.0 X2 = alXl + a2x2 + bjU , X 2(0) = 0.0 where = --2.0, a2 = -1.0, b1 1.0, and take Z = 4. system is of the form • X2 xl x2 ul X4 x3 x4 "al" U2 X6 = X5 x6 _a2_ + U3 bl uX8, X7 X8 U -^U4 _ • • • (2.37) where x_ = x,, xr = x,, x_ = xn. 3 4 5 6' 7 8 The following sampling intervals are used: (1) T = 0.1, (2)- T =• 1.0-. The-cor-responding- values of N are: (1) (1) N = 21, (2) N = 3. The value of R is taken as R =1. The initial parameter estimate is c = 0 and u(t) = 0.1 + sint + sin 3t. Figures (2.7), (2.8) and (2.9) give the results for quantization noise only. Figures (2.10), (2.11) and (2.12) give the results for quantization noise and additive state dependent noise where «1 = a2 = 0.1, E^g^ =0, i = 1, 2. Comparing these time responses with those obtained in Reference [3] for the noise-free case, it is evident that identification is more rapid. The speed of response of parameter tracking systems seems to be adversly effected by the interaction between parameter generating subsystems which is an inherent feature of the steepest descent search procedure used. Such subsystems are not required for the method developed here. Consequently, in the noise free case, one-step identification is achieved. FIG. (2.2) EXAMPLE 1. ESTIMATION OF WITH QUANTIZATION NOISE 0.4 FIG. (2.4) EXAMPLE 1. ESTIMATION OF a2 WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE. 0.0 t 2 4 6 8 10 12 U 16 IS 20 FIG. (2.61 EXAMPLE!. ESTIMATION OF a2 WITH QUANTIZATION AND BIASED . MEASUREMENT NOISE^ 0.4r 1 1 1 1 1 1 2 3 4 5 5 7 0 9 70 ®-\ N \ ' V FIG.(2.7) EXAMPLE 2. ESTIMATION OF a, WITH QUANTIZATION NOISE. -i 1 1 1 1 1 1 i 1 1— 123456789 10 FIG. (2.8) EXAMPLE 2. ESTIMATION OF a2 WITH QUANTIZATION NOISE. 1.6 1.2 0.4 0.0 or* I 2 3 10 -0.4L FIG. (2.9) EXAMPLE 2. ESTIMATION OF b. WITH QUANTIZATION NOISE. FIG. (2.10) EXAMPLE 2. ESTIMATION OF a, WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE. 04. 2.0 . -2.4 . FIG. (2.11) EXAMPLE 2. ESTIMATION OF WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE. 2.4 2.0 1-6 1.2 0.8 0.4 0.0 Or** w FIG. (2.12) EXAMPLE 2. ESTIMATION OF b, WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE. 27 2.6 Comparison with a Limited Memory Filter In this section a comparison is made between the proposed digital identifier using augmented states (I > 1) and the standard limited memory identifier ( I = 1). The standard limited memory identi fier (or filter) is of the type discussed in Section 2.3. The filter operates in a block data processing mode where R old data blocks are deleted and R new data blocks introduced resulting in a limited (finite) memory filter. Example 2 is used to illustrate the superiority of the augmented state identifier for various values of T and N (R = 1). Figures(2.13), (2.14) and (2.15) show the results for the identification of the parameters a^, a^ and b^, respectively. It is evident that the standard limited memory filter fails to track the parameters for values T = 1., N = 3, R=l and exhibit poor tracking accuracy for values T = 0.1, N = 21, R = 1. The failure of the standard limited memory filter to give acceptable estimation of the parameters is due to the fact that it is not an optimal filter for the proposed identification problem. Consequently, excessive quantization noise and correlation between samples, as is the case here, can have a serious effect on its performance. In conclusion it can be stated that the use of augmented states and controls results in a more rapid and accurate identification of continuous time linear system parameters than can be achieved by either parameter tracking systems or the standard limited memory filter. FIG. (2.13) EXAMPLE 2. ESTIMATION OF a, WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE USING A STANDARD LIMITED MEMORY FILTER. FIG (2.14) EXAMPLE 2. ESTIMATION OF o2 WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE USING A STANDARD LIMITED MEMORY FILTER. 29 28. 2.4 . 2.0. 16. FIG. (2.15) EXAMPLE 2. ESTIMATION OF b, WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE USING A STANDARD LIMITED MEMORY FILTER. 30 2.7 Some Statistical Properties of the Digital Identifier In the previous section, the performance of the digital identi fier was evaluated by comparing estimated values with known values. In a practical situation, however, the parameters are not known and statistical means must be used to establish confidence in the estimates. In reference [5] Kishi made use of extensive studies made by Linnik [17] to determine confidence intervals for parameter estimates in the case of a simple single-input single-output system. The relations given by Linnik will be used in this section to obtain approximate confidence limits for the parameter estimates of the digital identifier. However, Linnik's results can be applied only if the different components of the noise vector N (see (2.10)) can be represented by a statistical noise model. Such a task is not an easy one. However, it will be shown that under certain simplifying assumptions it can be argued that the different noise com ponents will be uncorrelated. 2.7.1 Correlation of Noise Samples In Section (2.2), it has been assumed that the noise samples of the different models are uncorrelated. Hence, equation (2.16) is assumed to be true. In the following, it will be shown that if the noise samples are uncorrelated then the samples of the integrated noise are also un correlated. For simplicity, we will consider a single noise component in order to prove the validity of the previous statement. Consider the noise component n(t) which can be represented by the following statistical model E[n(t)] = 0 (2.38) E[n(t)n(x)] = ae-01!^' If a is of infinite value then n(t) will be a type of white-noise com ponent. (In the white-noise case the amplitude a will be equal to in finity) . For a physical realization of the noise component n(t) it will be assumed that a is a very large but finite quantity and that a' is a very large but finite quantity and that a is of finite value. Let Nx(t) = n(T1)dx1 (2.39) It follows from (2.38) and (2.39) that E[N1(t)n(T)] = E[n(x)n(T1)]dx1 1 o 1= a /t e-alT-Tl' dx, (2.40) o 1 To evaluate the above integral we will consider two cases: (i) Case 1, x ^ t E[N1(t)n(x)] = a £ e~a(T_Tl) dx^ = a [e-a(x-t) _ e-ax a (ii) Case 2, o sc x «_ t T-TUT / \ i rT -a(x-x.) , , ,t -a(x..-x) , E[N,(t)n(x)J = a J e 1 dr, + a f e 1 dx, 1 o 1 x 1 a -axr ax , a axr -ax -at, = — e [e — 1J H e [e -e J a a a ro -ax -a(t-x), ,0 = — [2-e-e] (2.42) a 32 For fixed time t and a very large value of a, equations (2.41) and (2.42) can be approximated by f = 0 , x > t I (t)n(x)] J U^,0<x< E[NlVW'n"J ** (2.43) As seen from (2.43), the noise components n(t) and N^(t) become un correlated as a -*- 00. Using equation (2.41) or (2.42) for the special case where x = t yields E[N1(t)n(t)] = f [1 - e"at] = — , for fixed t and a ->• °° The following non-dimensional ratio can be used as a measure of the degree of correlation of N^(t) and n(t): E.IMr(t)q(.t)] r = -r , for fixed t, a -> « (2.44) at Equation (2.44) shows that for fixed t and a large value of a the correlation between N^(t) and n(t) samples is much less than the correlation between the samples of n(t). This shows clearly the advantage of using the inte grators in order to generate the augmented states and controls. It is also interesting to study the degree of correlation between the samples of Nx(t). Equation (2.39) yields E[N1(t)N1(x)] = E[N1(t)n(x1)]dx1 (2.45) As before, we have two cases: (i) Case 1, x > t ElN.WN.d)] = f]2 - e aTl - e~a(T"Tl)] dx, 1 1 o a 1 -ax, -a(x-x,) . . « [2Tl I - £ a 1 a a o a ro -ax , -at -a(x-t) = —j [2at +e +e - e -1] a 2at = , for fixed t and a -*• 00 (2.46) a (ii) Case 2, o < x < t E[N (t)N (x)] = /T E[N (x)n(x1)] dx, + / E[N (x)n(x.)] dx, li ol llxl 11 fx a ro -ax -a(x-x1)1 , , = J — L2 - e 1 - e 1] dx, + o a 1 rt a r -a(x,-x) -ax, -, , / — Le 1 - e 1] dx, x aan ~ax , ~at -a(t-x) = —[2ax + e +e - e -1] a 2ax = , for fixed t and a •> 00 (2.47) For the case x = t equation (2.46) or (2.47) yields E[N-2(t)] =~ [at + e at - 1] a = — t (2.48) a Equations (2.48) and (2.38) yield E[N2(t)] —4 al (2.49) tE[n^(t)] a Equations (2.46), (2.47), (2.48) and (2.49) show that the correlation between the samples of the noise component N^(t) vanishes as a -»• 00. Equation (2.49) shows that for fixed t and large value of a the corre lation between the samples of N^(t) is much less than that between the samples of n(t). The same analysis can be directly extended to ^(t) where N2(t) = N1(x1) dx1 (2.50) Since the same approach is to be used it is sufficient to summarize the corresponding results E[N2(t)N1(t)] = ^ t2 (2.51) 2 at3 E[N2(t)] ~ fo~ (2"52) Efr-Cty^ct)]-, 9 = (2.53) t3E[n2(t)] E[N2(t)] ~ 1 (2.54) t4E[n2(t)] 3at Equations (2.51), (2.52), (2.53) and (2.54) show that the same conclusions that have been reached for N^(t) hold for N2(t). Suc cessive integrations of statistically independent noise samples result in statistically independent noise components. Also, successive inte grations of correlated noise samples result in less correlated noise components. The above simplified mathematical treatment shows that the independence property assumed by equation (2.16) is justifiable for many practical situations where n(t) is usually taken as white noise. Equation (2.16) can now be used to prove the unbiasness property of the digital identifier. Equation (2.10) has the form ZN(k + 1) = 6N(k)c + N(k) (2.55) Where the argument value denotes the sampling instant at which the various vectors and matrices are evaluated. Combining (2.55) and (2.12): E[c] = E[PN(k)9^(k)ZN(k + 1)] = E[PN(k)e^(k)(GN(k)c + N(k))] = c + E[PN(k)e^(k)N(k)] (2.56) The elements of PM(k) and 8 (k) depend on Z (k) and U (k), A = 1, n+m. Equation (2,55) shows that Z^(k) is only related to N(k - 1). Since N(k) is statistically independent of N(k - 1) as shown by the above analysis, then ^(k) and 6jj(k) are not statistically dependent on N(k) and the following relations hold E[PN(k)8^(k)N(k)] = 0 (2.57) E[£] = c (2.58) Equation (2.58) shows that the estimate obtained by using the digital identifier is unbiased. The results given in figures (2.1)-(2.4) and (2.7)-(2.12) would not be obtainable if (2.16) did not hold approximately. Equation (2.58) holds when 7^ = (9^9^) ^. This case has been referred to by Kishi as the observable case. The case where the inverse (9^9^) ^ does not exist is referred to as the non-observable case. In the non-observable case an estimate is obtained by use of the equation c - 9* ZN (2.59) where is the generalized inverse of 6^. That is, e^eNx = X (2.60) It follows from (2.16) and (2.60) that the estimator (2.59) is unbiased. 2.7.2 The Error Covariance Matrix Let V = 6Na - ^ (2.61) be the error vector representing the difference between a predicted measurement vector and the measurement vector. For the observable case the above equation can be written in the form (see (2.12)) V = <Wi - I)ZN (2'62) Using equation (2.55) yields (ENPN6N " I)S (2'63) For the non-observable case the expression for V is v = (eNe^ - I)N (2.64) where use is made of the property N N N N The error convariance is defined by 6*B = fl (2.65) Qw = E[W»] (2.66) Making use of equation (2.63) yields QW = E[(WN " I)M'(9NPNEN " I)] (2*67) For the non-observable case using (2.64) yields QW = E[(6N6N " ^'(VN " I)] (2'68) 2.7.3 Estimation of the Confidence-Interval Length Only the observable case will be studied in this subsection since for the non-observable case the length of the confidence interval will depend on the degree of redundancy in the measurements which, would change from one system to the other. Equation (2.63) yields V = (6NPN6N " I)S (2>63) The fol*lowing relations' hold for the observable case:" R(8N) = R(6^) = R(6^6N) = R(PN) = n(n + m) (2.69) where R(*) stands for the rank of the matrix written between the brackets. Let the matrix Y„ be defined as N YN - PN6N (2-70) Using the known lemma R(AB) 4 min(R(A), R(B)) (2.71) where A and B are any arbitrary matrices, it follows from (2.70) that, R<V = n(n + m) (2.72) Equation (2.70) yields 3N = (9N9N)YN (2'73) Hence R(8jJ) = n(n + m) < min(R(6^9N) , R(Y^)) < min(n(n.+m), R(Y )) (2.74) Hence R(YN) > n(n + m) (2.75) It follows from (2.72) and (2.75) that R(YN) = R(Pn<3N) = N<N + M> <2'76> Let the matrix be defined as N "HHW (2-77) Then WM9M = 9.. (2.78) N N N It follows from (2.78) that is a projection matrix. Using (2.78) and (2.71) yield n(n + m) = R(8N) < min(R(WN) , R(e^) Hence. <=min(R(WN), n(n + m)) (2.79) R(WN) ^ n(n + m) (2.80) Applying (2.71) to (2.77) it is seen that R(WN) 4 n(n + m) (2.81) It follows from (2.80) and (2.81) that R(WN) = n(n + m) (2.82) Considering now the computation of the number of degrees of freedom of the random vectors c and V, that is (as defined by Linnik [17]) the number of statistically independent components of c and V. Starting by computing the degrees of freedom of the random vector c we will make use of Theorem 2.3.1 given by Linnik [17] which states that if the vector U is expressed as and the matrix A is of dimension r * £ with r < Z and R(A) = r then the random vector U is of r degrees of freedom. From equation (2.12) £ = PM9'ZM (2 N N N Equation (2.76) yields Applying Theorem 2.3.1 given by Linnik [17] it is readly seen that the vector c is a random vector of n(n + m) degrees of freedom. To compute the degrees of freedom of the noise error vector V, (2.63) is expressed in the form U = AX (2.83) R(P.T0') = n(n + m) (2.84) V = (WM - I)N (2.85) where W is given by (2.77). Let 'N (2.86) The matrix is idempotent as can be seen from the following = W - 2W + I N N = I - N (2.87) Accordingly, there exists an orthogonal matrix 0^ such that (see reference [17]) where the matrix is a diagonal one. The matrix $^ has the property that $N*N " VN = 1 (2'89) The matrix is an idempotent matrix as can be seen from the following N N N N N N N N N N =z $'w $ = n (2.90) Hence: d2. = d.. i.e., d.. =0 or 1 (2.91) n ii n The number of d^s equals to unity is given by the rank of the matrix W^, that is, by n(n + m). From (2.89) it follows that R(I - WN) = R(^(I - WN)$N) = Rd - y = n(n + m)(N - 1) (2.92) The final result follows from the fact that the total number of diagonal elements of is n(n + m)N and of these n(.n + m) are equal to unity (see (2.91)). Applying the above result given by (2.92) to (2.85) and making use of Theorem 2.3.1 [17], it follows that the noise error vector V is a random vector of n(n + m)(N - 1) degrees of freedom. From the result reached in the previous statement it is seen that V'V can be expressed in the form n(n+m)(N-1) V'V = I (2.93) i=l Let QN = E[c c'] (2.94) cc With c is computed according to (2.12). The following results given by Kishi [5] and Linnik [17] are required. If £ and E,^ are statistically independent Gaussian random 2 2 — variables with distributed as % and having m degrees of freedom, then the t-distribution is formed by the following ratio t = g (2.95) m Let (see Appendix IV) c. -c. K = 1 1 (2.96) / (Q ).. CC IX where (Q ).. is the i-th element of the main diagonal of the covariance cc 11 matrix Q and c. stands for the i-th element of the vector c. Also cc 1 let (see Appendix IV) m = n(n + m)(N - 1) .2 £q = V'V (2.97) Then using (2.96), (2.97) into (2.95) yields c. - c. 1 1 (2.98) "n(n + m)(N - 1) ^ )>>(v,v/n(n + m)(N _ 1} cc 11 42 Using the t-distribution it is possible to compute the length of the interval about the parameter c. which would include c. with a certain i 1 probability. For example, let the assigned probability be 0.9, then ^nCn + DCN - 1)1 4Y> = 0.9 (2.99) where P{*} stands for the probability density of the expression between the brackets. The number y can be found from tabulated standard tables using the number n(n + m)(N - 1) as an index. Hence Ic. - c.I = Y/(Q ).. , , VsYM TT (2.100) 1 l x1 1 cc ii n(n + m)(N - 1) If the augmented states and controls are not used the denominator under the square-root sign in (2.100) will be n(N - n - m). Thus using the augmented states and controls decreases the estimation error | c\ - c^| and results in a better estimation accuracy. The range 2A of the confidence interval for a probability of 0.9 is given by 2A = [c. ± YAQ ).. , . VZ TV (2.101) l ' %cc ii n(n + m)(N - 1) 2.8 Discussion For the observable case equation (2.12) is used to estimate the unknown parameters. The matrix P^ is recursively computed without any need to invert a matrix. For the non-observable case the matrix P^ does not exist. Equation (2.59) is to be used instead of (2.12) in order to estimate the unknown parameters. The generalized matrix inverse 6^ is required. It is desirable to compute G^j in a way similar to that developed in this Chapter for computing P^ in the observable case. This requires a computation algorithm to evaluate 0^ recursively without any need for a matrix inversion. In the next chapter an algorithm will be developed in order to compute 6* recursively and no matrix inversion will be needed. This algorithm is computationally simple and is suitable for real-time compu tations . A simplified version of the algorithm that will be developed in the next chapter can be used in order to compute P^. This presents us with another method to compute P^. Due to the simplicity of this algorithm it can be used for state estimation applications in cases where the plant is represented by a discrete model. Since the generalized inverse has numerous applications besides those which occur in state and parameter estimation, the important topic of recursive computation of a generalized matrix inverse is dealt with in the next chapter in a general way. 3. RECURSIVE COMPUTATION OF THE GENERALIZED INVERSE OF A MATRIX 3.1 Introduction Because of its value in the solution of minimum norm problems in finite dimensional spaces, the generalized inverse (or pseudo-inverse) finds important applications in linear system theory. Applications of the generalized inverse to estimation and filtering problems are given in [18], [6], [19]. Applications to control problems are given in [20], [21]. A number of methods have been proposed for computing the generalized inverse A''" of a matrix A. The method proposed by Ben-Israel et. [22] al. requires the formation of C = A*A, where A* is the conjugate transpose of A, and a series expansion of C. The method proposed by Pyle [7] requires a two-stage Gram-Schmidt orthogonalization process and subsequent solution of sets of linear equations. A simpler algorithm has been pro posed by Rust et. al. [8]. The algorithm is based on a two stage Gram-Schmidt orthogonalization process. Unlike the other two methods men tioned above, it does not require the generation of powers of matrices or the solution of sets of linear equations. However, by eliminating the Gram-Schmidt orthogonalization process, an even simpler recursive algorithm can be derived. 3.2 Properties of the Generalized Inverse Let A be an m x n complex matrix. The generalized inverse of A, as defined by Penrose [23], is the unique n x m atrix A* which satisfies the following relations: AAIA (AA1) * (A1 A)* = A = AA AXA (3.1) (3.2) (3.3) (3.4) Following Rust et. al., suppose that A can be partitioned into the form A = (R,S) (3,5) where R is an m x n^ matrix with linearly independent columns and S is an m x n2 matrix whose columns are linearly dependent on the columns of R. The following relations hold: RIR = I R1 = (R*R)-1R* S -RU U = R S (I + UU*) 1RI U*(I + uu*)"^1 (3.6) (3.7) (3.8) (3.9) (3.10) (In (3.6) and (3.10) I is a unit n^ x n^ matrix). The generalized inverse is given in partitioned form by (3.10) where ^ and U are given by (3.7) and (3.9), respectively. Rust et. al. show that the inverses in (3.10) can be found by a two stage Gram-Schmidt orthogonalization process (the first stage generates R*)• The possibility of replacing both orthogonalization stages by a recursive algorithm is considered in the next two sections. 3.3 Recursive Evaluation of R The following relationships are required to develop the recursive algorithm. Taking the conjugate transpose of (3.1) and using (3.3) yields Let A*AA = A* A = (V,W) (3.11) (3.12) be a partitioning of A where V is composed of linearly independent columns. Equations (3.7) and (3.11) can be applied to V: v1 = (V*V)-1V* v*vv = V* (3.13) (3.14) From (3.1) and the partitioning (3.12) it is seen that AA V = V (3.15) Taking the conjugate transpose of (3.15) and using (3.3) (3.13) and (3.14) yields V^1 = (V*V)-1V* = V1 (3.16) Let r_., (j = 1, n^) denote the linearly independent columns of R. Consider the partitioning k+1 k+1 (3.17) where R, = r,, R = R, and where the b.'s are row vectors. From (3.17) 1 1 n1 J it follows that " W+i + rk+ibk+i (3.18) Applying (3.6) to yields = I (3.19) Applying (3.11) to 8^+1 111 the partitioned form (3.17) yields k+1 r* k+1 (3.20) Partitioning of (3.20) and solving for b^-fi gives -1 \+i= <WW rk+i(I " Vk+i> (3.21) Replacing A by and V by R^ in (3.16) gives ^W^+l " (3.22) Substituting (3.18) into (3.22) and using (3.19) yields w = V1 - WW' (3.23) Substituting (3.23) into (3.21) gives W = wW1 - \v (3.24) where °k+l = (rk+lrk+l " rk+l\<rk+l)_1 (3.25) is a scalar. The recursive algorithm for the computation of R is based on (3.17) (3.23) (3.24) and (3.25). It is desirable, however, to compute ,1 • the matrix product R^^k. in (3.24) from a recursive formula. Substituting (3.23) and (3.24) into (3.18) yields = Vk"+ vi(i - \^)rk+it+i(i - v£ (3-26) Since (see (3.3)) (3.27) it is seen that the second matrix on the right-hand side of (3.26) is the outer product of two vectors Summary of the Algorithm for Finding R"*" Initialize. Set k = 1, R][ = r , R* = (R*R )-1R* \+i" (rk+irk+i ~ t+A^W-1 bk+l = \+lrk+l(I " \RJ)  Bk+1 " ^(I " VlW vAi = \*i+ vi(I -A^wW1 -<+1 3.. k. = k+1 k+1 k+1 4. If k < n1 go to 2, The above algorithm replaces the first Gram-Schmidt orthogonali zation stage of Rust's algorithm, which, in the interest of accuracy, reorthogonalizes each column after it is first orthogonalized. An addi tional advantage of the above algorithm is the manner in which A is partitioned into the form (3.5). In the Gram-Schmidt process, the de-49 tection of a zero vector after orthogonalization (but before normali zation) indicates that the vector being orthogonalized is linearly dependent on the previous vectors. This vector is consequently moved into S and the next column vector of A is chosen. In the algorithm presented here, linear dependency is detected when CL ^ (see (3.25)) is K."t"-L zero. Since a^+^ is required in (3.24) no additional computations are required to test for linear dependency. If 0^+1 ^s zero» fc^e associated vector is moved into S and the next column vector of A is chosen. 3.4 Recursive Evaluation of A Consider the partitioning \+l. " (R»Sk+l?' \+l " k+1 k+1 (3.28) where S. = sn is the first column vector of S and S = S. The generalized 11 n? I z inverse R can be obtained by the recursive algorithm given in Section (3.3). Applying (3.10) to A^, it is seen that J 1 C. = U*B. J 11 (3.29) (3.30) U. = (U]L, u2, , u.) = R S. . J J (3.31) where D. = (I + U.U*) 1 11 (3.32) The matrix inversion lemma [9] \li - \' - Dk\+i(I + "kVA'W^A1 (3.33) is used to develop a recursive formula for B^. Setting j = k, k + 1 in (3.29) and using (3.33) gives where Bk+i= (I" ViVWk+i^k ek+1 = .(i + ug+1BkR uk+1)" (3.34) (3.35) The recursive algorithm for the computation of A* given is based on (3.28) (3.30) (3.31) (3.34) and (3.35). • Summary of the Algorithm for finding A* 1. Initialize. Set k = 1, B]. = D'V = [I - (1 + u*u1)~1u1u*]RI -1 Vi • (I - ek+i\R \+it+i)Bk 3. 4. 5. k = k + 1 If k ^ n^, go to 2, C = (RIS )*B n2 n2 n2 A1 = A1 nv. 3.5 Examples The first two examples are taken from the literature. This allows a comparison to be made showing the relative simplicity of the proposed algorithm. Example 1. This example is given in [22]. It is required to compute the generalized inverse of A = It is seen that (3.36) has the form (3.5) where 1 0 0 -1 -1 1 0 0 0 -1 1 0 0 0 -1 1 (3.36) R = 1 0 0 -1 1 0 0 -1 1 0 0 -1 s = (3.37) The algorithm of Section 3.3 is used to find R Step 1. Initialization. Step 2, Rl " rl 1 -1 0 0 a2 = 3 R* = |(1, -1, 0, 0) b0 = ±(1, 1, -2, 0) B, |(2, -1, -1, 0) R2R2 1 3 2 -1 -1 0 -1 2 -1 0 -1 0 -1 0 2 0 0 0 Step 2. (repeated) b„ = f (1, 1, 1, -3) = 1 f3 -1 -1 -l) "4 [2 2-2 -2J R = R = 3 -1 -1 -1 2 2-2-2 111-3 Having found R*, the algorithm of Section (3.4) is used to evaluate A*. Since S consists of a single column vector s^, only the initialization stage is required. ui " RISI = -1 -1 -1 Step 1. •i-i 3-3-1 1 1 3-3-1 -11 3-3 Step 5 Cl = UPl - | ( 3, 1, -1, -3) A1 = 1 Al 8 3 -3 1 3 -1 1 -3 -1 -1 1 -3 -1 3 -3 1 3 (3.38) Example 2. This example is given in [7]. It is required to compute the generalized inverse of 10 11 0 1-10 11 0 1 (3.39) It is seen that (3.39) has the form (3.5) where 53 1 0 1 l" R = 0 1 , s = -1 0 1 1 0 1 (3.40) The algorithm of Section (3.3) is used to compute R*. Step 1. Initialization Step 2. Rl " rl , RJ = \ (1, 0, 1) 0U = 2 3 f (-1, 2, 1) B2 = R1 = f (2, -1, 1) 1 = 1 , 2 -1 1 2 3 K-l 2 V Having found R^, the algorithm of Section (3.4) is used to evaluate A*, Step 1. Initialization U = (up u2) D"1 = I - ul(l + uju^u* 3 vl 2J 1 3 [0 1 IJ Step 2. B = -i- ( 3 ° 3) 2 15 v-l 5 h} Step 5. C2 = U*B2 - T5 (3 "o A2 " 15 3 -1 4 3 0 5 -5 0 3 4 -1 3 (3.41) Example 3. This example illustrates a case where A does not have the form (3.5) initially. It is required to compute the generalized inverse of A = 1 0 1 0 -1 1 0 0 0 -1 -1 0 0 0 0 1 (3.42) Note that the third column of A is the sum of the first two columns. Generally, the detection of linear dependency by simple observation is not possible. The algorithm of Section (3.3) is applied to (3.42). It is seen from (3.36) that there is no change in the computation given in Example 1 until k = 2. At this stage the algorithm gives = 0, hence the third column of A is shifted into S and the last column of A is chosen. This gives = 1 and b3 = (0, 0, 0, 1) 1 (2 -1 -1 0 3 3 1 1 -2 cr 2-1-10 1 1-2 0 0 0 0 3 The interchange of the last two columns of (3.42) is accomplished by a permutation matrix P given by 55 10 0 0 0 10 0 0 0 0 1 0 0 10 (3.43) The algorithm of Section (3.4) is now used to find A*' where A = AP, A = AP (3.44) The generalized inverse A''" is then given by (see [8]) A1 = PA1 (3.45) Step 1. Initialization u^ = R s^ = 1 -1 0 0 0 1 -1 0 0 0 0 3 Step 5. U1B1 = I (1» °» _1' 0) A1 1 3 -1 1 0 0 0 0 -1 0 0 3 -1 0 (3.46) The generalized- inverse A-^ is- found- from (-3-. 45) and- is-A = -1 1 0 0 0 0 -1 0 1 0 0 3 (3.47) 4. TWO-STAGE ESTIMATION OF TIME-INVARIANT AND TIME-VARYING PARAMETERS IN LINEAR SYSTEMS 4.1 Introduction Because of its importance in system engineering, considerable effort has been devoted to the problem of estimating the parameters of systems whose dynamic behaviour can be described by differential equations. Two distinct methods for estimating time-invariant parameters have been brought to a high degree of refinement. The parameter tracking method is based on a steepest-descent parameter adjustment of a quadratic error function or a Liapunov function. The error can be taken to be a response error or a generalized equation error. Detailed descriptions of para meter tracking systems are given by Lion [2], Pazera and Pottinger [3], Rose and Lona [4]. The second method for parameter estimation is based on optimum filtering. In its simplest form an optimum filter minimizes a mean-square error cost function subject to a prescribed dynamical constraint. Detailed descriptions of optimum filters and an extensive bibliography are given by Sage and Melsa [9]. The estimation of time-varying parameters is generally a much more difficult problem. Parameter tracking systems and optimum filters often give acceptable estimates provided that the parameters vary suffi ciently slowly with respect to time. Lessing and Crane [10] have made use of integral transforms to extend the capability of parameter tracking systems to track time-varying parameters. Practical implementation of the integral transform approach requires the assumption that the parameters can be represented by polynomials of known order. The estimation of time-varying parameters by optimum filters requires essentially the same assumption. It is assumed that the parameters are solutions of differential equations of known form with unknown time-invariant para meters. The problem is then reduced to the estimation of time-invariant parameters of a higher order augmented system. It is the purpose of this Chapter to develop a computationally simple and accurate method for estimating time-varying parameters which is not restricted by a priori assumptions about the dynamical behaviour of the parameters. 4.2 Problem Formulation Consider the continuous-time linear system x = Ax + Bu, x (t ) = x (4.1) p p p o po where x is an n x 1 dimensional state variable vector and u is an m x 1 P dimensional control vector which is a known function of time. The system matrix A is considered to be decomposable into a time-invariant component A^ and a time-varying component A^: A =• Af + Av (-4.2) The most general estimation problem associated with (4.1) is to estimate x , A and B given an observed continuous-time measurement vector defined P over an observation interval t 4 t <_ t^, and given complete statistical information about measurement noise and system noise. However, if emphasis is placed on computationally simple methods then a less ambitious problem must be considered. Suppose that B is known and that the state x^ is 58 measured during the observation interval. The problem to be considered is the estimation of A given x^ and given no information about noise statistics. The estimation is performed sequentially in two stages. During the first stage an estimate of the time-invariant A^ is made given an estimate of A^. During the second stage a new estimate of A^ is made using the estimate of A^ found during the first stage. This section deals with the second stage of identification. Let z. = x - x, (4.3) 1 p 1 x, = A.x. + Bu, x.(t ) = x (4.4) 1 f 1 1 o po It follows from equations (4.1), (4.3) and (4.4) that z., = A.Z.. + A x 1 f 1 v p = A.z, .+ 6(x )a, z,(t ) = 0 (4.5) t l p — 1 o In equation (4.5), a is a vector whose components are the time-varying elements of A . v The equation A x = 0(x )a (4.6) v p p -defines a matrix 0 whose elements are linear functions of x . Let P g = xp - xx h A g(t) -Tg(t-T) (4>y) In (4.7), x^ is defined by (4.4), x^ is assumed known (through measurements) and h. is an approximation to g, with the accuracy of the approximation depending on the fixed-time increment T. Consider the following quadratic cost index J=jflfo [||g-^1l|2Q1+ llh- ^||"2Q2]dT (4.8) where z^ and are defined by (4.5) where the unknown vector a. be re placed by an estimate a_. The weighting matrices and are taken to be positive definite. It follows from (4.3)-(4.8) that lim J = 0 pro-T->0 vided that the time-varying parameter vector ci is known exactly. Con sequently, (4.5) and (4.8) can be used to formulate a problem in optimum estimation. In the problem formulation the unknown time-varying para meter vector is considered to be a control vector for the system des cribed by (4.5). The optimum control minimizes J where g and h are known functions of time. 4.3 Estimation of A^. Algorithm I. The minimization of J (see (4.8)) subject to the dynamical con straint given by (4.5) is a standard problem in optimal control. The Hamiltonian for this problem is (prime denotes transposition). H = p'CA^ + 6(x )a)- -|(g-z1)'Q1(g-z1)- j(h-Af z±-Q (x ) a) ' Q2(h-Af z.j-6 (x ) a) (4.9) The costate equation is p = -Hz = -A^p + Q1(g-z1) - A^Q2(h-AfZ;L - 6(xp)a), p(tf) = 0 (4.10) and the gradient condition is 60 H = 9' (x )[p+ Q,(h - A,z - G(x )a)] = 0 (4.11) a p L r 1 p — Solving (4.11) for the optimum control a yields a = (9»(x )Q29(x )) V(xp)[p + Q2(h - Afz^ ] (4.12) Existence of an optimal control for the formulated problem can be viewed as a controllability condition in a function space S of elements (g, h), where h is either equal to g or an approximation to g. Taking (4.8) to define the norm, then the system described by (4.5) is defined to be controllable in S if a control ji exists such that the element (z^ - g, z^ - h) has minimum norm. Satisfaction of this con trollability condition requires the existence of the matrix inverse in (4.12). Equations (4.5), (4.10) and (4.12) constitutes a two-point boundary-value problem. Due to the linearity of the equations the solution is easily found using standard methods. Let L(xp) ^ (6'(xp)Q26(xp)) V(xp) (4.13) Substituting (4.12) into (4.5) and (4.10) it is seen that where C. I 9(x )L -r' 4 + r (4.14 C± = (I - 6(x ) LQ2) Af C2 = -Ql + A'Q^ 9(xp)LQ2h Q1§ - A^Q2(I - 0(x ) LQJh (4.15) 61 and where I is the unit matrix. The general solution of (4.14) has the form where P(tf) - *(tf, tQ) + p(tc) y2 A y = Ly2 J A tf = / 1 iKt, T) r(x)dx fco (4.16) (4.17) and where i|;(t, x) is the state transition matrix for (4.14). It is seen from (4.16) that y can be found by solving (4.14) taking z^(tQ) = 0 = p(tQ) for the initial values. The state transition matrix is found by solving (4.14) taking -Kt , t ) = I, r = 0. Substituting the boundary conditions into (4.16) and parti tioning ^(t^,, t ) in the form it is seen that *(tf, tQ) = "*11 |*12 (4.18) *21 j*22 p(tQ) - .-1 -*22 y2 (4.19) The initial costate given-by (4.19) in conjunction with (4.5), (4.10)' and (4.12) defines the optimal estimate of a. over the observation interval tQ <^ t ^ t^. The computational procedure will be referred to as algorithm I. Minimization of a cost index of the form given by (4.8) subject to the dynamical constraint of (4.4) has found extensive applications in state estimation (Sage [11]). There is, however, a fundamental difference in the problem formulations. In the state estimation formulation, (4.5) and (4.8) are replaced by x = Ax + Bu, x (t ) = x (4.20) P P P o po ~J = 1 t/tf [ll*p - *PM2q1+ ll*p - %\\2^ <4-21> The minimization of J subject to the dynamical constraint of (4.20) is a problem in optimum filtering. The problem can be considered as one of finding an approximation x to z so that the error (z - x ) has minimum P P P P norm subject to a cost penalty associated with the norm of the control effort. The solution of the two-point boundary-value problem defined by (4.20) and (4.21) can be obtained by invariant imbedding or by use of a Riccati transformation. Either approach results in the Kalman filter algorithm, for generating, the. optimal, estimate, x^ of the state. xp given. the measurements z over an observation interval t < T < t. Comparing p o = = the second terms in the integrands of (4.8) and (4.21) it is seen that there is a fundamental difference between the state and parameter esti mation problems. This difference manifests itself analytically in the controllability condition (the existence of L(xp), see (4.13)) which does not arise in the solution of the state estimation problem. 4.4 Estimation of a Periodic A . Algorithm II. , v S Several questions arise concerning the estimation accuracy of algorithm I. The first question concerns the effect of measurement noise on the optimal estimate. It is seen from (4.7) that h, since it approxi-mates g, may lose its significance if the data representing g is noisy. To some extent this loss of significance can be taken into account by decreasing the absolute value of the elements of the weighting matrix Q2 (see (4.8)). A second question arises concerning the satisfaction of the controllability condition. Algorithm I must be modified if practical answers to the above questions are to be found. When no information is available about noise statistics there is no unique best method for for mulating an estimation problem, since the choice of a cost index, such as the quadratic cost index given by (4.8), is an arbitrary one. Con sequently, estimation algorithms can only be judged on the basis of trade-off between computational complexity and estimation accuracy. With the above discussion in mind, an algorithm is developed for estimating time-varying periodic parameters. Let x3 = z - x2 (4.22) x2 =-Afx-2 ••+ e(-x^);a., x-2(tQ) = 0* * (4.2-3) It is seen from (4.5) and (4.20) that X3 = AfX3 + 9(x2 + X3)-' X3(to) = 0 (4.24) Equation (4.23) is taken to represent the dynamical constraint for the new problem formulation. It is seen from (4.1) and (4.23) that the dynamics are noise-free. A cost index of the form given by (4.8) must be specified to complete the problem formulation which is based on the following assumptions: (1). L(x^) is well-conditioned (see (4.13)). (2) an estimate I of a is available and an estimate x^ of x^ can be found by solving (4.24) using successive approximations. That is, x^ is defined x3 = Afx3 + Q(x2)a, x3(t0) = 0 • x0 = AjL + 0(x.)a, x_(t ) = 0 (4.25) Z J. I. 1 — I o (3) The cost index is (see (4.8) and (4.22)) 3 - X2H2q1 + MH " X3 " X2H\]DT (4.26) o Successive approximations converge if the norm of a is sufficiently small. This restriction is imposed to achieve convergence in the two-stage estimation of A., and A . (If A., is known, then the first estimation f v f stage is unnecessary; (2.4) could then be used to define x^ and no norm constraints need be imposed on a). It is always possible to avoid ill-conditioning of L(Xp) by modifying the values of its time-varying elements. The modifications, however, cannot be arbitrary. They must enter into the problem formulations in a mathematically consistent manner. In the present problem formulation it is seen from (4.22) and (4.26) that L(x^) is replaced by- L(Xj-).. Equations- (4. 23) , (4..25). and. (4..26). define-a--two-point boundary-value problem for the optimum estimate. be a fixed estimation interval. At the end t^ of an estimation interval a new measurement sequence is taken over an observation interval t, < t < tf + HT. That is, the next estimation interval is t + £T < t < Let T be the length of a sampling interval and let t,. - t = NT o (N + £)T (see Fig. (4.1)) . •OBSERVATION INTERVAL FOR NEW MEASUREMENTS 0IT NT(UN)T h Ha =H ESTIMATION INTERVAL OPERATION PERIOD flG.(4.l) THE RECURSIVE OBSERVATION - ESTIMATION INTERVALS. In order to solve Eqn. (4.25) over this interval it is necessary to extrapolate a^(t) . This can be done in a variety of ways. One method is to fit a polynomial to sampled values ji(to + kT) , k = 0, N, and to use the polynomial approximation for extrapolation. ..Alternatively, , if it is known that the length of the estimation interval NT is greater than the smallest period of the periodic components of a_, then the known periodicity of a_ can be used for extrapolation. The following is a summary of Algorithm II: 1. A measurement sequence g(tQ + kT), h(tQ + kT) and an estimate a_(tQ + kT) over an estimation interval k = 0, N are assumed known. Take a new measurement sequence over an observation interval k = N, N + I and define a_ over this interval by extrapolation. 2. Solve the two-point boundary-value problem associated with the new estimation interval for a new estimate (Equations (4.23)*, (4.25) and (4.26)). The two steps of the algorithm are repeated recursively. 4.5 Estimation of a Periodic A^. Algorithm III. Algorithm II makes recursive use of measurement data and pre-vaious estimates. It is based on the solution of an optimization problem having.-the noise free, dynamics, of. (4.. 23). and a cost index given by (4.26) where x^ is generated in a smoothed form by integrating (4.25). Conse quently, it can be anticipated that algorithm II will be superior to algorithm I if measurement noise is present. There is no guarantee, however, that the controllability condition is satisfied. That is, L(x^) may be ill-conditioned. As mentioned in the previous section, one method for overcoming ill-conditioning is to modify the time-varying elements of L by modifying the problem formulation. Let z^ be an estimate of z^ and let X2 = AfX2 + 9(xl + V-' x2(toJ = ° (4.27) represent the dynamics associated with a third problem formulation. The cost index is taken to have the quadratic form (see (4.8)). J = \ t/tf [||g - *2\\\ + ||h - x2||2Q2]dT (4.28) o Comparing (4.27), (4.28) with (4.5), (4.8) it is seen that z^ should be chosen so that x2 approximates z^. (Ideally, in the noise-free case and when L(x ) is well-conditioned, then the choice z, = z, can be made. P 1 1 This results in x2 = z^ and Algorithm III degenerates into Algorithm I). If L(x^) is ill-conditioned Algorithm II will fail. This ill-conditioning can be overcome in Algorithm III provided z^ is chosen so that L(x^ + z^) is well-conditioned. A further constraint is imposed on z^ by the require ment that the two-stage estimation procedure converges. This requirement, which was discussed in the previous section, led to the use of successive approximations and the assumption that the norm of a_ must be sufficiently small. It is evident that z^ is constrained by several conflicting re quirements. Some sort. of. mathematically tractable, approach must be taken, which accounts for the various trade-offs associated with ill-conditioning, dynamical approximations and two-stage convergence. Solving (4.5) by successive approximations using z^ = 0 as an initializing function and an estimate a. gives x2 - Afx2 + 6(XL)I, x2(tQ) = 0 x„ = A.x„ + 9(x. + x0)a, x0(t ) = 0 (4.29) J 1 J 1 Z — JO 6; for the next two approximating functions. Let z^ be defined by z1 = ax2 + 3x3 (4.30) The general structure given by (4.30) allows the trade-offs mentioned above to be investigated as a function of two scalar parameters a and $. The choice a = 0, B = 1 results in z^ being defined by successive approxi mations given by (4.29). Numerous algorithms for the solution of linear algebraic equations are given in [24]. Abramov's method for accelerated convergence results when a = 2, 3 = -1. (The justification for using these algorithms is easily seen by the discretization and subsequent rearrangement of (4.5) into a linear algebraic equation). Algorithm III is based on the solution of the two-point boundary-value problem defined by (4.27), (4.28), (4.29) and (4.30). The scalar parameters a. and p can. be. determined by numerical. experimentation., or, they can be chosen on the basis of successive approximations (a = 0, 3 ='. 1) or Abramov's method (a = 2, 3 = -1). 4.6 TwO-Stage Estimation of A^ and A^. During the first stage an estimate of A^ is found given an estimate of A^. Three different algorithms have been developed in the previous sections for the second stage during which an estimate of A^ is found given an estimate of A^. The two stages are performed sequentially. Together they comprise a composite estimator for A. For reasons of simplicity, the discussion is limited to the development of a first stage algorithm which is associated with second stage Algorithm I. The modifi cations required if second stage Algorithm II or III is used are straight forward. Equation (4.4) defines the dynamics associated with the first stage of estimation. If x^ were available through measurements, (4.4) could be discretized and estimated by use of standard minimum mean-square error estimation procedures. However, x^ is not directly available. Furthermore, a two-stage estimation is proposed and the interaction between stages must be accounted for if a convergent sequence of esti mates is to be obtained. Chapter 2 has suggested augmenting (4.4) by additional states and controls to reduce the effect of quantization noise which is intro duced by discretization. It is shown that the use of augmented states and controls results in rapid convergence of the estimates. The augmented states and controls are defined by successive integrations of (4.4): y. = Ay. + B u. Jx fx i yi(to) =i 1 I 0 , i > 1 (4.31) A ft A yi = o yi-ldT> yl = Xl u. = •/"t u. ^dx, u, = u (4.32) l o i-l 1 For' composite estimation" of" A, successive- integrations has' the beneficial effect of reducing the influence of the time-varying components of A^ on the augmented states. The augmented system of equations (4.31) are dis cretized and rearranged into a set of linear algebraic equations. The optimum estimate for A^ if taken to be the solution with minimum mean-square error (For details see Chapter 2). In the case of composite estimation the above procedure must be slightly modified since x^ is not directly available. Consequently, x^ must be estimated from existing data. From (4.3) and (4.5) it is seen that XJL = xp - (4.33) where z- = Acz, + 9(x )a, z.(t ) = 0 (4.34) 1 fl p —' 1 o In (4.34), and a. are estimates of A^ and a_ respectively. The first stage of estimation consists of a minimum mean-square error solution of the discretized equations (4.31) using (4.33) and (4.34) to define x^. The new estimate of A^ is used in place of the old estimate provided that it results in improved prediction of the state of the sys tem. The following cost index can be used as a basis for making this decision: Jl = I t/f'l*l'xp " Xl " SiM2Q3d"T (4'35) In (4.35), is a positive definite weighting matrix. The rate of-change of between estimation intervals is useful in making decisions concerning the displacement length (£T) of the observation interval. The displacement can be increased when the rate of decrease of goes below a preset threshold. 4.7 Examples In all examples the following computational procedures were adopted. The two-point boundary-value problem was solved using a fourth-order Runge-Kutta method with an integration step-size of 0.025. The time-varying estimates a.(t) were generated in the computer by storing the values at grid points and taking a piecewise constant approximation 70 a(t)=a(tQ+kT), kT < t<(k + 1)T). In the case of algorithms II and III, extrapolation of _a was performed using a very simple constant approxi mations, _a(t^ + kT) = a_(t^) , k = 0, , and the initializing functions for successive approximations were taken to be identically zero. 4.7.1 Example 1. Consider a first order system x = -(1 + b sino)t)x + u (4.36) where the control u is taken to be a step function of magnitude equal to 3 and x(0) = 0 with T = 0.025, I = 4, N = 80 (4.37) The length of the whole operation period for all the examples is taken to be 10 seconds. The length of the estimation interval, as seen from (4.37) is NT = 2 seconds and the same has been considered for all the subsequent examples. The maximum percentage error of the estimate of a(t) = b sin wt for different b and w values are given in Table (4.1). Table (4.2) shows the results when a(t) is a square wave: f 0.5,. 0 < L< 0..5. a(t) =\ L-0.5, 0.5 4 t < 1.0 (4.38) In the case of algorithms II and III, two estimation intervals were re quired to obtain the results given in the tables. Algorithm III was tried with successive approximations (a = 0, 6=1) and Abramov's method (a = 2, g = -1). 4.7.2 Example 2. Consider the second order system x^ = -(1 + b^sin io^t)x2 + u^ x2 = -(2 + t>2 sin o)2t)x1 - x2 + u2 (4.39) where is taken to be a step function of magnitude equal to 2 and {t , 0 < t < 1 , t > 0  <_ 3 3 (4.40) In addition to (4.37) the following values are used: x^(0) = 0 = x2(0), b^ = 0.4, b2 = 0.6, co^ = 3, = 2. Table (4.3) gives the maximum per centage estimation errors for the parameters a^(t) = b^ sin co^t and a2(t) = b2 sin o)2t, respectively. In the case of algorithms II and III, two estimation intervals were used. 4.7.3 Example 3. Consider the third order system x^ = (1 + b^ sin to^t) x2 + u^ x2 = -(2 + b2 sin ui^t) x3 - x2 4- u2 x3 = -(3 + a3(t))x1 - 2x3 + u3 (4.41) where the controls u^ and u3 are taken to be step functions of magnitudes 0.1 and 0.5, respectively and where ft , t < "2 =ll.«> 3 3 (4.42) The parameter ^(t) is periodic and given by a3(t) -1 , 0 <t < 0.25 Ul , 0.25 < t <_0.5 (4.43) In addition to (4.37) the following values are used: x^(0) = ^(0) = x3(0) = 0, b± = 0.3, b2 = 0.5, o>1 = 3., OJ2 = 2. Table (4.4) gives the maximum percentage estimation errors. In the case of algorithms II and III , four estimation intervals were used. 4.7.4 Example 4. The next two- examples illustrate the composite two-stage estimator. In both examples the first-stage of estimation was initialized by taking a. = 0 as the initial estimate for ci and only one set of aug mented states and controls was used (i = 1, 2, see (4.31) and (4.32)). Consider the second order system x^ = (1 + b^ sin to^t)x2 + u^ x2 = (b-2- sin o)2t.)x1- + a2x2 + U2 (4.-44)-where the controls u^ and u2 are taken to be step functions of magnitude equal to 2. and 3., respectively. In (4.37), £ is taken equal to two and the following values are used: x^(0) = = ^» ^1 = 0.3, t>2 = -0.5, co^ = 1, co2 = 2, a3 = 1. After two estimation intervals the maximum per centage estimation error for a^ was 1.2% and after ten estimation in tervals, the maximum percentage estimation errors for a2 and a3 were 1.8% 72 and 1.0%, respectively. 4.7.5 Example 5. Consider the second order system x^ = (1 + sin ui^t)x^ + X2 = a2Xl + a3X2 + U2 (4-45) where b^ = 0.4, a^ = -2, a^ = -1, co^ = 2 and where all remaining data is the same as in Example 4. After the second estimation interval, the maximum percentage estimation errors for a^, a^, and a^ were 0.9%, 0.4% and 0.3%, respectively. 4.7.6 Example 6. In the previous examples the measurements were noise-free. Example 1 has been tried with zero mean gaussian measurement noise with unity standard deviation and the random noise generator has been randomly initialized. A peak r.m.s. signal to noise ratio of 10 was taken. The following values were chosen: b = 0.4, co = 1, a = 2, (3 = -1. The maximum percentage estimation errors for algorithms. I.,.. II and. III. were 6%., 3.2% and 2%, respectively. In the case of algorithms II and III, seven and four estimation intervals were used, respectively. The examples show that the proposed algorithms and the two-stage estimation method can give accurate estimates of linear system parameters. If the measurements are noise free, Algorithm I gives the most accurate estimates. In the case of measurement noise, Algorithms II and III give better estimates. In the above examples, no problem 11 was encountered with ill-conditioning in any of the above cases. 75 b 1 W 0.5 1.0 1 / 0.6 % 0.75% JT 3-0 10 % 15 % 2.0 0.7% 0.85% 777 °(=2, 3=-1 3.0 0.9 % 1.0% 777 3.0 0.93% 0.95% TABLE (4.1) TABLE (4.2) TABLE (4.3) TABLE (4.4) 5. OPTIMUM CONSTANT FEEDBACK GAINS FOR LINEAR SYSTEMS AND THE PROBLEM OF TRAJECTORY SENSITIVITY REDUCTION 5.1 Introduction The optimum solution of the linear state regulator problem with quadratic cost has assumed a central position in the development of modern control theory. The optimum gain, however, is in general time-varying and this has severely restricted practical applications. Practical considerations generally specify constant or piece-wise con stant feedback gains. Considerable effort has gone into the development of design procedures for constant feedback gains and in the design of specific optimal controllers [25-29]. There are, however, computational as well as practical diffi culties associated with optimum constant gains. Specifying a constant gain results in a constrained optimization problem and the solution of the associated two point boundary value problem is generally more diffi cult. Furthermore, the optimum constant gain depends on the initial conditions. One method that has been proposed to overcome these diffi culties is to introduce the trace of a cost matrix which then allows ini tial conditions to be eliminated in the problem formulation by arbitrarily assuming a uniform distribution of initial states over a sphere [25-27]. From a design standpoint, however, it is important to investi gate the effect of initial conditions on optimum feedback gains before computing average suboptimum gains. Furthermore, in many practical situ ations, state trajectories are controlled, and the effect of initial conitions on optimum gains can become important. This occurs, for example, in certain schemes for combined estimation and control. During a finite observation interval, measurements of the state are taken and used to determine parameter estimates of a linear model. The linear model is used to compute a control strategy which is applied in a subsequent control interval [11], [5]. This scheme seems particularly attractive for adaptive systems. Its practical implementation, however, requires an on-line computational algorithm for determining optimal or suboptimal gains. Trajectory sensitivity minimization is an important practical problem that has received considerable attention [30-37], A widely adopted method is based on augmenting the state vector with the sensi tivity vector and formulating a linear state regulator problem for the augmented state. However, if time-varying gains are permitted, the problem, is ill-posed [35, 38]. A well-posed, traj.ectory. sensitivity mini mization problem requires that the gain be constrained to be time-invariant. In general, the sensitivity vector is affected by initial conditions. Consequently, if feedback is to minimize trajectory sensi tivity, the effect of initial conditions on the optimum constant gain must be evaluated. It may be possible to choose a suboptimal average gain. A unified treatment of several problem formulations dealing with time-varying gains and constant gains is presented. Extensive" use is made of the known results for the linear state regulator problem to develop a simple algorithm which results in rapid convergence. An approximate method is developed to replace the nonlinear two point boundary value problem by an algebraic optimization problem. This method appears to be suitable for on-line computation of suboptimum constant feedback gains. 5.2 Problem Formulation Consider a linear time-varying system x = Ax + Bu, x(0) = Xq (5.1) where the state x and the control u are n x 1 and m x 1 vectors, res pectively. The optimal control which minimizes the quadratic cost index (prime denotes transposition) 1 tf J = ± / (x'Q-x + u'Ru)dt, (5.2) o 2 o 1 where is positive semidefinite and R positive definite, is given in feedback form by u = -R •''B'K.jX. (5.3) The time-varying gain matrix is found by solving a matrix Riccati differential equation and has the property that = K|. To ease the notational burden, argument values of functions are omitted when they appear unessential and only one parameter q in the matrix A is considered to have a significant effect on sensitivity. Let z = represent the dq closed-loop trajectory sensitivity vector. A feedback structure of the form u = -R-1B'(K x + K2z) (5.4) is postulated, where = K| and = K^. This specific structure is postulated because it resembles the feedback structure of the optimal control resulting from solving the state linear regulator problem (see (5.3)). Substituting (5.4) into (5.1) and neglecting the second s' order sensitivity function —j yields 9q z = 1^ x + (A - FKjz, z(0) = 0 3q J. (5.5) where F = BR "*"B'. By defining an augmented vector y' = (x', z'), (5.1) and (5.5) can be combined into a composite vector differential equation y = (A - FK)y = Cy, y(0) = yr where (5.6) A = 3q 1J F 0 0 0 - A , K = Kl K2 LK2 K3 ' yo = x (5.7) Equation (5.6) represents (5.1) and (5.5) when = K^. This constraint is relaxed in (5.6) in order to obtain a unified treatment of the linear and nonlinear formulations. The design objective is to choose gain matrices and so that (5.2) is minimized subject to the additional constraint that the feedback (5.4) reduces, in some sense, the closed-loop trajectory sensitivity vector z. A mathematically tractable approach that has been widely adopted in sensitivity studies is to neglect second-order sensitivity terms so that (5.6) is valid and to introduce a cost index J(Kr K, tf) = j / f y'(Q + KFK)y dt (5.8) where Q is an augmented positive semidefinite weighting matrix [32-38], The optimum gain matrix minimizes (5.8) subject to the differential constraint (5.6). It will be seen that the sensitivity problem formu lated by (5.6) and (5.8) is well-posed only if the gain matrices are constant. The argument list in the cost index (5.8) is used to distinguish between different problem formulations. Let I(K^, K, t^) and I(K^, K, t^), where is arbitrary, represent the linear and nonlinear formulations, respectively (If (5.6) is rewritten in the form y = Ay + Fu, then for the linear formulation the matrices A and F are completely known and inde pendent of K^. This results in linear dynamic constraints. For the nonlinear formulation the parametric matrix A depends on the unknown gain matrix K^. This results in nonlinear dynamic constraints). The Hamiltonian for either formulation is H = p'Cy - |- y'(Q + KFK)y (5.9) where the costate vector p is defined by p = -C»p + (Q + KFK)y, p(tf) = pf = 0 (5.10) Gradient matrices (see Appendix III) are used to derive the gradient condition for a minimum. In order to use gradient matrices all matrix elements must be independent and cannot be constrained by the symmetry condition K = K'. This is easily accomplished by taking K = K + K', where K is a matrix with arbitrary elements. In order to introduce a compact notation which proves useful in subsequent discussions, it is assumed that F is time-invariant. The gradient condition for the linear formulation I(K^, K, t^) can then be expressed in the form (see Appendix V) 8J(K K,t ) t i——= /J-_-5-dt=$ KF + FK$ +F$ +<S> F = 0 3K o 3K yy yy py yp (5.11) where $ k yv» dt (5.12) yv o J 81 (Equation (5.11) remains valid in the general case of time-varying F provided that KF and F are moved under the respective integral signs. Keeping this in mind, it is seen that the form (5.11) is actually not a restriction but rather a notational convenience). The gradient condition for the nonlinear formulation I(K^, K, t^) differs from (5.11) in an additional term which arises from the dependency of A on K^. Let p' = (p|, pp be a partitioning of the costate vector which corresponds to the partitioning of y and let K^ = L + L', where L is arbitrary. The additional term which is associated with K^ is a tf ~ f Tr(-FKlZp!)dt = F$ + $ F 3L o 1 r2 p z zp„ (5.13) Introducing the augmented matrix E = F$ + $ F ' 0 P2z Zp2 ! 0 1 0 (5.14) allows the gradient condition for the nonlinear formulation to be ex pressed in the form 8K " = E + $ KF + FKo + F'f + $ F = 0 yy yy py yp (5.15) Minimizing J(K^, K, t^) is a constrained optimization problem. Conse quently J(K*, K*, tf) = J(K1} K, tf) (5.16) where K* is the optimum gain matrix for the nonlinear formulation. If K^ = K*, it is seen from (5.16) that the minimum for the linear formulation 82 is achieved when K = K*. The condition E* = 0, which is satisfied when F$ + $ F = 0 (5.17) p2z zp2 is an optimality condition. It is interesting to note that if time-varying gains are permitted, the optimality condition (5.17) becomes a singularity condition which is impossible to satisfy [32, 37]. This indicates that the sensitivity problem formulated with time-varying gains is ill-posed and cannot, in general, meet the design objectives. Equations (5.11) and (5.15) can be simplified and unified by introducing a matrix S defined by FS$ + $ SF = -F$ - $ F - E (5.18) yy yy py yp (Several methods are available for solving an equation of the form (5.18) for S. See Appendix II). Substituting (5.18) into (5.15) yields 9J(K K,t ) \v = $ (K - S)F + F(K - S)$ =0 (5.19) 9K yy yy Equations (5.6), (5.10) and (5.19) define a nonlinear two-point boundary value problem which must be solved in order to determine the optimum gain matrix. 5.3 Unified Treatment of Different Formulations Because it is a constrained optimization problem, the computation of optimum constant gain matrices is generally more difficult than the computation of unconstrained time-varying gains. It is only in the special case of time-invariant systems and infinite tf that the problems are computationally equivalent. From an engineering point of view, however, constant gain matrices are easier to implement and therefore preferable, despite the possible computational difficulties. It is evidently desirable to retain the comparative computational simplicity of the unconstrained time-varying gain case, as far as this is possible. This requires developing a unified treatment for the various formulations and then making maximum use of the known results for the linear state regulator problem with time-varying gain and quadratic cost. Substituting the Riccati transformation p = -Sy (5.20) into (5.10) and using (5.6) yields ' S* + SC + C'S = -Q - KFK, S = 0 (5.21) Consider the general case of a time-varying system with time-varying gains. In the linear formulation I(K^, K, t^), E = 0. Substituting (5.20) into (5.18) yields tf / (F(S - S)yy' + yy'(S - S)F)dt = 0 (5.22) Equation (5.22) is satisfied when S = S and (5.19) is satisfied when K - S =• S. Substituting this- constraint into (5.21)- gives- the- standard matrix Riccati differential equation for the optimum time-varying gain. In the case of a constant gain matrix the above procedure must be modified. Heuristic considerations indicate that some sort of time average of S should result in a good suboptimum constant gain (For a time-invariant system, the steady-state solution of (5.21), taking K = S, often gives a good suboptimum result). The correct time average to use is given by (5.22). Replacing Sy by p, (5.22) can be solved for S by use of Kronecker products (See Appendix III). The vector representation (S) of the matrix S satisfies the linear equation (5.23) A method of successive approximations could be considered for solving the two-point boundary value problem. Using a nominal gain K = K^, (5.6) is integrated in the forward direction and (5.16) in the backward direction. Equation (5.23) can then be solved for S and (5.19) used to determine the gain increment. The updated gain is and the procedure is repeated iteratively. While successive approximations' seems straightforward, they generally fail to converge in the case of nonlinear two-point boundary value problems. The special structure of the formulated problem, however, allows some insight to be obtained concerning the convergence of successive approximations. 5.4 Gradient Descent and Successive Approximations 3 J Let (6K) and (—) be the vector representations for a matrix v 3K v gain increment and the gradient matrix, respectively. The incremental cost associated with 6K is given by (5.24) (5.25) where (5.26) is positive semidefinite (see Appendix III). Consequently, if 6K = w(S - K), 0 < w = 1 (5.27) then 5J = -w(S - K)'W(S - K) = 0 (5.28) The gain is updated by \+l " \ + W(\ " V (5.29) where the subscript k is used to indicate an iteration stage. For the problem formulation under discussion, successive approximations (5.24) is a special case (w = 1) of. gradient descent (5.29) with a weighting matrix (5.26). Consequently, the proposed successive approximations method is equivalent to a deflected gradient method. In the case of quadratic cost indices it is.known that some deflected gradient methods (such as conjugate gradients) result in significantly improved convergence rates. Having established the property that (5.24) is a deflected gradient descent, it is of interest to see whether it has the properties of improved convergence, such as quadratic convergence. It is known that successive approximations results in quadratic convergence in the case of unconstrained time-varying gains. Quadratic convergence can be proven by considering (5.21) for two successive iterations. Introducing k+1 S, , the linear matrix differential equation for AS, can be solved using standard procedures, giving -ASk(0) o /fcf K+±(T, 0)^ - - V^Ct, 0)dx + Cf *k+l(T>°)C(\+l ~ ~ W + (\ " \+l)F"(\+l " Sk)]W(T'°)dT (5.30) ' In (5.30) ^(t, x) is the state-transition matrix of (5.6) at the k-th iteration stage. It is easily shown that the cost index at stage k is given by J, = — y'S.y . The incremental cost AJ, = J, ,, - J, is found k 2 o k o  k+1 k from (5.30) and can be expressed in the form -AJk - TrtlyyCkKK^ - \)H\+± ~ VI + TrK^)^ - K^) ] (5.31) where - Cf [Wk+i(\+i - VF + F(\+i - VWk+i^ (5-32) and where the notation $ (k) is used to indicate that (5.12) is evaluated yy at stage k. The gradient condition (5.19) can be obtained from (5.31) by noting that Tr[(A + A')B] = 0 for arbitrary symmetric matrix B only when A + A' =0. Equating (5.32) to zero yields (5.19). The remarkable properties of the case with unconstrained time-varying gain are seen from (5.32). By setting Kk+^ = Sk> (5.32) vanishes and the decrease in cost (5.31) is quadratic [38]. This is accomplished by solving (5.21) .using successive approximations. The same quadratic convergence would result in the case of constant gain if could be chosen to make (5.32) vanish. In principle this may be possible. However, y^-^ depends on K^-^* Consequently, any attempt to satisfy the condition (^)^ = 0 would result in poor convergence. It is reasonable to look for an approximation by replacing 87 y^+^ in (5.32) by y^. Let (^r)^ represent the matrix arising from this substitution. It is seen that <M>k • Vk)(V-l " VF + f(\+l " V$yy(k> (5i33) where S^ is defined by (5.22) evaluated at stage k. It follows from (5.33) that (^O^ = 0 if is determined by successive approximations (set w = 1 in (5.29)). Since (7^). is an approximation to (7^), it is reasonable to anticipate that successive approximations results in rapid convergence. A rigorous proof of quadratic convergence would require shoxtfing that the second term of (5.31) is either never negative or, if it is negative for a certain gain sequence, that its absolute value is less than a fixed fraction of the first term, which for Kjc+-^ f1 is always positive. From a practical point of view such a proof is not essential. If the choice w = 1 should happen to result in divergence, then a choice w < 1 can be made which will always give a convergent sequence (see (5.28)). The proposed computational algorithm is the following. Choose a nominal K = K^. Integrate (5.6) in the forward direction and (5.10) in the backward direction. Solve (5.18) for S^ and update the gain by use of (5.29). The rate of convergence is controlled by w. Computationally, the proposed method is similar to the conventional steepest descent procedure. The. essential, difference is in the use. of a variable weighting matrix (5.26) to give a deflected gradient. The associated change in the gain matrix is actually given by the simple successive approximations formula (5.29). Limited computational experience to date indicates that this variable gradient deflection allows a large step (w = 1) to be taken resulting in rapid convergence. Note that the linear and nonlinear formulations of the sensi tivity problem are treated in a unified manner. This arises from the fact that in both formulations a nominal gain K is chosen to perform the computations. The above algorithm is suitable for off-line design studies. Since the optimum constant gain depends on initial conditions, there is a problem of practical on-line implementation. A method for eliminating initial conditions and deriving suboptimal feedback control laws is given in.[39]. A method which is suitable for combined on-line esti mation and control is discussed in Section (5.6). 5.5 Optimum Constant Gain: Trajectory Sensitivity Reduction Because of its engineering importance, open-loop and closed-loop system sensitivities have been extensively investigated (see Reference [32] for a bibliography). However, for the closed-loop sensitivity formulation given by (5.6) and (5.8), a paradox arises if unconstrained time-varying gains are permitted. This paradox has been discussed in the literature but its effect on sensitivity reduction does not seem to have been fully appreciated [33, 37]. A simple example can be- used to illustrate- the paradox. Consider-x = -ax + u, x(0) = x (5.34) o z = -x - az - aKjZ, z(0) = 0 (5.35) u = -Hi K - K2z = 0 (5.36J = -| /°°(x2 + Q2z2 + Ru2)dt (5.37) and suppose that and R are both chosen large compared to unity. This would seem a reasonable choice to make in order to keep both sensitivity and control effort small. The optimum time-varying gains, however, satisfy the constraint (5.36) and are infinite. This can be seen by taking K -> + » in (5.35). Equations (5.34)-(5.37) then yield. 1 aK x = 0 (5.38) K2 = K2 (5.39) (5.40) The optimum (infinite) gain can result in a very significant smaller cost than that associated with other suboptimal (but physically realizable) gains (see [37]). It is evident, however, that minimization of (5.37) has not succeeded in the design objective, which is to decrease system sensitivity. The system (5.34) with the optimum control (5.36) operates open-loop. The same difficulty is experienced with a general n-th order system (see (5.1), (5.4) and (5.5)). The sensitivity problem defined by (5.6) and (5.8) is consequently ill-posed if unconstrained time-varying gains are permitted. From both a computational and an engineering point of view the simplest constraint is that of time-invariant gains. It is interesting to note that this choice was made by Kreindler in his design studies which dealt exclusively with linear time-invariant systems and t^ = + <» [3]. It seems natural to choose constant gains in this case. However, if t^ is finite and/or the system is time-varying, then the choice of constant gains is not obvious. 90 5.6 On-Line Computation of Gains From an engineering point of view a very attractive feasible on-line adaptive control scheme is to identify a process as a linear time-invariant system of the form (5.1) over a finite observation interval and then apply a control of the form (5.3) over a subsequent finite control interval. In such a combined estimation and control scheme it is often important to reduce the system sensitivity to incorrect parameter estimates. It is, however, not generally computationally feasible to solve a nonlinear two-point boundary value problem on-line. Some form of suboptimal control must be introduced. It is known that the steady-state solution of the matrix Riccati equation (take K = S in (5.21)) often results in a good suboptimal control. With this in mind let P = -My - Ny + e (5.41) where M and N are constant matrices. Substituting (5.41) into (5.10) and choosing M to be defined by MC + CM = -Q -KFK, (5.42) which is the steady - state form of (5.21) gives e + Ce = -(NC + CN)y, ef = (M + N)yf (5.43) A For reasons of notational convenience the linear formulation I(K^, K, t^) is discussed. Substituting (5.41) into (5.11) yields F(S - M - N)$ + $ (S - M - N)F = -F$ - $ F (5.44) yy yy ey ye The decomposition (5.41) does not uniquely define N and e. Consider the possibility of defining N in terms of observational data so that the right hand side of (5.44) has a negligible effect on S. The updated gain, as given by (5.24), is then S = K2 = M + N (5.45) It is seen from (5.43) that e is "small" when N and e^ are "small". Besides the obvious trade-offs involved in keeping N and e^ simultaneously "small", there is also the requirement for on-line computation using observational data. The following cost index accounts for the various factors which enter into an "optimal" decomposition (5.41): Je = ^ / f (e + C'e)*W. (e + C'e)dt + ^ e/W.e. (5.46) 2 o 1 2 r 2 f In (5.46) W^ and are positive definite weighting matrices. Sustituting (5.43) into (5.46) it is seen that Je is an algebraic function of N, $ and y^y1^' Setting N = L + L', where L is arbitrary, allows gradient matrices to be used to determine the minimum. This yields 1^- = TC1 + CT + W0(M + N)i|) + iK(M + N)W_ = 0 (5.47) oL 2 t r 2 where T = W1S$ + $ SW, 1 yy yy 1 S = NC + C'N <J>f = yfy^ (5.48) Equation (5.47) can be solved for N by use of Kronecker products. It is seen that (S) = G. (N) v 1 v (T)v = G2(S)v where (TC* + CT) = G.(T) v 3 v [G1G2G3 + G4](N)v = ~G4(M)v (5*49) G1 = C,©I + I©C' =G3 2 yy 1 1 yy G4 = <|»f ® W2 + W2 O ^f (5.50) Efficient algorithms are available for evaluating (5.50) [40]. The solution of (5.49) for N involves one matrix inversion. The system parameters... in. the. C. matrix- are- determined, by the identifier from, ob.-servational data. The matrices $ and ii, can be found by a forward yy f integration of (5.6). Alternatively, it may be possible to use measured states or state estimates to evaluate these matrices. The gain compu tation (5.45) is then an algebraic problem and replaces the two-point boundary value problem. On-line use of (5.42), (5.49) and (5.45) to update the gain is computationally equivalent to solving these equations by successive approximations. The first iteration must result in an improved sub-optimum gain. It is desirable, but no essential, that successive itera tions result in a convergent suboptimum sequence. Since and W2 are arbitrary positive definite weighting matrices, it is not possible to reach general conclusions concerning the convergence of successive approximations. It is easy to show, how-ever, that given a nominal K, weighting matrices exist which result in one-step convergence. Suppose that the optimum gain K* is determined by off-line computations. Equation (5.42) can be solved for M. Let (see (5.45)) N = K* - M and take W^ = I. Equation (5.47) is a Liapunov type equation which can be solved for a positive definite l^. With the weighting matrices so definedo there is one-step convergence to the optimum gain. It should therefore be possible to choose weighting, matrices so that the first iteration results in an improved suboptimum gain for small variations about nominal parameter values. If the system is time-invariant and if N is small in the sense that terms involving products of the elements of N can be neglected, then it is seen by sub stituting (5.45) into (5.42) that M is the steady-state gain. In this case (5.42) can be solved for the steady-state gain and N evaluated by (5.44). 5.7 Examples To give a non-trivial illustration of the successive approxi mations as well as the approximation method discussed in Section (5.6) t^ should be chosen smaller than the settling time of the open-loop response for the nominal gain, and the nominal gain should differ signi ficantly from the optimum gain. These requirements are satisfied in the case of the second order system Cl = "qX2 + U ' xl(°) = 1,0 (5.51) x2 = xx , x2(0) = -1.0 t J = \ f f (x2 + 3x2 + u2)dt o I o 1 2. 94 by taking q = 1, tf = 2.0, and = (-10, 10). (The system is an oscillator for u = 0 and is unstable for the chosen nominal gain). Successive approximations ((5.6), (5.10), (5.23) and (5.24)) required eight iterations to converge. No attempt was made to determine weighting matrices which would result in rapid convergence for the on-line approximation method ((5.42), (5.49) and (5.45)). An arbitrary choice = W2 = I was made. Six iterations were required. The results are the following: 1. Successive approximations K^** = (-1.06, - .215), J ** = 1.518 2. On-line approximation method. K ** = (-1.713, - .952), Jq** = 1.565 (5.52) The minimum cost for unconstrained time-varying gain is J* = 1.503. It is seen that both methods give very good suboptimal constant gains. The increases are only 1% and 4% for case (1) and (2), respectively. To illustrate the sensitivity problem consider (5.51) to be augmented by (5.5) and take J = 5 ftf (z2 + z2)dt + J (5.53) o 12 o A comparison is made with Lee's et al. method which uses a time-varying gain and the choice t^ = 10 is made. The following results are obtained (a) Lee's method [35] (augmented state feedback; t = t ) u = -2x± - 1.5x2 + .25z1 - .5z2, J* = 1.94 (5.54) (b) Constant gain (augmented state feedback) u = -6.8xi - 6.6x2 + 4.66Z;L + 20z2> J* = 1.70 (5.55) (c) Constant gain (Non-augmented state feedback) u = -1.812X - 1.22x2, J = 2.04 (5.56) Figures (5.1)-(5.3) illustrate the improvements that result by-use of constant gain and augmented state feedback. The magnitude of the control is smaller, the cost index is smaller and there is a significant reduction in trajectory sensitivity. That the design objective of tra jectory sensitivity reduction is achieved can be seen by comparing the sensitivity functions of case (b) and case (c) (case (c) is based on (5.51) taking tf = 10). The optimality condition (5.17) was checked and gave (-.0002, .0004). It should be noted that Lee's method uses a time-varying gain which is a suboptimum solution of an ill-posed sensitivity problem. (The gain values given in (5.54) are for t = t ). It is therefore not surprising that an optimum constant gain can give significantly better results. Note that increasing t^ from 2 to 10 and using sensitivity augmentation does not result in an appreciable increase in the cost index (compare (5.55) and (5.52). The convergence of successive approximation (5.29) was investi gated by choosing K± = (-1.5, -1.0), K2 = (-0.1, -0.5) (5.57) as the initial gain. This choice is significantly different from the optimum gain (see (5.55)). Only five iterations were required to converge to the optimum. The following example is used to illustrate the on-line approxi mation method when it is applied to a realistic system of moderate complexity. Figure (5.4) illustrates a block diagram for the feedback control of a power generator connected to an infinite bus [41]. The matrices in (5.1) are A = 0.0 1.0 0.0 0.0 -0.676 -.25 25.0 0.0 0.0 0.0 -2.0 2.0 0.0 -0.06 0.0 -2.0 , B = 0.0 0.0 0.0 2.0 (5.58) The components of the state vector represent the following variables; x^ = power angle deviation, x2 = angular frequency deviation, x^ = mechanical power deviation, x^ = governor position variation, u = governor position control. The initial state is taken to be x^ = (0., 3.0, 2.0, 0.0) and in (5.2) and (5.46) the choice Q± = W = W2 = I and R = 1 is made. The choice t^ = 1 sec. is made to accentuate convergence problems, in case they should exist. The minimum cost using the optimum time-varying gain S(t) (see (5.21)) is J* = 33.215. Successive approximations for the optimum con stant gain is initialized using K = S(0). The results for the cases (a) K = M, (b) K = M + N, are as follows. k = 1 2 3 4 (a) J = 45.404 39.167 36.636 36.121 (b) J = 38.035 36.267 35.601 34.682 (a) u* = -.66x± - 1.22x2 -8.44x3 -3.35x4 (b) u* = -.069x -1.47x2 -10.67x3 -5.47x4 (5.59) (5.60) The suboptimal choice K = M in case (a) amounts to taking the solution of the infinite time case for the constant gain. It is known that (5.42) converges quadratically to the optimum gain for the infinite time case. The suboptimum choice K = M + N in case (b) is based on the approximation technique discussed in Section (5.6). It is seen to have a quadratic-type convergence and results in a smaller cost index than for case (a). In fact, the first iteration results in a significant decrease. A one-step improvement of this kind could be of practical significance in an on-line application. a FIG.(5.1) EXAMPLE I CONTROL AS A FUNCTION OF TIME FOR: (a) LEE'S METHOD, (b) CONSTANT GAIN AND AUGMENTED STATES, (c) CONSTANT GAIN AND NON- AUGMENTED STATES. 0.24 0.161 ^ 0.081 to 0.0 co-0.08 ^~ o UJ £ -0.24 FIG. (5.2) EXAMPLE 1. TRAJECTORY SENSITIVITY FUNCTION 37 FOR-(a) LEE'S METHOD- (b) CONSTANT GAIN AND AUGMENTED STATES, (c) CONSTANT GAIN AND NON - AUGMENTED STATES. VO VO FIG. (5.3) EXAMPLE I TRAJECTORY SENSITIVITY FUNCTION Z2 F0R: (a) LEE'S METHOD, (b) CONSTANT GAIN AND AUGMENTED STATES, (c) CONSTANT GAIN AND NON- AUGMENTED STATES. SPEED -GOVERNOR TURBINE GENERATOR FIG. (5.4) EXAMPLE 2. BLOCK DIAGRAM OF A SYNCHRONOUS GENERATOR VOLTAGE CONTROL SYSTEM. 6. OPTIMAL ADAPTIVE CONTROL 6.1 Introduction In this chapter we will discuss one possible strategy for optimal adaptive control which is particularly attractive from an engineering point of view. The strategy is to identify system parameters at the end of an observation interval and then to use the parameters to derive an "optimal" control for a subsequent control interval. This chapter develops one version of this strategy which is based on the identifiers and optimum controls developed in Chapters 2, and 5. 6.2 Optimum Control During the First and Subsequent Optimization Intervals The particular optimization technique to be used during the first optimization interval is of importance since the initial information about the system is usually insufficient to derive a good control. If parameter estimates are poor then the optimum constant gains may not be sufficiently close to the correct values. The problem of convergence can become serious during the first optimization interval. Hence, it is re commended to use the successive approximations method with a step size constant w < 1 (see equation (5.27)). Figures (6.1-a, b, c, d) show the time intervals associated with a recursive observation-optimization strategy which would be suitable for the optimal adaptive systems. The system is considered to start operating at time t and to terminate its operation at time t^. In Figures (6.1-a, d) (T) is the length of the sampling interval. The length of each observation interval is (NT) seconds and the start of the i-th INITIALIZATION INTERVAL OPERATION PERIOD -sa-f«a R J L t0-NT t o FIG. (6.1-a) THE INITIALIZATION INTERVAL AND THE OPERATION PERIOD. rT rT NT t t0-NT t0 tj° t°+NT FIG. (6.1-b) THE SUCCESSIVE OBSERVATION INTERVALS. FlG.(6.1-c) THE SUCCESSIVE OPTIMIZATION INTERVALS (l = r). IT , Tc FIG. (6.1-d) THE SUCCESSIVE OPTIMIZATION INTERVALS (l = 2r). observation interval is denoted by the instant (t^). The displacement between two successive observation intervals is (rT) seconds (see Fig. (6.1-b)) where r is an integer number. The length of each optimization interval is (T ) seconds and the start of the i-th optimization interval is denoted by the instant (t^) (see Fig's (6.1-c, d)). The displacement between two successive optimization intervals is (&T) seconds. The para-, meter I is a multiple of r, i.e., I = mr. The value of m need not be the same for every two successive optimization intervals. The choice of the parameter (m) is based on the rate of variation of certain perfor mance indices to be subsequently discussed (see subsection(6.3.4)). During the initialization interval (tQ-NT < t < t ) special input test signals can be used to drive the system. Measurements of the input and output at successive sampling instants are used to determine initial estimates of the unknown parameters of the system. The choice of the initial conditions at the start of the initialization interval (t = t -NT) o is completely arbitrary. At the end of each observation interval new estimates of the plant parameters are computed. At the beginning of each optimization interval new estimates of the "optimum" constant gains are made using the latest estimates of the plant unknown parameters. These new gains are used for the coming control interval of IT seconds. 6.3 Examples In this section a linear, continuous system with time-invariant parameters is considered. Some of the system parameters are unknown and an optimal adaptive control is computed for the system. The recursive observation-optimization strategy described in the previous section is followed. The identification of the unknown parameters is performed using the digital identifier developed in Chapter 2. The number of augmented sets of states for Examples 1 and 2 are one and two, respectively. The test inputs used during the initialization interval are taken to be step functions with unity magnitude, in both examples. The initial conditions at the start of the initialization interval are taken to be identically zero in both Example 1 and 2. It should be noted that use of a digital identifier results in quantization noise. Consequently, even though the examples treat deter ministic cases, the system models are perturbed by quantization noise. This noise also affects the value of the optimum constant feedback gains which are computed using the identified parameters. It should be further noted that the optimum constant gains depend on the system state at the start of an optimization interval. The examples illustrate convergence of" the proposed" adaptive"strategy in the presence of quantization noise. 6.3.2 Example 1 The system considered in this example is represented by the dynamic model Xl=alx2+U x2 =- a^x^ (6-. 1) where the nominal values for a^ and a^ are -1.0 and 1.0, respectively. The initial conditions are x^CO) = ^(0) = 1.0. The essential parameters of the recursive observation-optimization strategy (see Figs. (6.1-a, d)) are taken to be: T = 0.2, t, = 10., t = 0., r = I = 2, N = 10 f o T = NT = 2.0 (6.2) It is required to compute the optimum constant feedback gains and such that the performance index 2+t? J. = \ f 1 (x2 + x2 + u2)dt (6.3) l L • C 1 L I is minimized using the new estimates of a^ and a^ which are obtained at the end of the i- th observation interval. The optimum constant gains are computed using the successive approximations method developed in Section (5.2) with w = 1. The optimum solution is assumed to be obtained when the inequality 'Ji " Ji_1| = 10"5 ' (6,4) is satisfied, where J. is the value of J. at the end of the k-th iteration. l l The time t^ (see Fig.'s (6.1-c, d)) denotes the beginning of the i-th, optimization interval. During the first optimization interval the initial guess for the feedback gains and is taken to be = Y.^ = -10. For the subsequent optimization intervals, the initial guess for the i-th optimization interval is taken to be the optimum values obtained for the (i-l)-th optimization interval. The stopping rule (6.4) was satisfied after six iterations for the first optimization interval. No more than two iterations were required to satisfy the stopping rule given by (6^4) for the subsequent optimization intervals. The performance index of the optimal adaptive system is taken to be 25 J* = E AJ* (6.5) • a i=l 1 whe re the optimum value of J^ (see (6.3)) is denoted by J* and AJ* is i i the value of J* for an integration interval t^ < t < t^ + 0.4. Since the optimum constant gains computed at the start of the i-th optimization interval are only used for a period of (JIT) seconds starting from the c instant t., then J* as given by (6.5) is the total operation cost. For X cL the sake of comparison the optimum value of the performance index J = \ f (x2. + xl + u2)dt (6.6) o I o 1 l is computed using (6.1) with a^ and a^ being equal to their respective nominal values and the control feedback gains are taken to be the optimal time-varying gains. The optimum value of Jq is found to be 2.046 while the value of J* (see (6.5)) is equal to 2.202, i.e., the increase in the performance index value does not exceed 8%. Figures (6.2) and (6.3) show that the identification errors do not exceed 5%. 6.3.3 Example 2 The system considered in this example is represented by the dynamical model x^ = a^c^ + 2u x2 = a2x3 + 2u x3 = a3x2 - u (6.7) where the nominal values for a^, a2, and a3 are -1.0, 2.0, and -2.0, respectively. The initial conditions are x^(0) = x2(0) = x^(0) = 1.0. The essential parameters of the recursive observation-optimization strategy are as given by (6.2). It is required to compute the optimum constant gainsK , K2 and K3 such that the performance index 108 2+t° J.=|J 1 (x2 + x2 + 4 + u2)dt (6.8) i is minimized during the i-th optimization interval. The optimum solution is assumed to be obtained when the inequality (6.4) is satisfied. In this example, using the successive approximations method with w = 1 (see equation (5.27)) during the first optimization interval results in a diverging sequence of iterations. Thus w is taken to be 0.1 and the initial guess for the feedback gains are taken to be = -1, = -2, and = -3. Twenty iterations are required to satisfy the stopping rule (see (6.4)). For the subsequent optimization intervals, the initial guess for the i-th optimization interval is taken to be the optimum values obtained for the (i-l)-th optimization interval. No more than two iterations, are required to satisfy. (6.4) for the subsequent optimi zation intervals and w is taken equal to unity. For the sake of comparison the optimum value of the performance index 1 10 2 2 2 2 J = i /1U (xf + x^ + x^ + u )dt (6.9) o 2 o 1 2 3 « is computed using (6.7) with a^, a^, and a^ being equal to their respective nominal values and the control feedback gains are time-varying. The optimum value, J*. of J was found to be 0.655 while the value of J* ^ ' o  a (see (6.5)) was equal to 0.705, that is, the increase in the performance index value does not exceed 7.5%. Figures (6.6), (6.7), and (6.8) show the identification results for the unknown parameters a^, a^, and a^, respectively. The maximum identification errors do not exceed 5%. The rapid rate of identification is evident in both this example and the previous one since all the unknown parameters have been identified 109 to a good degree of accuracy after only one observation interval. In Figures (6.4), and (6.5), the piece-wise constant curves which are drawn as solid lines represent the optimum constant gain K* and K*, respectively. The parameters a^, and (see (6.1)) are updated at the end of each observation interval and a new estimate of the feed back gains is made at the start of the following optimization interval. Since the optimum values of the constant feedback gains depend on the plant state at the start of the optimization interval, these optimum values are not constant during the whole operation period (see Figure (6.1-a)). The curves drawn as dashed lines in Figures (6.4), and (6.5) represent the optimum time-varying feedback gains K^, and respectively. These gains minimize the performance index given by equation (6.6) subject to the dynamical constraints given by (6.1) with a^ and a^ being equal to their respective nominal values. Figures (6.9), (6.10), and (6.11) show the optimum feedback gains K*, K* and K*, respectively. The curves drawn as solid lines represent the optimum constant gains. The parameters a^, a^, and a^ (see (6.7)) are updated at the end of each observation interval and a new estimate of the feedback gains is made at the start of the following optimization interval. The optimum constant gains are not constant along the whole operation period but unlike the situation in Example 1 the feedback gains for Example 2 (see Figures (6.9), (6.10), (6.11)) change less frequently. The curves draxm as dashed lines in Figures (6.9), (6.10), and (6.11) represent the optimum time-varying feedback gains K^, K^, and respectively. These gains minimize the performance index given by equation (6.9) subject to the dynamical constraints given by (6.7) with a^, a^, and a^ being equal to their respective nominal values. 110 The recursive observation-optimization strategy described in Section (6.2) is characterized by its flexibility. This flexibility is due to the fact that it is possible to adjust the displacement distance (£T) between the successive optimization intervals (see Figures (6.1-b, c, d)) according to the behaviour of the system under consideration. The shortest time interval during which the gains remain constant is equal to (2rT) in Example 1 (see Figures (6.4), (6.5)) and equal to (4rT) in Example 2 (see Figures (6.9), (6.10), (6.11)). Hence, using I = 2r and 4r for Examples 1 and 2, respectively, results in a satisfactory performance with less computational effort. The value of % used in both Example 1 and 2 was % = r. In the following discussion two performance indices will be discussed. These performance indices can be used to decide whether or not it is necessary to update the identified parameters and/or the con stant feedback gains at the end of each observation interval. 6.3.4 Discussion The following performance index can be used to reach decisions concerning the necessity for updating model parameters: NT 1 2 o p m p m In Eqn. (6.10) Q is a positive definite weighting matrix, xp is the measured plant state during the latest observation interval, and x is the corresponding model state during the same interval using the m latest identified values of the plant unknown parameters. The following performance index can be used together with given by (6.10) to reach decisions concerning the necessity for updating feedback gains: x'(t?)x (t?) - x*(tC ..)x (t? ) j = P 1 P 1 P X-X P ^ (6.11) In Eqn. (6.11), x (t^) Is tne plant state at the start of the i-th optimization interval (the present interval), and x (t^- -j) is the plant state at the start of the (i-l)-th optimization interval. The assignment of values for the upper bounds of and J^y respectively, would depend upon the degree of familiarity with the parti cular system under investigation (as can be seen from Figures (6.4), (6.5), (6.9), (6.10), and (6.11)). Taking IT = 2rT would be a reasonable choice for the system considered in Example 1. For Example 2 taking £T = 4rT would be a reasonable choice. Parameter or gain updating would be performed whenever or exceeded their respective bounds. 0-4. 0.2 0-0. 1 23456789 Wt 0.2 'I 0.4 . -0.5. -0.8 . -10. FIG. (6.2) EXAMPLE 1. IDENTIFICATION OF PARAMETER o, WITH QUANTIZATION NOISE. 1.0. 0.8. 0.6 . 0-4V 2 0.2 0.0 0.2 0.4 10 t FIG.(6.3) EXAMPLE I. IDENTIFICATION OF PARAMETER a2 WITH QUANTIZATION NOISE. 0.2 0.0 -0.2\ -0-4 «1 -0.6 -0-0 -1.0 -1.2 -1.4 OPTIMAL SOLUTION ADAPTIVE SOLUTION i i 1 1 1 1— 1 2 3 4 5 6 —I 1 r-7 8 9 FIG. (6.4) EXAMPLE 1. THE OPTIMUM PIECEWISE CONSTANT FEEDBACK GAIN -0.61 FIG.(5.5) EXAMPLE J. THE OPTIMUM PIECEWISE CONSTANT FEEDBACK GAIN K2. -0.2 I -0.4 -OS -as -t.o FIO.(SJS) EXAMPLE 2. IDENTIFICATION OF PARAMETER a, WITH QUANTIZATION NOISE 2.0 IS 1.2 0 8 0-4 -0-4 -0-8 10 I FI61S-7) EXAMPLE 2. IDENTIFICATION OF PARAMETER aj WITH QUANTIZATION NOISE. as. 0-4. ao. -0.4 "3 -as. -1-2 . -IS -2X1. FIGIS-S) EXAMPLE 2. IDENTIFICATION OF PARAMETER o, WITH QUANTIZATION NOISE. 02 0.0 -0.1 -0.4 OPTIMAL SOLTN . ADAPTIVE SOLTN FIG. (5.9) EXAMPLE 2. THE OPTIMUM PIECEWISE CONSTANT FEEDBACK CAIN Kt 0.4 0.2 0.0 -0.2 -0-4 -OS -08 -t-0 OPTIMAL SOLTN — — — ADAPTIVE SOLTN I234SS7S9 FIG.t6.IO) EXAMPLE 2. THE OPTIMUM PIECEWISE CONSTANT FEEDBACK GAIN Kj 1.0. OS OS 0.4 3 0.2. OO -0.2 -0.4 10 t OPTIMAL SOLTN - -ADAPTIVE SOLTN FIG.(S-U) EXAMPLE 2. THE OPTIMUM PIECEWISE CONSTANT FEEDBACK GAIN Kj. 7. CONCLUSIONS In Chapter 2 the generalized equation error has been shown to be applicable to a mean-square method of rapid digital estimation of linear system parameters. Due to the imposed structure of the estimator the manipulation of high order matrices is avoided. The matrices used by the estimator have a block diagonal structure. Consequently only the lower order matrix within a block need be considered. It has been shown that the matrix can be recursively evaluated and that matrix inversions are not required. Examples illustrate the effectiveness of the estimator for a variety of cases dealing with quanti zation noise as well as measurements noise. In Chapter 3 a simple algorithm for computing the generalized inverse of a matrix has been developed. The advantage of this algorithm over other methods is that it eliminates the need for Gram-Schmidt ortho gonalization and its associated (in the interest of accuracy) reorthogon-alization as new vectors are introduced. Inherent in the method is a very simple check on the linear dependence of the different matrix columns. In Chapter 4 a two-stage method for estimating, time-invariant and time-varying parameters in linear systems has been developed. During the second stage, the time-varying parameters are considered as unknown control inputs to a linear subsystem of known dynamics. This method for estimating time-varying parameters is computationally simple in that the estimates are obtained by the solution of a linear two-point boundary-value problem. Standard methods of solution result in one-step convergence A method of stage and control augmentation in conjunction with discreti-zation and mean-square error minimization has been shown to be effective for the first stage of estimation. This approach reduces undesirable coupling between the two-stages which could result in convergence problems. Numerous examples illustrate the effectiveness of the method. In Chapter 5 method of successive approximations for solving the two point boundary-value problem for optimum constant gain matrices has been developed. The method is shown to be computationally equivalent to a de flected gradient method. Convergence can always be achieved by choice of a scalar step-size parameter. Rapid convergence is important in off line design studies when the effect of various weighting matrices on optimal control performance is investigated. A close relationship be tween successive approximations and the solution technique for the linear time invariant case with infinite t^ is exploited so that in many cases rapid convergence is achieved. The problem of trajectory sensitivity reduction by augmented state feedback is shown to be well-posed if gain matrices are constrained to be constant. The problem of on-line implementation of a constant gain matrix whose elements depend on initial conditions is discussed. The simplest approach is to take a weighted average of the computed optimum gain matrices. This results in an average suboptimal linear feedback control law. A second approach is to introduce nonlinear feed back control [39]. An on-line approximate method is developed which appears suitable for systems whose parameters must be identified. The two point boundary-value problem is replaced by an algebraic problem, the solution of which gives a suboptimal constant gain. It does not appear feasible to give a formal proof of convergence or a proof of improved response. By proper choice of weighting matrices, however, it does seem possible to achieve an improvement with the first iteration. Examples and para meter values were chosen to deliberately accentuate situations where convergence problems could arise. No difficulties, however, were en countered and rapid convergence was achieved in all cases investigated. A natural development of the present research is to augment the sensitivity vector with the sensitivity functions with respect to the initial conditions. This will reduce the sensitivity of the optimum constant gains with respect to the state of the plant at the start of the optimization interval. Accordingly, the optimum constant gains will be updated less frequently. Further research is needed concerning the "best" choice of the constrained structure of the feedback gains. One possibility has been investigated in this thesis where the gains are taken to be time-invariant. A more general choice is to use specific time-varying gains. Let the feedback gains be constrained by the specific structure K(t) = FG(t), where the time-varying matrix G(t) is to be specified before starting the optimization process. The time-invariant gain matrix F is to be computed according to the algorithm developed in Chapter 5. The problem is how to choose G(t) systematically such that the resulting time-varying gain matrix K(t) is "optimum". A specific time-varying gain could result in better results concerning sensitivity function reduction. APPENDIX I INITIAL CONDITIONS LINEAR LEAST-SQUARES ESTIMATOR The values of the optimum constant feedback gains as computed by the algorithm developed in Chapter 5 depend on the plant state at the start of the optimization interval. Accordingly, a good estimate of the plant initial conditions is required. In a noise free situation, it may be possible to obtain the initial conditions by measuring the state of the plant. However, if the measurements of the plant state are conta minated by noise then an initial conditions estimator is required. One possible choice of such an estimator is the linear least-squares one. This estimator does not require a priori knowledge about the noise statistics and is characterized by its simplicity. Considering-the- i-ttv optimization interval (see Figures- 6vl-a, b, c, d), the instant (t^) indicates the start of this interval as well as the completion of the i-th observation interval. c c Let the plant be represented (for t^ - NT < t < t^) by x = A±x + B± u (1.1) where A^, and are the identified values of the plant parameteric matrices A, and B, respectively. These values are obtained according to the measurements taken during the i~th observation interval. The second term on the right hand side of equation (1.1) can be expressed as B.u = -F.KjX (1.2) 1 1 i where the definition of the matrix F. is as given in Chapter 5 (F = BR "*"B') and represent the piecewise constant optimum feedback gain matrix being used during the i-th observation interval. Substituting (1.2) into (1.1) yields x = (A. - F.K.)x i 11 (1.3) The transition matrix i^(t^, x) of the dynamic system represented by (1.3) can be obtained by solving the following matrix differential equation <L = -(A. - F .K. ) ' ii; iKt?, t?) = I r \ x llT T 1 1 (1.4) at Let Xj, represent the noisy measurements of the plant-state C /s t = t^ - KT, and let x^ be an estimate of x^ defined by *K = *K x(ti> (1.5) where x(t^) is the estimate of x(t^) which minimizes the quadratic error function e = -r (x - x ) ' (x - x ) 2 a a a a (1.6) where *K-*<^ 'i"1^ Xa = (1.7) The state vector x is given by a b J x = ij, x(t?) a a l (1.8) where N (1.9) The value of x(t^) which minimizes the error function e given by (1.6) is known to be x(tC) = W fJ'Vx, (1.10) 1 a a a a Equation (I.10) gives the best estimate of the plant state at t = t^ which is the initial condition for the i-th optimization interval (the present interval). These initial conditions are required in order to compute new values for the optimum constant feedback gains. APPENDIX II ANALYTICAL SOLUTION FOR THE LIAPUNOV MATRIX EQUATION [42] This algorithm is characterized by its simplicity since no matrix inversion is needed and the solution is obtained in a closed form and not in a numerical or iterative way. Consider the matrix equation S'X + XS = -Q (II.1) where S, X and Q are all of dimension (n x n) Let A k k-i S, = Z X.S (II.2) k i=0 1 where the parameters X^(i = 0, n-l) are the coefficients of the chracteristic equation of the matrix S. These parameters can be computed easily using Faddeev's method. X. = - ~ Tr[SS. .], S. = SS. . + X.I, i = 1, n 1 x i-l ' l i-l x ' ' ' XQ = 1 SQ = I, S = 0 (II.3) The last relation is used as a check. We may write down the solution to equation (II.1) as 1 2n"2 k k-l X = jj-i E A (S, -S) E s;Q(-l)k *S (II.4) (-l)nR(S, -S) k=0 K 1 £=0 % k % 123 where: R(S, -S) = (-l)n2n T°- H2(S) n (II.5) Hn(S) = A, A 1 o A A A0 Ai A 3 2 1 o 0 0 n (II.6) The determinant H (S) is known as the Hurwiz determinant, r o AK+1(S, -S) k odd _ H (S) l» 2n —7 G, (S) k even A k n (II.7) where G^, for k even, is the (- + l)th. Cofactor of the first column of H (S). n Hence equation (II,. 4) can be written in a simplified form as X = 1 n_1 21 I E G2±(S) E (-iyS\QSn4 2*^ (S) i=0 n £,=0 2i-£ (II.8) From (II.8) it is clear that the condition for a unique solution is H (S) ^ 0 which- is equivalent, to. the condition, that, the dynamical system. (X = SX) is asymptotically stable. For the non-trivial case where the matrix Q is symmetrical, the term 21 I E (-D S;QS A=O (II.9) 124 will be symmetrical and since G„.(S) and H (S) are scalers, X is also J - 2i n symmetrical and we need only to solve n(n + l)/2 entries of X instead 2 of n . Thus for this case we have XU=-rT-i VG (S)Z <-!>*<SjQS ) ' (11.10) 11 2N H (S) 1=0 L Si=0 % lX % U where is the upper triangle of X and (s1pQS2i-X?]J is trie uPPer triangle for (SlQS„. „•) for the case i = £, while for i ^ I it is a triangular matrix with its main diagonal equal to twice the main diagonal of (SlQS ) and every element (i, j) is equal to the sum of the corresponding element (i, j) in the upper triangle of (S£Qs2i_£) Plus the (i + 1, j - 1) element in the lower triangle. 125 APPENDIX III The trace properties and gradient matrices used in deriving the gradient conditions are the following [43], Tr[AB] = Tr[BA], Tr[ABf] = Tr[BA1] Tr[AKB] = A'B' Tr[AK'B] = BA Tr[AKBK] = A'K'B' + B'K'A' dK Tr[AK'BK'] = BK'A + AK'B dK. ^ TrfAKBK'] = A'KB' + AKB (III.l) dK The following properties of Kronecker matrix products have been utilized in the algebraic manipulation of matrix equations [44], (AB)v = (10 A) (B)v = (B' © I) (A)v (A ® B) (C © D) = (AC) © (BD) (A + B) ® (C + D) = A®C + A®D + B®C + B®D (III. 2) If Aj,, i = 1, .. ., n, and u , j = 1, n, are the eigenvalues of A and B, respectively, then u^. are the eigenvalues of A@ B. Because of the partitioning (5.7) of F, the eigenvalues of F can be found from the eigenvalues of F. Consider the following eigenvalue equations: Fz. = X.z. x x x (111,3) $ V. = U.V. yy J ii It is seen from the defining equations for F and $ that (see (5.5) and (5.12)) (B'z.)'R 1 (B'z.) = Xizi' z. = 0 (III.4) 2 / (v.'y) dt = u.v.'v. = 0 o 1 111 The matrix W is therefore positive semi-definite. Since ("j^) = W(S - K)^ it follows that the equality sign in (5.28) holds only when the gradient condition is satisfied. 127 APPENDIX IV The confidence interval length deduced in Subsection(2.7.3) (see (2.101)) is based on the assumptions that the random vector (c - c) is normally distributed and that the random vector V is statistically r. independent of (c - c) and is normally distributed. Also the random 2 variable V'V is assumed to be distributed as % . The above mentioned assumptions can only be justified if the matrix 6^ (see (2.11)) is deterministic. In the following we will prove the above assumptions for the case where 6^ is a deterministic quantity: Using (2.12) and (2.38) yields c - c = PM0'ZM - c N N N = PN6N(6NC + N) " C = PN6^N (IV.1) where, for convenience, the random vector N is taken to be normally distributed with E[N] = 0 and E[N N'] =1. It follows directly from (IV.1) that (c - c) is a normally distributed random vector. Using (2.63) and(EV.l) yields . E[(c - c)V] = E[PN0^NN'(6NPN6^ - I)] = PN6N(9NPNeN ~ I} =0 (IV.2) Thus the random vectors (c - c) and V are statistically inde pendent. The vector V is normally distributed as can directly be seen from (2.63) where 128 V = (0NPN N " I)S = (WN " ^ (2'63) Since 8JJ and P^ are deterministic matrices then V is normally distributed as N. Using (2.63), (2.64), and (2.65) yields W = N'(WN - D*(WN - I)N = N'(WN - I)2N = Nr(I - WN)N (IV.3) Considering the orthogonal matrix <J>^ (see (2.66)), equation (IV.3) can be written as V'V = N$ $'(I - W )$ $'N V V VNV V N N N'$N(I - DN)^N = EyJ (IV.4) where y. is a linear combination of all the elements of the random vector (the matrix (I - D^) is a diagonal one with the main diagonal elements being equal to 1 or 0 (see (2.69))). Thus the random vector y^ is distributed as the vector $' N. The vector $'N is distributed as N N the random vector N which is a direct result of Fisher's Theorem [17]. Accordingly, the random vector is normally distributed like N and 2 2 Ey^ is distributed as % . The above proofs are based on a deterministic 0^. If 0^ is a random variable, it is no longer possible to prove that the random vectors (c - c) and V are normally distributed along with the assumption that VrV is distributed is % . There appears to be no mathe matically tractable approach to derive confidence intervals in the case 129 where 0^ contains quantization noise. To retain the basic simplicity of Equation (2.101) the following procedure is suggested. In equations (IV.1 -4), it is possible to replace ZN(k) by ZN(k) where ^N(k) is obtained by use of the best estimate c at stage k, that is by use of (see (2.10)). zN = ?N; dv.5) The "smoothed" values 0^ of 0^ are used in equations (IV.1-4). Since the randomness in 6^ tends to be eliminated, it seems justifiable to consider 0^ as a deterministic quantity. APPENDIX V The gradient condition as given by equation (5.11) can be derived as follows: Let H denotes the Hamiltonian, then using (5.6) and (5.9) yields H = p'(A - F(K + K'))y - | y»(Q + (K + K')F(K + K'))y (V.l) = p'Ay - p'FKy - p'FK'y - | y'Qy - | y'KFKy - j y'KFK'y - -| y'K'FKy - j y'K'FK'y Using (III.l) the Hamiltonian can be written as H = p'Ay - TrfFKyp'] - Tr[FK'yp'] - ~ y'Qy - |- Tr[yy'KFK] - | Tr[yy'KFK'] - | Tr[FKyy'K'] - |- Tr[yy'K'FK'] (V.2) The following gradient matrix operations are required g| Tr[FKyp'] = F(py') ^ Tr[FK'yp'] = (yp')F " \ ~k Trtyy'KfK] = " f(yy')K'F - \ FK' (yy') - \ g| Tr[yy'KFK'] = -(yy')KF - Tr[FKyy'K'] = -FK(yy') (v.3) " \ a! Tr[yy'K'fK'] = - \ ^'(yy') - f(yy')K'F The above operations allow ~ to be written as 13 where 3H 1 1 — = - F(py») - (yp')p - |(yy«)K'F - -| FK»(yy*) - (yy')KF - FK(yy') - | FK'(yy') - | (yy')K'F = - F(py') - (yp')F - (yy*)K'F - FK'(yy') - (yy')KF - FK(yy«) = - F(py') - (yp*)F - (yy')(K + K')F - F(K + K,)(yy') = - F(py') - (yp')F - (yy')KF - FK(yy*) (V.4) K = K + K? From (V.4) and for time-invariant F and K equation (V.4) gives t aH _ t t t t £ - 9K = F / py'dt + / yp'dt F + / yy'dt KF + FK / yy'dt = 0 000 o (V.5) Using the definition A t & = / yv'dt y v J O Equation (V.5) can be written as t gxr _ _ _ / - -r^r = F<£> + $ F + $ KF + FK $ =0 (V.6) q 9K py yp yy yy This completes the proof of equation (5.11). 132 REFERENCES 1. P. V. Kokotovic', J. B. Cruz, Jr., J. E. Heller and P. Sannuti, "Synthesis of Optimally Sensitive Systems", Proc. I.E.E.E., Vol. 56, No. 8, August, 1968. 2. Lion, P., "Rapid Identification of Linear and Nonlinear Systems", Joint Automatic Control Conference, August, 1966, pp. 605-615. 3. Pazde ra, J., and Pottinger, H., "Linear System Identification Via Liapunov Design Techniques", Joint Automatic Control Conference, August,1969, pp. 795-801. 4. Rose, R., and Lona, G., "An Approximate Steepest Descent Method for Parameter Identification", Joint Automatic Control Conference, August, 1969, pp. 483-487. 5. Kishi, F. H., "On-Line Computer Control Techniques and Their Appli cation to Re-entry Aerospace Vehicle Control", C. T. Leondes, ed. Advances in Control Systems Theory arid Applications, Vol. 1, Academic Press, New York, 1964. 6. C. Price, "The Matrix Pseudo-Inverse and Minimal Variance Estimates',' SIAM Rev. 6, 1964, pp. 115-120. 7. L. D. Pyle, "Generalized Inverse Computations Using the Gradient Projection Method", Journal of ACM, 11, 1964, pp. 422-428. 8. B. Rust, W. R. Burrus, C. Schneeberger, "A Simple Algorithm for Computing the Generalized Inverse of a Matrix", Journal of ACM, 9, 1966, pp. 381-387. 9. Sage, A., and Melsa, J., Systems Identification, Academic Press, 1971. 10. H. C. Lessing and D. F. Crane, "The Use of Integral Transforms in the Estimation of Time Variable Parameters", NASA Techn. Note, NASA TN D-5194, May 1969. 11. Sage, A. P., Optimum Systems Control, Prentice-Hall, Inc., 1968, pp. 507. 12. Mayne, D. 0., "Optimal Non-Stationary Estimation of the Parameters of a Linear System with Gaussian Inputs", Journal of Electronics  and Control, Vol. 14, 1963, pp. 101--112. 13. Lee, R. C. K., Optimal Estimation and Control, M.I.T. Press, Cambridge, Mass., 1964. 14. Lee, R. C. K., and Ho, Y. C, "Identification of Linear Dynamic Systems", Information arid Control, Vol. 8, No. 1, February, 1965. 15. Leathrum, J., "On-Line Recursive Estimation", Joint Automatic Control Conference, August, 1969, pp. 173-178. 16. Bakke, R., "Case Study in the Paper-Making Industry", Case Studies in System Control, 1968, pp. 113-192. 17. Y. V. Linnik, Method of Least Squares and the Principles of the  Theory of Observation, Pergamon, New York, 1961. 18. M. Aoki, Optimization of Stochastic Systems, Academic Press, 1967, pp. 318-324. 19. B. J. Prochaska, "Applications of Generalized Inverses to Esti mation in a Dynamic System", JACC, 1971, pp. 820-829. 20. D. Luenberger, Optimization by Vector Space Methods 9 John Wiley and Sons, 1969. 21. J. S. Meditch, "A Class of Suboptimal Linear Controls", Trans, on Automatic Control, July 1966, pp. 433-439. 22. A. Ben-Israel, A. Charnes, "Contributions to the Theory of Genera lized Inverses", SIAM Rev. 11, 1963, pp. 667-699. 23. R. Penrose, "A Generalized Inverse for Matrices", Proc. Cambridge Philosophical Society, 51, 1955, pp. 406-413. 24. D. K. Faddeev and V. N. Faddeeva, Computational Methods of Linear  Algebra, W. H. Freeman and Co., 1963, 25. D. L. Kleinman and M. Athans, "The Design of Suboptimal Linear Time-varying Systems", I.E.E.E. Trans., Automatic Control, Vol. AC-13, pp. 150-159, April 1968. 26. D. L. Kleinman, T. Fortmann and M. Athans, "On the Design of Linear Systems with Piecewise Constant Feedback Gains", I.E.E.E. Trans., Automatic Control, Vol. AC-13, pp. 354-361, August 1968. 27. W. Levine and M. Athans, "On the Determination of the Optimal Constant. Output. Feedback..Gains for Linear Multivar.iab.le.. Systems", I.E.E.E. Trans., Automatic Control, Vol. AC-15, pp. 44-48, February 1970. 28. M. M. Newmann, "Specific. Optimal Control of the Linear Regulator Using a Dynamical Controller Based on the Minimal-Order Luenberger Observer", Intern. J. Control, Vol. 12, pp. 33-48, No. 1, 1970. 29. R. L. Kosut, "Suboptimal Control of Linear Time-Invariant Systems Subject to Control Structure Constraints", Proc. JACC, pp. 820-828, 1970. 30. T. Hendericks and H. D'angelo, "An Optimal Fixed Control Structure Design with Minimal Sensitivity for a Large Electric Booster", Proc. of the 5th Annual Allerton Conference, pp. 142-148, 1967. 134 31. E. Kreindler, "On Minimization of the Trajectory Sensitivity", Intern. J. Control, Vol. 8, No. 1, pp. 89-96, 1968. 32. E. Kreindler, "Synthesis of Flight Control Systems Subject to Vehicle Parameters Variations", Grumman Aircraft Engg. Corp., Bethpage, N. Y., Techn. Rept. AFFDL-TR-66-209, April 1967. 33. E. Kreindler, "Formulation of the Minimum Trajectory Sensitivity Problems", I.E.E.E. Trans., Automatic Control, Vol. AC-14, pp. 206-207, April 1969. 34. J. F. Cassidy Jr. and I. Lee, "On the Optimal Feedback Control of a Large Launch Vehicle to Reduce Trajectory Sensitivity", Proc. of 1967 JACC, pp. 587-595. 35. H. J. Dougherty, I. Lee and P. M. DeRusso, "Synthesis of Optimal Feedback Control Systems Subject to Parameter Variations", 1967 Proc. JACC, pp. .125-133. 36. P. Sannuti and J. B. Cruz Jr., "A Note on Trajectory Sensitivity of Optimal Control Systems", I.E.E.E. Trans., Automatic Control, Vol. AC-13, pp. 111-112, February 1968. 37. A. J. Bradt, "Sensitivity Functions in the Design of Optimal Controllers", I.E.E.E. Trans., Automatic Control, Vol. AC-13, pp. 110-111, February 1968. 38. D. L. Kleinman, "On an Iterative Technique for Riccati Equation Computations", I.E.E.E. Trans., Automatic Control, Vol. AC-13, . pp. 114-115, February 1968. 39. A. G. Longmuir and E. V. Bohn, "The Synthesis of Suboptimal Feed back Control Laws", I.E.E.E. Trans., Automatic Control, Vol. AC-12, pp. 755-758, December 1967. 40. C. F. Chen and L. S. Shiek, "A Note on Expanding PA + A'P = -Q", I.E.E.E. Trans., Automatic Control, Vol. AC-13, pp. 122-123, February 1968. 41. Ii. K. Kirchmayer, Economic Control of Interconnected Systems , New York, John Wiley and Sons, 1959. 42. P. Mueller, "Solution of the Matrix Equations AX + XB = -Q and S'X + XS = -Q*", SIAM J. Appl. Math., Vol. 18, No. 3 pp. 682-687, May 1970. 43. M. Athans and F. C. Schweppe, "Gradient Matrices and Matrix Calculations", M.I.T. Lincolin Lab., Lexington Mass. Techn. Note 1965-53, November 1965. 44. H. Neudecker, "A Note on Kronecker Matrix Products and Matrix Equation Systems", SIAM J. Appl. Math., Vol. 17, No. 3, pp. 603-606, May 1969. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

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>
                        
                    
IIIF logo 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-0101279/manifest

Comment

Related Items