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

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

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
UBC_1972_A1 H35.pdf [ 5.83MB ]
Metadata
JSON: 1.0101279.json
JSON-LD: 1.0101279+ld.json
RDF/XML (Pretty): 1.0101279.xml
RDF/JSON: 1.0101279+rdf.json
Turtle: 1.0101279+rdf-turtle.txt
N-Triples: 1.0101279+rdf-ntriples.txt
Original Record: 1.0101279 +original-record.json
Full Text
1.0101279.txt
Citation
1.0101279.ris

Full Text

LINEAR SYSTEMS IDENTIFICATION AND OPTIMIZATION WITH APPLICATION TO ADAPTIVE CONTROL  ADEL ABDEL-RAOUF HANAFY B.Sc,  U n i v e r s i t y o f C a i r o , C a i r o , Egypt, 1965  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE  REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  i n the Department of Electrical  Engineering  We a c c e p t t h i s t h e s i s as conforming t o the required  THE  standards  UNIVERSITY OF BRITISH COLUMBIA J u l y , 1972  In p r e s e n t i n g t h i s t h e s i s  in p a r t i a l  f u l f i l m e n t o f the r e q u i r e m e n t s  an advanced degree at the U n i v e r s i t y of B r i t i s h C o l u m b i a , the L i b r a r y  shal1 make i t f r e e l y  available for  I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e  I agree  r e f e r e n c e and copying of t h i s  representatives.  It  is understood  that copying or  thesis  t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my  written  permission.  The U n i v e r s i t y o f B r i t i s h Columbia Vancouver 8, Canada  Date  26-1-7 Z  or  publication  of t h i s  Depa rtment  that  study.  f o r s c h o l a r l y purposes may be g r a n t e d by the Head of my Department by h i s  for  ABSTRACT  T h i s t h e s i s i s concerned w i t h the problem of i d e n t i f y i n g c o n t r o l l i n g l i n e a r continuous for  systems.  Algorithms  r e a l - t i m e d i g i t a l computations are developed  the c o n t r o l l e r and The  the  which a r e  feasible  f o r the d e s i g n o f  Due  g e n e r a l i z e d e q u a t i o n e r r o r i s shown to be a p p l i c a b l e t o a  to the imposed s t r u c t u r e of the e s t i m a t o r the  of h i g h o r d e r m a t r i c e s  i s avoided.  Examples i l l u s t r a t e  ness of the e s t i m a t o r f o r a v a r i e t y of cases n o i s e as w e l l as measurements n o i s e .  the  effective-  dealing with quantization  In some c a s e s ,  this d i g i t a l  identi-  matrix.  a l g o r i t h m f o r computing the g e n e r a l i z e d i n v e r s e of a m a t r i x i s  developed.  T h i s a l g o r i t h m e l i m i n a t e s the need f o r Gram-Schmidt  g o n a l i z a t i o n and  i t s a s s o c i a t e d ( i n the i n t e r e s t of a c c u r a c y )  g o n a l i z a t i o n as new  v e c t o r s are  time-varying  ortho-  reortho-  introduced.  A two-stage e s t i m a t o r i s developed and  para-  manipulation  f i e r r e q u i r e s the computation of the g e n e r a l i z e d i n v e r s e of a A simple  both  identifier.  mean-square method of r a p i d d i g i t a l e s t i m a t i o n of l i n e a r system meters.  and  for estimating time-invariant  parameters i n l i n e a r systems.  During  the second  stage,  the t i m e - v a r y i n g parameters are c o n s i d e r e d as unknown c o n t r o l i n p u t s to a l i n e a r subsystem of known dynamics. identifier  F o r the f i r s t  s t a g e , the  digital  i s shown t o be e f f e c t i v e i n i d e n t i f y i n g the t i m e - i n v a r i a n t  parameters.  Numerous examples i l l u s t r a t e  the e f f e c t i v e n e s s o f the method.  To d e s i g n a feedback c o n t r o l l e r a method of s u c c e s s i v e a p p r o x i mations f o r s o l v i n g the two constant  gain matrices  p o i n t boundary v a l u e problem f o r optimum  i s developed.  The method i s shown to be  t i o n a l l y e q u i v a l e n t to a d e f l e c t e d g r a d i e n t method.  ii  Convergence  comptacan  always be a c h i e v e d by c h o i c e o f a s c a l a r s t e p - s i z e parameter. l i n e approximate method i s developed whose parameters must be i d e n t i f i e d .  An  on-  which appears s u i t a b l e f o r systems The  two p o i n t boundary v a l u e problem  i s r e p l a c e d by an a l g e b r a i c problem, the s o l u t i o n o f which g i v e s a optimal constant gains.  The  sub-  problem o f t r a j e c t o r y s e n s i t i v i t y r e d u c t i o n  by augmented s t a t e feedback i s shown to be w e l l - p o s e d s t r u c t u r e of the g a i n m a t r i c e s i s taken.  The  i f constrained  simple s t r u c t u r e o f  constant  gain matrices i s considered. Based on the developed  i d e n t i f i e r and  c o n t r o l l e r , one p o s s i b l e  s t r a t e g y f o r o p t i m a l a d a p t i v e c o n t r o l which i s p a r t i c u l a r l y from an e n g i n e e r i n g p o i n t of view i s s t u d i e d .  The  attractive  s t r a t e g y i s to  i d e n t i f y system parameters at the end of an o b s e r v a t i o n i n t e r v a l then to use  the parameters to d e r i v e an " o p t i m a l " c o n t r o l f o r a sub-  sequent c o n t r o l i n t e r v a l . illustrate  and  Two  examples are c o n s i d e r e d i n o r d e r t o  the e f f e c t i v e n e s s of t h i s s t r a t e g y .  iii  TABLE OF CONTENTS Page ABSTRACT  i i  TABLE OF CONTENTS  iv  LIST OF TABLES  vi  LIST OF FIGURES  v i i  ACKNOWLEDGEMENT  ix  1.  INTRODUCTION  1  2.  RAPID DIGITAL IDENTIFICATION OF LINEAR SYSTEMS 2.1 2.2 2.3 2.4 2.5 2.6 2.7  2.8 3.  8 9 14 17 18 27 30  2.7.1 2.7.2 2.7.3  30 36 37  Correlation o f Noise Samples The Error Covariance Matrix Estimation of the Confidence Interval Length  Discussion  42  RECURSIVE COMPUTATION OF THE GENERALIZED INVERSE OF A MATRIX. 3.1 3.2 3.3 3.4 3.5  4.  Introduction Problem Formulation Recursive and Block Data Processing Measurement Requirements Applications Comparison with a Limited Memory F i l t e r Some S t a t i s t i c a l Properties of the D i g i t a l I d e n t i f i e r . . .  Introduction Properties of the Pseudo-Inverse Recursive Evaluation of R^ Recursive Evaluation of A Examples  44 44 46 49 50  1  TWO-STAGE ESTIMATION OF TIME-INVARIANT AND TIME-VARYING PARAMETERS IN LINEAR SYSTEMS 4.1 4.2 4.3 4.4 4.5 4.6 4.7  Introduction.. Problem Formulation Estimation of A . Algorithm 1 Estimation of a P e r i o d i c \ Algorithm I I Estimation of a P e r i o d i c A . Algorithm I I I Two-Stage Estimation of Af and A Examples  56 57 59 62 65 67 69  4.7.1 4.7.2  70 71  v  v  v  Example 1 Example 2  Page 4.7.3 4.7.4 4.7.5 4.7.6 5.  71 72 73 73 THE  Introduction Problem F o r m u l a t i o n U n i f i e d Treatment o f D i f f e r e n t F o r m u l a t i o n s G r a d i e n t Descent and S u c c e s s i v e Approximations Optimum Constant G a i n : T r a j e c t o r y S e n s i t i v i t y Reduction. On-Line Computation o f Gains Examples  76 78 82 84 88 90 93  OPTIMAL ADAPTIVE CONTROL 6.1 6.2 6.3  Introduction Optimum C o n t r o l D u r i n g the F i r s t Optimization Intervals Examples 6.3.2 6.3.3 6.3.4  7.  3 4 5 6  OPTIMUM CONSTANT FEEDBACK GAINS FOR LINEAR SYSTEMS AND PROBLEM OF TRAJECTORY SENSITIVITY REDUCTION 5.1 5.2 5.3 5.4 5.5 5.6 5.7  6.  Example Example Example Example  102 and Subsequent 102 104  Example 1 Example 2 Discussion  .  CONCLUSIONS  105 107 110 116  APPENDIX I  119  APPENDIX I I  .'  .  122  APPENDIX III.-  125  APPENDIX IV  127  APPENDIX V  130  REFERENCES  •  132  v  LIST OF TABLES  Pagi  T a b l e (4.1)  75  T a b l e (4.2)  75  T a b l e (4.3)  75  T a b l e (4.4)  75  LIST OF FIGURES Figure  Page  1.1  The Proposed General Optimal Adaptive System  2.1  Example 1.  Estimation of a^ with Quantization Noise  21  2.2  Example 1.  Estimation of  21  2.3  Example 1. Estimation of a^ with Quantization and Unbiased Measurement Noise  22  Example 1. Estimation of a^ with Quantization and Unbiased Measurement Noise  22  Example 1. Estimation of a-^ with Quantization and Biased Measurement Noise  23  Example 1. Estimation of a^ with Quantization and Biased Measurement Noise  23  2.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  Example 2. Estimation of a^ with Quantization and Unbiased Measurement Noise  26  Example 2. Estimation of b^ with Quantization and Unbiased Measurement Noise  26  Example 2. Estimation of a., with Quantization and Unbiased Measurement Noise Using a Standard Limited Memory F i l t e r . . .  28  Example 2. Estimation of a~ with Quantization and'Unbiased' Measurement Noise Using a Standard Limited Memory F i l t e r . . .  28  Example 2. Estimation of b with Quantization and Unbiased Measurement Noise Using a Standard Limited Memory F i l t e r . . .  29  4.1  The Recursive Observation-Estimation Intervals  64  5.1  Example 1. Control as a Function of Time For: (a) Method. (b) Constant Gain and Augmented States. (c) Constant Gain and Non-Augmented States  2.4  2.5  2.6  2.11  2.12  2.13 2.14 2.15  5.2  2  with Quantization Noise  1  Lee's 98  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 Example 1. Trajectory S e n s i t i v i t y Function z.^ f o r : (a) Lee's Method. (b) Constant Gain and Augmented States. (c) Constant Gain and Non-Augmented States  100  Example 2. Block Diagram of a Sychronous Generator Voltage Control System  101  6.1-a  The I n i t i a l i z a t i o n Interval and the Operation Period  103  6.1-b  The Successive Observation Intervals  103  6.1-c  The Successive Optimization Intervals (X, = r)  103  6.1-d  The Successive Optimization Intervals (I = 2r)  103  6.2  Example 1. Noise  I d e n t i f i c a t i o n of Parameter a^ with Quantization., 112  Example 1. Noise  I d e n t i f i c a t i o n of Parameter a^ with Quantization  6.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. Noise  I d e n t i f i c a t i o n of Parameter a  Example 2. Noise  I d e n t i f i c a t i o n of Parameter a^ with Quantization  Example 2. Noise  I d e n t i f i c a t i o n of Parameter a^ with Quantization  6.9  Example 2.  The Optimum Piecewise Constant Feedback Gain K^...  6.10  Example 2.  The Optimum Piecewise Constant Feedback Gain  6.11  Example 2.  The- Optimum Piecewise Constant Feedback Gain K,.-..- 115  5.3  5.4  6.3  6.7 6.8  112  with Quantization 114 114  114  viii  115  l^...  ACKNOWLEDGEMENT  I am indebted to my supervisor, Dr. E. V. Bohn, for h i s  continuous  guidance and assistance throughout the period of t h i s research. Thanks are due to Dr. V. G. Haussmann f o r h i s valuable comments on the d r a f t . i s duely  Reading the draft by Dr. M. S. Davies and Dr. G. F. Schrack  appreciated. 1 also wish to thank Miss Linda Morris for her patience i n  typing the manuscript. The f i n a n c i a l 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 i s a control that w i l l cause a plant to operate i n an "optimum" manner without a P r i o r i knowledge about the constraints on e i t h e r the magnitude or the rate of v a r i a t i o n of the plant uncertainties.  These uncertainties may  be i n the plant parameters,  i n i t i a l conditions, and/or noisy measurements.  One  possible mathematical  design approach for an adaptive control i s to use s t a t i s t i c a l methods. Plant uncertainties are then described by p r o b a b i l i t y density Such s t a t i s t i c a l methods are generally not computationally  functions.  f e a s i b l e even  when a priori knowledge of the s t a t i s t i c s of the plant uncertainties e x i s t . In most p r a c t i c a l s i t u a t i o n s , s i m p l i f y i n g assumptions must be and these detract considerably  from the usefulness  introduced  of a mathematically  precise s t a t i s t i c a l approach. An a l t e r n a t i v e approach i s the optimally sensitive control [1] where no a p r i o r i knowledge of s t a t i s t i c a l parameters i s required.  However,  the range of plant parameter and/or i n i t i a l condition variations i s restricted. This thesis i s concerned with the development of an  optimally  s e n s i t i v e feedback control which does not require detailed s t a t i s t i c a l data for i t s r e a l i z a t i o n . However, i n order to increase the range of optimal control action, parameter estimation i s used.  Algorithms which  are f e a s i b l e for real-time d i g i t a l computations are developed for the design of both the c o n t r o l l e r and the Figure (1.1) i l l u s t r a t e s  identifier.  the general structure of the optimal  adaptive control system which has been developed i n t h i s t h e s i s .  The  -n LEAST Sq. I.C ESTIMATOR  A-  NOISE 11  PLANT  OUTPUT  DIGITAL IDENTIFIER FEEDBACK Y- - 1 SENSITIVITY L_ M?Z?£Z. ^ GAINS ^  ADJUSTABLE C = 3 MODEL ~—7>  /7/.7F/P  11  r 11  'L  -SCOMPOSITE T ^JDENTIFI^ _p  FIG. (I.I)  THE PROPOSED GENERAL OPTIMAL ADAPTIVE  SYSTEM,  3  blocks drawn i n dashed l i n e s represent elements which may not be required i n many p r a c t i c a l a p p l i c a t i o n s . generally be divided into  This optimal adaptive control system can  two main parts; the estimator and the c o n t r o l l e r .  A d e t a i l e d discussion of the components which comprise these two parts together with the thesis layout i s considered In the following two sections.  1.1  The Estimator  The estimator i s used to obtain a "best" estimate for the unknown plant parameters and i n i t i a l conditions.  I t i s also used to obtain a "best"  estimate f o r the noisy and/or the unmeasured plant state v a r i a b l e s . These estimates are needed i n order to construct a reasonably accurate deterministic l i n e a r model of the plant.  A deterministic l i n e a r model i s  necessary i f use i s to be made of the powerful tools that have been developed from optimal control theory. this thesis.  Linear continuous models are considered i n  The choice of l i n e a r models i s made because of t h e i r r e l a t i v e  s i m p l i c i t y from the instrumentation and computation point of view.  This  s i m p l i c i t y i s largely due to the high degree of development of the theory of l i n e a r systems.  Even when the plant i s nonlinear, however, the use of  q u a s i l i n e a r i z a t i o n techniques can result i n a l i n e a r i z e d model of the plant. Over a fixed period of time-, such, a linear... model is. of ten. adequate. to. describe the dynamic behaviour of a nonlinear plant. The estimator i s composed of four major elements (see Figure (1.1)). These four elements are:  a- I n i t i a l Conditions Estimator  This estimator i s a standard l i n e a r least-squares estimator which minimizes  4  a quadratic error function based on a system of l i n e a r algebraic equations. It i s necessary  to obtain a good estimate of the plant i n i t i a l conditions  since the values of the c o n t r o l l e r constant feedback gains depend on i n i t i a l conditions.  This estimator i s 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 i t s s i m p l i c i t y .  Also, no  a p r i o r i knowledge of noise s t a t i s t i c s i s needed.  b - Filter  The f i l t e r gives a continuous  estimation of the plant states.  I f some  of the states are unmeasurable then a Luenberger observer i s a good choice to reconstruct the unavailable states.  On the other hand a Kalman  f i l t e r can be used for best estimation of the noisy states.  c - Digital Identifier A continuous-time  parameter i d e n t i f i c a t i o n method that has been extensively  investigated i s the parameter-tracking  technique  on a steepest descent parameter adjustment.  [2,3,4],  These  This i s based  parameter-tracking  i d e n t i f i e r s have been largely developed for use with analog computers. P r a c t i c a l on-line i d e n t i f i c a t i o n , however, requires the use of d i g i t a l computers and the d i s c r e t e version of parameter-tracking several d i f f i c u l t i e s .  introduces  The optimum gain or optimum step-size depends on  the input signals and i s d i f f i c u l t to determine on-line [5]. In Chapter 2 a d i g i t a l i d e n t i f i e r i s developed which gives rapid and accurate estimation of the unknown, time-invariant parameters of the plant.  The i d e n t i f i e r uses data obtained from augmented input  controls and output states which r e s u l t s i n p a r a l l e l data processing.  5  A d i s c r e t i z e d version of the continuous augmented model i s used. simulation r e s u l t s have shown that the i d e n t i f i e r developed has very good rates of convergence.  Digital  i n this thesis  This d i g i t a l i d e n t i f i e r requires  the computation of a matrix inverse which may not exist i n some cases. In such cases a generalized-inverse of a matrix i s required.  A number  7  of methods have been proposed f o r 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 l i n e a r algebraic equations. In Chapter 3 an algorithm i s developed which replaces these operations and computes- the generalized inverse from simple recursive formulas.  Inherent i n the method i s a simple check on the l i n e a r inde-  pendency of the column vectors of the matrix. d - Composite I d e n t i f i e r A two-stage estimator i s developed  i n Chapter 4.  In the f i r s t stage the  d i g i t a l i d e n t i f i e r i s used to obtain an estimate f o r the time-invariant parameters.  In the second stage the unknown time-varying parameters are  considered as c o n t r o l inputs and a l i n e a r state regulator quadratic cost function dynamic optimization problem i s formulated. associated two-point  The solution of the  boundary-value problem f o r the optimum control"  r e s u l t s i n an estimate for the time-varying parameters.  This two-stage  i d e n t i f i c a t i o n algorithm i s s u i t a b l e for the cases where only a few of the unknown parameters of the plant are time-varying. complicated techniques developed  The use of the  for the i d e n t i f i c a t i o n of time-varying  parameters i s not j u s t i f i e d f o r such cases.  The usual approach taken  to i d e n t i f y time-varying parameters i s to consider them as unknown augmented states described by d i f f e r e n t i a l equations of known form and then use some form of state estimation [9]. transforms have been used [10]  Alternatively, integral  to extend the c a p a b i l i t y of parameter-  tracking system to track time-varying parameters.  The implementation  d i f f i c u l t i e s for such techniques make t h e i r use questionable i f not completely u n p r a c t i c a l .  1.2  The Controller  In Chapter 5 the design of a c o n t r o l l e r with optimum constant feedback gains i s discussed.  A method of successive approximations i s  developed and i t i s shown to be equivalent to a deflected gradient descent approach with a variable weighting matrix.  This property i s u t i l i z e d  to develop an algorithm f o r o f f - l i n e design studies which r e s u l t s i n rapid convergence.  An approximate method i s 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 i n Figure (1.1)).  Trajectory  s e n s i t i v i t y function- minimization i s an important requirement... The i n troduction of the s e n s i t i v i t y functions i n the feedback structure of the adaptive c o n t r o l l e r makes i t possible to up-date the constant feedback gains less frequently and makes the plant c o n t r o l l e r less s e n s i t i v e to modelling e r r o r s .  The s e n s i t i v i t y functions are generated by the sensi-  t i v i t y model shown by f i n Figure (1.1).  The parameters of the s e n s i t i v i t y  model are up-dated each time a new estimate of the unknown plant parameters i s made.  In Chapter 6 a recursive observation-optimization strategy i s described.  During a f i n i t e observation I n t e r v a l , measurements of the  state are taken and used to determine estimates of the plant parameters. The l i n e a r model parameters are adjusted according to these new  estimates.  The l i n e a r model i s used to compute a control strategy which i s applied i n a subsequent control i n t e r v a l .  This scheme seems p a r t i c u l a r l y a t t r a c t i v e  for adaptive systems. In summary, the objective of this thesis i s to investigate p r a c t i c a l methods f o r designing an optimal adaptive control system based on a reasonably large domain of plant uncertainties.  Emphasis i s placed  on techniques which can be implemented on real-time d i g i t a l 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 d i g i t a l i d e n t i f i e r f o r continuoustime systems and by developing simple algorithms for the computation optimally s e n s i t i v e feedback gains.  of  8  2.  RAPID DIGITAL IDENTIFICATION OF LINEAR SYSTEMS  2.1  Introduction  D i g i t a l control of adaptive systems i s f e a s i b l e only  if  parameter i d e n t i f i c a t i o n and evaluation of control laws can be p e r formed i n r e a l time.  To reduce the complexity of the computations,  s i m p l i f y i n g assumptions must be made.  One f e a s i b l e method i s to  identify  the system as an approximate l i n e a r system within a given time i n t e r v a l and then use a l i n e a r optimal control [5,11],  For such a control policy  to be e f f e c t i v e i n an adaptive system the i d e n t i f i c a t i o n should be rapid. A continuous-time parameter i d e n t i f i c a t i o n method that has been extensively investigated i s the parameter tracking technique [2,3,4].  This i s 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 e r r o r .  The generalized error approach i s  attractive  since i t results i n an exact continuous-time steepest descent adjustment. Furthermore, by use of state variable f i l t e r s ,  i t allows considerable  f l e x i b i l i t y i n the generation of additional input-output additional- equation e r r o r s .  signals and  These signals play an important:- role- in the  estimator that w i l l be discussed i n 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 i n parameter space which then allows very rapid i d e n t i fication.  P r a c t i c a l o n - l i n e i d e n t i f i c a t i o n , however, requires the use  of d i g i t a l computers and the discrete version of parameter tracking introduces several d i f f i c u l t i e s .  9  The optimum gain or step-size depends on the input signals and i s d i f f i c u l t to determine on-line [5].  Furthermore,  a steepest  descent search f o r a minimum does not generally r e s u l t i n rapid i d e n t i fication. A popular method f o r on-line i d e n t i f i c a t i o n by d i g i t a l computers i s the mean-square technique.  In the l i n e a r system case a l l  state variables must be measured and the estimated parameters are biased. Furthermore,  i f rapid i d e n t i f i c a t i o n i s attempted,  a r i s e with i l l - c o n d i t i o n e d matrices.  d i f f i c u l t i e s can  A d i g i t a l estimation method i s  developed here which can overcome these d i f f i c u l t i e s by use of augmented input controls and output states.  Due to the imposed structure of the  estimator the manipulation of high-order matrices, which i s generally the case i n the mean-square approach, i s avoided.  2.2  Problem Formulation Consider a completely observable and controllable continuous-  time l i n e a r system x = Ax + Bu + w  (2.1)  where x i s an n-state vector, u i s an m-control vector and w i s additive plant noise.  Let T be the sampling i n t e r v a l and l e t t , t ^ , t^,  a set of sampling instants.  Q  be  For the moment, to ease the n o t a t i o n a l  burden, l e t z^(k) = z^(t^) be the sampled measurement at t ^ where z ^ t ) = x(t) + v(t) and where v ( t ) i s additive measurement noise.  (2.2) (The estimator to be  developed does not require that a l l states be measured as implied by (2.2).  The general measurement requirements  are discussed i n Section (2.4)  Let Z  i+l  -  ( t )  f  o\^ '> d  U  £ l  ( t )  i \ ™ * >  =  +  W  = l \ ^  fc)  d  ^  <'> 2  3  The vectors z^, u^, (£ = 2), are defined as augmented states and augmented controls, respectively.  I t i s seen from (2.1), (2.2) and (2.3)  that z„ = Az„ + Bu„ + w. % I I I A  A  •  where u^ = u, w^ = t  = x = t  k  k  +  1  (2.4)  v - Av + w.  Integrating (2.4) over the i n t e r v a l  and l e t t i n g Z (k) = z ( t > , U (k) - u ( t ) , y i e l d s Z (k+1) = AZ^(k) + BU (k) + N (k) £  £  k  £  £  £  £  k  (2.5)  £  where A = I + TA, B = TB and  V  k )  " [k t ( < ) k+1  5  z  T  " £< z  £  t  ) )  +  B ( u  k  Jl  ( T )  " £ u  ( t  k  } )  +  V  T ) ] d T  (2  '  6)  The discrete set of equations forms the basis f o r a l l further discussion. To ease the notational burden a l l n(n+m) elements of A and B are considered unknown.  The problem i s 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  z (k) 4  A  V - -  (2...Z)  c =  k)  b' n where prime denotes transpose and where a^ and b^ are the j - t h 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 + N (k) £  £  4  (2.8)  where  Yj(k) 0  0  ( Y  l  ( k ) )  (2.9)  Bn  •Y k) i(  (Due to the frequent occurrence of matrices of block-diagonal form i t i s convenient to introduce a s p e c i a l notation.  The f i r s t subscript B i s  used to indicate that a block diagonal matrix i s 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) i s an nx(n  + nm) matrix).  JO  Consider the sampling instants k-N+1,  k, so that N blocks of data  are available and l e t N(k-l)  Z(k) Z £ ^N  T  (2.10)  c + N(k-N)  Z (k-N+1) J where  z (k)  e(k-i)  x  K  N  , 6(k) =  •  mm  "  (2.11)  , Z(k)  e(k-N) n+m  n+m< >  Z  k  The matrix 9(k) i s n(n+m) x n(n+m) and the vector Z(k) i s n(n+m) x 1. The matrix 6^ i s of dimension nN(n+m) x n(n+m) as can be seen from (2.11). It i s 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 i s to determine an estimate of the parameter  vector c given 9^ and Z^, where these quantities are related by (2.10).  The simplest approach i s to define the optimum estimate to be the estimate which minimizes the quadratic form (Z^ - 8 ^ 0 ) ' ( Z ^ - Q^c). The solution i s the minimum mean-square (also c a l l e d least-squares) estimate  [6] c = P 6*Z. N N N M  (2.12)  T  where  h  1  = <w  (2.13)  Since the basic structure of a l l l i n e a r mean-square estimators i s the same i t i s of i n t e r e s t to compare the estimator (2.12) with previously proposed estimators [5, 12-16]. approach.  Mayne has used the Kalman f i l t e r  This requires that (2.10) be interpreted as a measurement  equation with 6^ as a deterministic measurement matrix.  This i s the  case only i f there i s no measurement noise and i f a l l states are measured. A further and very s i g n i f i c a n t difference i s i n the d e f i n i t i o n of 0^ and N(k).  These i n v e s t i g a t o r s take 9 (k-l)  N (k-1) x  1  ,  0 00  N (k-n-m)  1  (References case. (see  [2, 13-15], furthermore,  (2.14)  N(k)  0 (k-n-m)  treat only the scalar measurement  Their r e s u l t s do not d i r e c t l y apply to the above case).  Since  (2.5)) E [ Z ( j + l ) N ^ ( j ) ] = E [ N ( j ) N j ( j ) ] •* 0 1  1  i t follows that the estimate i s biased. viously [13].  (2.15)  L  This has been pointed out pre-  In order to eliminate the bias, s u f f i c i e n t time must  elapse between samples.  This can make rapid i d e n t i f i c a t i o n d i f f i c u l t .  Even when the bias i s small enough to give acceptable i d e n t i f i c a t i o n , further d i f f i c u l t i e s can a r i s e with data i n the form [2.14].  I f a small  sampling  i n t e r v a l i s chosen to achieve rapid i d e n t i f i c a t i o n there may  not be a s i g n i f i c a n t difference between adjacent rows of 6(k) r e s u l t i n g i n an i l l - c o n d i t i o n e d matrix.  These d i f f i c u l t i e s tend to be avoided i f  data i n the form (2.10) i s 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.  I t follows  from this assumption and (2.5) that E[Z (k)N'(k)] = 0, (Z, 0  Consequently 6.. at sampling _L  (2.16)  time t, i s independent of N (k) and the estimate ic p  (2.12) i s unbiased when N = 1. noise i s an important  p = 1, .... n + m)  From (2.6) i t i s seen that quantization  component of N^(k) making an analytic i n v e s t i g a t i o n  into the dependency of the noise samples i n t r a c t a b l e . observations can, however, be made.  Some general  Successive integrations performed  on w (T) not only reduce the e f f e c t i v e amplitude  of the noise but also  reduces the c o r r e l a t i o n between noise samples N^(k), (£ ~ 1 , ..., n + m). The truth of this statement for a certain class of noise samples i s shown i n Section (2.7). Consequently, for a given k, the quantization noise terms of (2.6) for d i f f e r e n t I tend to be uncorrelated.  Con-  sequently, l i t t l e i f any d i f f i c u l t y i s to be expected with estimation bias.  A l l examples t r i e d v e r i f i e d this expectation.  The e s s e n t i a l  difference between the two estimators can be seen from (2.15) and (2.16). If £ = 1, sequential data processing must be used.  I f T i s chosen large  to overcome i l l - c o n d i t i o n i n g then quantization noise can become excessive. The combination  of heavy noise and estimation bias (2.15) can r e s u l t i n  f a i l u r e of the conventional estimator to give acceptable estimation.  These  d i f f i c u l t i e s tend to be avoided i f augmented states and augmented controls are chosen (£ = n + m) so that p a r a l l e l data processing can be used.  H e u r i s t i c engineering j u s t i f i c a t i o n s can also be given for the use of augmented states and controls i n l i n e a r system parameter i d e n t i f i cation.  I f combined parameter estimation and optimal control i s attempted  signal-to-noise problems can a r i s e [12],  Optimal control attempts to  keep state deviations from a zero reference state small which decreases the signal-to-noise r a t i o . signal-to-noise r a t i o .  Integrators are low-pass f i l t e r s which improve  The data used i s then less noisy.  Such a s i t u a -  t i o n has been reported by Bakke and he suggests using a hybrid i d e n t i f i cation method [16].  B a s i c a l l y h i s hybrid method i s to use data i n the  form (2.14) and then replace Y^, 6^ by X^  y  82*  Two d i f f e r e n t  are taken and a decision i s then, made on the best estimate.  estimates However, a  f u l l set of augmented states i s not considered, nor i s the data combined i n a mean-square sense.  2.3  Recursive and Block Data Processing  The estimator (2.12) operates i n a block data processing mode. A new estimate i s made a f t e r every N a d d i t i o n a l data blocks.  Other  processing modes can be used such as the recursive method suggested by K i s h i where old data blocks are deleted and new data blocks introduced.  -  Let,  N-R  9  -  .> '  N  9  N-R N-R  ZZ  z- X T= N  R  (2.17) zR  It follows from (2.13) that P  N' =N-R R R P  +  9  9  (2.18)  Suppose that a l l quantities i n (2.12) are known and i t i s desired to determine  ^ from these quantities  (This amount to determining the  best estimate when R data blocks are deleted).  Making use of ( 2 . 1 8 ) i t  follows from ( 2 . 1 2 ) that N-R N-R - ^  9  " W  Z  -N-R  "  ( 2  1 9 )  By p a r t i t i o n i n g ( 2 . 1 2 ) i n the manner indicated by ( 2 . 1 7 ) , i t i s seen that °N = W R  N N-R N-R  +  P  6  Z  ( 2  '  2 0 )  substituting ( 2 . 1 9 ) into ( 2 . 2 0 ) y i e l d s  C  VR"  ( I  "  V R V ^ N  V R V  "  (2  ' > 21  With the aid of the i d e n t i t i e s <  "N R R " N R =N R  X  P  <  6  6  )  ~ VRV"  X  lp  1  e  =  P  1  +  6  (  "  I  ~N R R  ( I  P  e  0  W i ^  ) _ l p  N R R e  e  ( 2  '  2 2 )  ( 2 . 2 1 ) can be reduced to the standard from for recursive estimation  C  VR  = *N N R +  P  9  ( I  " R N R e  P  9  ) _ 1  (  V  " R> Z  N  ( 2  '  2 3 )  applying the matrix inversion lemma to ( 2 . 1 8 ) y i e l d s N-R = N  P  P  +  Vi  ( I  "  W R ^ V N  (  2  '  2  4  )  Equation ( 2 . 2 3 ) gives the best estimate i f the l a s t R data blocks are deleted and ( 2 . 2 4 ) gives the associated P-matrix. If R new data blocks are sampled, a new Z and new values of c^ and P^ can be computed.  and 6  D  are available  The equations giving these  values are obtainable from ( 2 . 2 3 ) and ( 2 . 2 4 ) by interchanging N with N - R, replacing 6_ by K  -8_ and leaving 0' unchanged. K  is.  With these changes ( 2 . 2 3 )  and ( 2 . 2 4 ) gives the Kalman f i l t e r algorithm (since the measurement matrix i s a random variable the Kalman f i l t e r i s not optimal). proposed i n References  The estimators  [ 5 , 1 2 - 1 5 ] are based on a scalar measurement. The  bracketed quantity i n (2.24) i s then a scalar and no matrix inversion i s required.  However, many scalar samples are required to estimate the  parameters within a given confidence  interval.  This can r e s u l t i n slow  identification. The matrix to be inverted i n (2.24) i s of dimension Rn(n + x Rn(n + m).  m)  This i s a consequence of the use of a vector measurement  of augmented states.  In e f f e c t , rapid i d e n t i f i c a t i o n i s attempted by  p a r a l l e l processing of augmented state samples.  the  The p r i c e paid for  p a r a l l e l processing i s i n the high order matrix to be inverted i n the Kalman f i l t e r algorithm.  An a l t e r n a t i v e recursive approach i s therefore  desirable. It follows from (2.9), and n+m  (2.14) that  k-1  SL=1 J=k-N where n+m  k-1  Q = £  £  N  V  j  )  -  Y  i  ( j )  ( 2  '  2 6 )  1=1 j=k-N consequently P  N = %\n  '  (2  27)  Equation (2.27) gives the i n t e r e s t i n g r e s u l t that only an (n + m) x (n + m) matrix need be inverted to find V^. n  2  + nm parameters, an (n  2  + nm)  x (n  2  + nm)  (In general, to determine matrix has to be inverted).  The structure of (2.27) allows a recursive algorithm to be developed which gives  iV  I  n  11  terms of Q"''. Consider (2.26) and N  = Y„(k  I  - N), S, = Q_ - V.V' ' 1 N 11 T  S  n+m  = Q  N-1 1  define  S. = S. - V. V'., j J-1 J j' (2.28)  applying the matrix inversion lemma y i e l d s  s- = 1  (Ql)  - v^p-  - Q J + Q" v c i + v ^ r ^ o "  1  1  1  1  1  1  (2.29)  Note that the quantities i n brackets on the right hand side of (2.29) are  scalars.  for  finding  The set of equations (2.29) define a recursive algorithm given  Repeating the algorithm R times allows  to be expressed i n terms of  (see (2.27)).  P^^  This procedure replaces  (2.24) i n the data deletion step. Substituting (2.18) into the f i r s t i d e n t i t y i n (2.22) y i e l d s  P  W  N-R R = 6  " WP"'  1  <' > 2 30  which allows (2.23) to be expressed i n the form ^N-R " ^N  +  P  N-R R VN 6  (  "V  The data deletion i s performed using (2.29) and (2.31). version i s required.  . ' (2  31)  No matrix i n -  Insertion of new data i s performed using (2.29)  i n reverse and by changing (2.31) to the appropriate form.  2.4  Measurement Requirements  1  Due to the use of state variable f i l t e r s not a l l 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 i s represented by x x  2  x  = a x  3  x  3  =  X  3  + a x 2  2  + a  l X ] L  + ...  3  = x  2  (2.32)  To i d e n t i f y  (2.32) i t i s s u f f i c i e n t t o generate the complete s t a t e s of one  the augmented models.  The f i r s t  augmented model t h a t can be g e n e r a t e d  c o m p l e t e l y i s the t h i r d model (Jl = 3) .  C  Z  ^ o  3 =  ^ o  f  t  o  It  i  x  d  i  T  ^  2  0  x„ d r , d x  0  1  3  Z  t  1  1  d r dx„  1 1  2  dx, l  (2.33)  2  i s seen from (2.33) t h a t o n l y x^ need be measured.  augmented model (l > 3) can be g e n e r a t e d u s i n g  2.5  x  o  •^ x o i  2  1  f  o  x„ d r , d x  2  o  d x  T h i s i s seen by t a k i n g  Any a d d i t i o n a l  (2.3) and (2.33).  Applications  Example 1. C o n s i d e r the l i n e a r  system x.^ =  X  2  a - | X  + a x 2  2  that only  simplification.  X  x  X  4  x  6  _8 x  =  = 0, t ^ = 16., and take I = 4 (see  x  l "  3  •  5  x  Q  5» 8 x  =  (2.34)  i n some  augmented system i s o f the form  X  3  where x^ = x^, g  0.0  2  "2  = X  , (°)  two elements need be i d e n t i f i e d r e s u l t s  l  x  1.0  x  The f u l l  x  + u, x^(o)  l  X  where a^ = -2.0, a^ = -3.0, t The f a c t  1  x  x  5  X  7  7'  " "  a  l "  l "  U  U  +  2 (2.35)  U  3  _4_ U  Equation (2.35) i s d i s c r e t i z e d and expressed appropriate d e f i n i t i o n of Z^(k), ^ ( k ) & an  c  '  i n the form of (2.8) by =  ( p a  »  The estimation of c i s performed by a block processing mode with R = 1. ( i f no old data i s deleted, the estimation algorithm i s e s s e n t i a l l y the Kalman f i l t e r algorithm).  The following sample i n t e r v a l s  and c o n t r o l 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 + s i n 3t + s i n lOt i s taken.  Figures  (2.1) and (2.2) give the r e s u l t s with zero additive measurement noise (n^ = o = n ) . 2  However, d i s c r e t i z a t i o n of (2.35) introduces quantization  errors which can be considered as additive correlated quantization noise. It i s seen from the figures that i d e n t i f i c a t i o n i s rapid even for large quantization noise.  Figure (2.3) and (2.4) i l l u s t r a t e the r e s u l t s f o r  state dependent noise of the form x,(k) n.(k) = a. • ( - i ) ng., ( i = 1, 2) m  (2.36)  where ng^ i s a pseudo-random computer generated noise sequence which has a maximum amplitude made.  of n . The choice a. = a„ = 0.1, E(ng.) =0. i s m 1 2 l  Figures (2.5) and (2.6) i l l u s t r a t e the results for the choice  = 0.1, a  2  = 0.25, E (ng^) = 0.1, E (ng ) =0.14. 2  I t i s seen from  these r e s u l t s that the estimator i s i n s e n s i t i v e to correlated quantization noise, state dependent measurement noise and to bias i n the measure ment noise. Example 2. Consider the l i n e a r system (See Reference [3])  X  X  i 2  = --2.0,  where  = x  2  = a  l X l  a  2  =  ,  x3^(0) = 1.0  + a x + bjU ,  X ( 0 ) = 0.0  2  2  -1.0, b  2  1.0, and take Z = 4.  1  system i s of the form • X  2  X  4  X  6  x  l  x  2  x  3  x  4  " l"  X  5  x  6  _ 2_  X  7  X  8  =  X u  8,  U •  3  4  l  U  2  U  3  +  b  a  (2.37)  l  ^ 4 _ U  •  •  where x_ = x,, x  a  u  5  r  = x,, x_ = x .  6'  7  8 n  The following sampling i n t e r v a l s are used: (1) T = 0.1, (2)- T =• 1.0-. (1) N = 21, (2) N = 3.  The-cor-responding- values of N are:  The value of R i s taken as R = 1 .  (1)  The i n i t i a l  parameter estimate i s c = 0 and u(t) = 0.1 + s i n t + s i n 3 t . Figures (2.7), (2.8) and (2.9) give the results f o r quantization noise only.  Figures (2.10), (2.11) and (2.12) give the results f o r  quantization noise and additive state dependent noise where  «  1  = a  2  = 0.1, E ^ g ^ = 0 , i = 1, 2.  Comparing these time responses with those obtained i n Reference [3] f o r the noise-free case, i t i s evident that i d e n t i f i c a t i o n i s more rapid. The speed of response of parameter tracking systems seems to be adversly effected by the i n t e r a c t i o n between parameter generating subsystems which i s an inherent feature of the steepest descent search procedure used.  Such subsystems are not required f o r the method developed here.  Consequently, i n the noise free case, one-step i d e n t i f i c a t i o n i s achieved.  FIG. (2.2) EXAMPLE 1. ESTIMATION OF  WITH QUANTIZATION  NOISE  0.4  FIG. (2.4) EXAMPLE 1. ESTIMATION OF a WITH QUANTIZATION AND 2  UNBIASED MEASUREMENT  NOISE.  0.0 2  4  t 6  8  10  12  U  FIG. (2.61 EXAMPLE!. ESTIMATION OF a BIASED . MEASUREMENT NOISE^  2  16  IS  20  WITH QUANTIZATION AND  0.4  r  1  1  1  1  1  1  2  3  4  5  ®-\  N  \  V  5  0  7  9  70  '  FIG.(2.7) EXAMPLE 2. ESTIMATION OF a, WITH QUANTIZATION NOISE.  -i 1  1  2  1  3  1  4  1  5  1  6  i  1  7  FIG. (2.8) EXAMPLE 2. ESTIMATION OF a  8  2  1  9  1—  10  WITH QUANTIZATION NOISE.  1.6 1.2  or* 0.4 0.0  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 UNBIASED MEASUREMENT NOISE.  WITH QUANTIZATION AND  2.4 2.0 1-6 1.2 0.8  Or**  0.4 0.0  w  FIG. (2.12) EXAMPLE 2. ESTIMATION OF b, WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE.  27  2.6  Comparison with a Limited Memory F i l t e r  In this section a comparison i s made between the proposed d i g i t a l i d e n t i f i e r using augmented states (I limited memory i d e n t i f i e r ( I = 1). fier  (or f i l t e r )  > 1) and the standard  The standard limited memory i d e n t i -  i s of the type discussed i n Section 2.3. The f i l t e r  operates i n a block data processing mode where R old data blocks are deleted and R new data blocks introduced r e s u l t i n g i n a limited memory f i l t e r .  (finite)  Example 2 i s used to i l l u s t r a t e the superiority of the  augmented state i d e n t i f i e r f o r various values of T and N (R = 1). Figures(2.13), (2.14) and (2.15) show the results for the i d e n t i f i c a t i o n of  the parameters a^, a^ and b^, respectively.  standard limited memory f i l t e r f a i l s  I t i s evident that the  to track the parameters f o r values  T = 1., N = 3, R = l and exhibit poor tracking accuracy for values T = 0.1, N = 21, R = 1. The f a i l u r e of the standard limited memory f i l t e r to give acceptable estimation of the parameters i s due to the fact that i t i s not an optimal f i l t e r f o r the proposed i d e n t i f i c a t i o n problem.  Consequently,  excessive quantization noise and c o r r e l a t i o n between samples, as i s the case here, can have a serious e f f e c t on i t s performance. In conclusion i t can be stated that the use of augmented states and controls r e s u l t s i n a more rapid and accurate i d e n t i f i c a t i o n of continuous time l i n e a r system parameters than can be achieved by either parameter tracking systems or the standard limited memory f i l t e r .  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 o WITH QUANTIZATION AND UNBIASED MEASUREMENT NOISE USING A STANDARD LIMITED MEMORY FILTER. 2  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 S t a t i s t i c a l Properties of the D i g i t a l I d e n t i f i e r  In the previous section, the performance of the d i g i t a l i d e n t i f i e r was  evaluated by comparing estimated values with known values.  a p r a c t i c a l s i t u a t i o n , however, the parameters are not known and means must be used to establish confidence i n the estimates.  In  statistical  In reference  [5] K i s h i made use of extensive studies made by Linnik [17] to determine confidence i n t e r v a l s f o r parameter estimates i n the case of a simple single-input single-output system.  The relations given by Linnik w i l l  be used i n this section to obtain approximate confidence l i m i t s for the parameter estimates of the d i g i t a l i d e n t i f i e r .  However, Linnik's results  can be applied only i f the d i f f e r e n t components of the noise vector N (see (2.10)) can be represented by a s t a t i s t i c a l noise model. task i s not an easy one.  Such a  However, i t w i l l be shown that under certain  simplifying assumptions i t can be argued that the d i f f e r e n t noise components w i l l be uncorrelated.  2.7.1  C o r r e l a t i o n of Noise Samples In Section (2.2), i t has been assumed that the noise samples of  the d i f f e r e n t models are uncorrelated. to be true.  Hence, equation (2.16) i s assumed  In the following, i t w i l l be shown that i f the noise samples  are uncorrelated then the samples of the integrated noise are also uncorrelated.  For s i m p l i c i t y , we w i l l consider a single noise component  in order to prove the v a l i d i t y  of the previous  statement.  Consider the noise component n(t) which can be represented by the following s t a t i s t i c a l model  E[n(t)] = 0 (2.38) E[n(t)n(x)] = a e  - 0 1  !^'  If a i s of i n f i n i t e value then n(t) w i l l be a type of white-noise component.  (In the white-noise case the amplitude a w i l l be equal to i n -  finity) .  For a p h y s i c a l r e a l i z a t i o n of the noise component n(t) i t  w i l l be assumed that a i s a very large but f i n i t e quantity and i s a very large but f i n i t e quantity and that a i s of f i n i t e  that a'  value.  Let N (t) x  I t follows from (2.38) and  =  n(T )dx 1  (2.39)  1  (2.39) that  E[N (t)n( )] = E[n(x)n(T )]dx 1 o 1 1 1  T  1  = a / o  To evaluate  (i)  e- l - l' a  T  dx, 1  T  the above i n t e g r a l we w i l l consider two  (2.40)  cases:  Case 1, x ^ t  E[N (t)n(x)] 1  = a £  =  (ii)  t  1  Case 2,  a a  [e  e~  a  (  T  -a(x-t)  _  T  l  _  dx^  )  e  -ax  o sc x «_ t  / \ i r E[N,(t)n(x)J = a J 1 o T-TUT  T  e  -a(x-x.) , 1 dr,  1  , ,t -a(x..-x) , + a f e 1 dx, x 1  a -ax ax , a ax -ax -at, = — e [e — 1J H e [e -e J a a r  a = — [ a r o  2  r  -ax - e  -a(t-x), - e ]  , (2.42) 0  32  For fixed time t and  a very large value of a, equations (2.41) and  (2.42)  can be approximated by  E  [  N  f=0 , x > t I (t)n(x)] J ' " ** U ^ , 0 < x < l  V  W  n  (2.43)  J  As seen from (2.43), the noise components n(t) and N^(t) correlated  become un-  as a -*- . 00  Using equation (2.41) or (2.42) for the s p e c i a l case where x = t yields E[N (t)n(t)] 1  = f [1 = —  The  following  e" ] at  ,  for fixed t and a ->• °°  non-dimensional r a t i o can be used as a measure  of the degree of c o r r e l a t i o n of N^(t)  E. (t)q(.t)] IMr  and  n(t):  r  = - r , for fixed t, a -> «  (2.44)  at Equation (2.44) shows that f o r fixed t and a large value of a the between N^(t)  correlation  and n(t) samples i s much less than the c o r r e l a t i o n between  the samples of n ( t ) .  This shows c l e a r l y the advantage of using the  grators i n order to generate the augmented states and  controls.  inte-  It i s  also i n t e r e s t i n g to study the degree of c o r r e l a t i o n between the samples of  N (t). x  Equation (2.39) y i e l d s E[N (t)N (x)] = 1  1  E[N (t)n(x )]dx 1  1  1  (2.45)  As before, we have two cases:  (i)  Case 1, x > t  ElN.WN.d)] = 1  1  f]2 - e a  o  a  l  T  -  e  -ax, .  «  [  2  T  I  l  1  a  ~  a  (  "  T  l  T  )  dx,  ]  1  -a(x-x,) . - £  a  a  o  a -ax , -at -a(x-t) = — j [2at + e +e - e a r o  =  (ii)  2at a  -1]  , f o r fixed t and a -*•  (2.46)  00  Case 2, o < x < t  E[N  l  (t)N (x)] = / i o  E[N (x)n(x )] dx, + / l l l x  T  E[N (x)n(x.)] dx, l 1 1  1  x a -ax -a(x-x ) , , = J — L2 - e 1 - e 1 ] dx, + o a 1 f  r o  1  1  t a -a(x,-x) -ax, -, , / — Le 1 - e 1] dx, x a 1 r  r  n ~ = —[2ax + e a a  a x  , ~ +e  a t  - e  -a(t-x)  -1]  2ax =  ,  For the case x = t equation  f o r fixed t and a •>  00  (2.46) or (2.47) y i e l d s  E[N- (t)] = ~ a  [at + e  2  = — a Equations  (2.47)  t  a t  - 1]  (2.48)  (2.48) and (2.38) y i e l d E[N (t)] —4 a l tE[n^(t)] 2  a  (2.49)  Equations (2.46), (2.47), (2.48) and  (2.49) show that the c o r r e l a t i o n  between the samples of the noise component N^(t) vanishes as  a -»• . 00  Equation (2.49) shows that f o r fixed t and large value of a the correl a t i o n between the samples of N^(t) i s much less than that between the samples of n ( t ) . The same analysis can be d i r e c t l y extended to ^ ( t ) where N (t) 2  =  N ( x ) dx 1  1  (2.50)  1  Since the same approach i s to be used i t i s s u f f i c i e n t to summarize the corresponding r e s u l t s E[N (t)N (t)] = ^ t 2  E [ N  (2.51)  2  1  2 2  at ( t ) ]  3  ~ fo~  ( 2  Efr-Cty^ct)], = t E[n (t)]  5 2 )  (2.53)  9  3  "  2  E[N (t)] 2  ~ t E[n (t)] 4  2  Equations (2.51), (2.52), (2.53) and same conclusions  (2.54)  1 3 a t  (2.54) show that the  that have been reached for N^(t) hold for N ( t ) . 2  Suc-  cessive integrations of s t a t i s t i c a l l y independent noise samples r e s u l t i n s t a t i s t i c a l l y independent noise components.  Also, successive i n t e -  grations of correlated noise samples r e s u l t i n less correlated noise components. The above s i m p l i f i e d mathematical treatment shows that the independence property assumed by equation (2.16) i s j u s t i f i a b l e for many p r a c t i c a l situations where n(t) i s 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 Z ( k + 1) = 6 ( k ) c + N(k) N  (2.55)  N  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 [ P ( k ) 9 ^ ( k ) Z ( k + 1)] N  N  = E [ P ( k ) e ^ ( k ) ( G ( k ) c + N(k))] N  N  = c + E[P (k)e^(k)N(k)]  (2.56)  N  The elements of P ( k ) and 8 (k) depend on Z (k) and U (k), A = 1, M  n+m.  Equation (2,55) shows that Z^(k) i s only related to N(k - 1).  Since N(k) i s s t a t i s t i c a l l y independent above analysis, then ^ ( k ) and 6jj(k)  a  r  of N(k - 1) as shown by the e  not s t a t i s t i c a l l y dependent on  N(k) and the following relations hold E[P (k)8^(k)N(k)] = 0  (2.57)  N  E[£] = c  (2.58)  Equation (2.58) shows that the estimate obtained by using the d i g i t a l i d e n t i f i e r i s unbiased.  The results given i n figures (2.1)-(2.4) and  (2.7)-(2.12) would not be obtainable i f (2.16) d i d not hold approximately. Equation (2.58) holds when 7^ = (9^9^) ^. to by K i s h i as the observable case.  This case has been referred  The case where the inverse (9^9^) ^  does not exist i s referred to as the non-observable  case.  In the non-  observable case an estimate i s obtained by use of the equation  c - 9* Z  where  (2.59)  N  i s the generalized inverse of 6 ^ .  That i s ,  e^e x = X  (2.60)  N  It follows from ( 2 . 1 6 )  2.7.2  and ( 2 . 6 0 )  that the estimator ( 2 . 5 9 )  i s unbiased.  The Error Covariance Matrix  Let V = 6a N  -  ^  (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 i n the form (see ( 2 . 1 2 ) )  =  V  Using equation ( 2 . 5 5 )  -  I ) Z  N  (  2  '  6  2  )  (  2  '  6  3  )  yields  ( E  For the non-observable  <Wi  N N N " P  6  I  )  S  case the expression f o r V i s v = (e e^ N  (2.64)  - I)N  where use i s made of the property 6*B = N N N  N  fl  (2.65)  The error convariance i s defined by Q  w  = E[W»]  (2.66)  Making use of equation (2.63) y i e l d s  Q  W  =  E  [  (  W N  "  I  )  M  '  (  N N N  9  P  "  E  I  )  ]  (  2  *  6  7  )  For the non-observable case using (2.64) y i e l d s  Q  2.7.3  W  =  E [ ( 6  " ^ ' ( V N  N N 6  "  I  )  ]  (  '  2  6  8  )  Estimation of the Confidence-Interval Length  Only the observable case w i l l be studied i n t h i s subsection since f o r the non-observable case the length of the confidence i n t e r v a l w i l l depend on the degree of redundancy i n the measurements which, would change from one system to the other.  =  V  N N N "  ( 6  P  6  I  )  Equation (2.63) y i e l d s  S  (  2  >  6  3  )  The fol*lowing relations' hold f o r the observable case:" R(8 ) = R(6^) = R(6^6 ) = R(P ) = n(n + m) N  N  (2.69)  N  where R(*) stands f o r the rank of the matrix written between the brackets. Let the matrix Y„ be defined as N  Y  N  - N N P  6  (  2  -  7  0  )  Using the known lemma R(AB) 4 min(R(A), R(B))  (2.71)  where A and B are any a r b i t r a r y matrices, i t follows from (2.70) that,  R  <V  =  n  (  n  +  m  N N  ) Y  N  )  (2.72)  Equation (2.70) y i e l d s 3  N  =  ( 9  9  (  2  '  7  3  )  Hence R(8jJ) = n(n + m) < min(R(6^9 ) , R(Y^)) N  < min(n(n.+m), R(Y ))  (2.74)  Hence R ( Y ) > n(n + m)  (2.75)  N  It follows from (2.72) and (2.75) that  R(Y ) N  Let the matrix  N  = R(P <3N) = < N  N  +  M  n  >  < ' > 2  76  be defined as  "HHW  -  (2  77)  Then W9 M  M  N N  It follows from (2.78) that  = 9..  (2.78)  N  i s a projection matrix.  Using (2.78) and  (2.71) y i e l d n(n + m) = R(8 ) N  < min(R(W ) , R ( e ^ ) N  < min(R(W ), n(n + m)) =  N  (2.79)  Hence. R(W ) N  ^ n(n + m)  (2.80)  Applying (2.71) to (2.77) i t i s seen that R(W ) N  4 n(n + m)  (2.81)  It follows from (2.80) and (2.81) that R(W ) N  = n(n + m)  (2.82)  Considering now the computation of the number of degrees of freedom of the  random vectors c and V, that i s (as defined by Linnik [17]) the  number of s t a t i s t i c a l l y independent components of c and V.  Starting by  computing the degrees of freedom of the random vector c we w i l l make use of Theorem 2.3.1 given by Linnik [17] which states that i f the vector U i s expressed as U = AX  (2.83)  and the matrix A i s of dimension r * £ with r < Z and R(A) = r then the random vector U i s of r degrees of freedom.  From equation (2.12)  £ = P 9'Z N N N M  (2  M  Equation (2.76) y i e l d s R(P. 0')  = n(n +  T  m)  (2.84)  Applying Theorem 2.3.1 given by Linnik [17] i t i s readly seen that the vector c i s a random vector of n(n + m) degrees of freedom. To compute the degrees of freedom of the noise error vector V, (2.63) i s expressed i n the form V = (W  M  where W  - I)N  (2.85)  i s given by (2.77). Let 'N  The matrix  (2.86)  i s idempotent as can be seen from the following  = W  N  - 2W  = I -  N  + I  N  (2.87)  Accordingly, there exists an orthogonal matrix 0^ such that (see reference  [17])  where the matrix  i s a diagonal one.  The matrix $^ has the property  that $  The matrix  N*N " V N  =  1  ( 2  '  8 9 )  i s 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:  d . = d.. 2  n  ii  i.e., d.. =0  or 1  n  (2.91)  The number of d ^ s equals to unity i s given by the rank of the matrix W^,  that i s , by n(n + m). From (2.89) i t follows that R(I - W ) = R ( ^ ( I N  = Rd -  - W )$ ) N  N  y  = n(n + m)(N - 1)  (2.92)  The f i n a l r e s u l t follows from the fact that the t o t a l number of diagonal elements of  i s n(n + m)N and of these n(.n + m) are equal to unity  (see (2.91)).  Applying the above r e s u l t given by (2.92) to (2.85) and  making use of Theorem 2.3.1 [17], i t follows that the noise error vector V i s a random vector of n(n + m)(N - 1) degrees of freedom.  From the r e s u l t reached i n the previous statement i t i s seen that V'V can be expressed i n the form n(n+m)(N-1) V'V  =  I i=l  (2.93)  Let Q  N  = E[c c']  (2.94)  cc With c i s computed according to (2.12). The following r e s u l t s given by K i s h i [5] and Linnik [17] are If £ and E,^ are s t a t i s t i c a l l y independent Gaussian random 2 2 — variables with d i s t r i b u t e d as % and having m degrees of freedom, required.  then the t - d i s t r i b u t i o n i s formed by the following r a t i o  t =  (2.95)  g  m Let  (see Appendix IV) c. -c.  K =  1  / (Q  where (Q matrix Q  (2.96)  1  CC  ).. IX  ).. i s the i - t h element of the main diagonal of the covariance cc 11 cc  and c. stands f o r the i - t h element of the vector c. 1  Also  l e t (see Appendix IV) m = n(n + m)(N - 1) .2 £q = V'V  (2.97)  Then using (2.96), (2.97) into (2.95) y i e l d s c. - c. 1  "n(n  + m)(N - 1)  ^  ) > > ( v  cc  11  (2.98)  1  ,  v / n ( n  +  m  )  (  N  _  1 }  42 Using the t - d i s t r i b u t i o n  i t i s possible to compute the length of the  i n t e r v a l about the parameter c. which would include c. with a certain i  probability.  1  For example, l e t the assigned probability be 0.9,  then  ^ n C n + DCN - 1)1 4 Y > = 0.9 where P{*}  stands for the p r o b a b i l i t y  the brackets.  density of the expression between  The number y can be found from tabulated standard tables  using the number n(n + m)(N - 1) as an index.  Ic. 1  If  (2.99)  - c.I = / ( Q  l  Y  x  1  ).. ,  Hence  , sY V  M  cc i i n(n + m)(N -  1  (2.100)  TT  1)  the augmented states and controls are not used the denominator  under the square-root sign i n  (2.100)  w i l l be n(N - n - m).  Thus  using the augmented states and controls decreases the estimation error | c\ - c^| and r e s u l t s i n a better estimation accuracy. The range 2A of the confidence i n t e r v a l for a p r o b a b i l i t y 0.9  i s given by  of  2A = [c. ± YAQ ).. , . VZ TV l ' c c i i n(n + m)(N - 1)  (2.101)  %  2.8  Discussion For the observable case equation  the unknown parameters.  Equation  It  computed without  For the non-observable case the matrix  (2.59)  i s to be used instead of  order to estimate the unknown parameters. 6^ i s required.  i s used to estimate  The matrix P^ i s recursively  any need to invert a matrix. P^ does not e x i s t .  (2.12)  (2.12)  The generalized matrix  in  inverse  i s desirable to compute G^j i n a way s i m i l a r to that  developed i n this Chapter for computing P^ in the observable case.  This  requires a computation  algorithm to evaluate 0^ recursively without any  need f o r a matrix inversion. In the next chapter an algorithm w i l l be developed i n order to compute 6* r e c u r s i v e l y and no matrix inversion w i l l be needed.  This  algorithm i s computationally simple and i s suitable f o r real-time computations . A s i m p l i f i e d version of the algorithm that w i l l be developed i n the next chapter can be used i n order to compute P^. with another method to compute P^.  This presents us  Due to the s i m p l i c i t y of this algorithm  i t can be used f o r state estimation applications i n cases where the plant i s represented by a discrete model.  Since the generalized inverse has  numerous applications besides those which occur i n state and parameter estimation, the important topic of recursive computation  of a generalized  matrix inverse i s dealt with i n the next chapter i n a general way.  3.  RECURSIVE COMPUTATION OF  THE  GENERALIZED INVERSE OF A MATRIX  3.1  Introduction Because of i t s value i n the solution of minimum norm problems  i n f i n i t e dimensional  spaces, the generalized inverse (or pseudo-inverse)  finds important applications i n l i n e a r system theory.  Applications of  the generalized inverse to estimation and f i l t e r i n g problems are given i n [18], [6], [19]. [21].  Applications to control problems are given i n [20],  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] a l .  requires the formation of C = A*A, of A, and a series expansion of C.  where A* i s the conjugate  transpose  The method proposed by Pyle  [7]  requires a two-stage Gram-Schmidt orthogonalization process and subsequent solution of sets of l i n e a r equations. posed by Rust et. a l . [8].  A simpler algorithm has been pro-  The algorithm i s based on a two stage Gram-  Schmidt orthogonalization process.  Unlike the other two methods men-  tioned above, i t does not require the generation of powers of matrices or the solution of sets of l i n e a r 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], i s the unique n x m a t r i x A* which s a t i s f i e s the following r e l a t i o n s :  AA A  (3.1)  I  (3.2)  = A (AA ) * = AA  (3.3)  (A A)*  (3.4)  1  AA  1  X  Following Rust e t . a l . , suppose that A can be partitioned into the form A = (R,S)  (3,5)  where R i s an m x n^ matrix with l i n e a r l y independent columns and S i s an m x n R.  2  matrix whose columns are l i n e a r l y dependent on the columns of  The following relations h o l d :  RR = I  (3.6)  I  R  1  S  = (R*R) R*  (3.7)  -1  (3.8)  -RU  U = RS (I + UU*)  (3.9) 1  R  I  U*(I + u u * ) " ^ (In  (3.6)  and (3.10)  (3.10)  1  I i s a unit n^ x n^ matrix).  The generalized  inverse i s given i n p a r t i t i o n e d form by (3.10) where ^ by (3.7) and ( 3 . 9 ) , respectively.  and U are given  Rust et. a l . show that the inverses  i n (3.10) can be found by a two stage Gram-Schmidt orthogonalization process (the f i r s t stage generates R*)•  The p o s s i b i l i t y of replacing  both orthogonalization stages by a recursive algorithm i s considered i n 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) y i e l d s  A*AA  = A*  (3.11)  Let A =  (V,W)  (3.12)  be a p a r t i t i o n i n g of A where V i s composed of l i n e a r l y independent columns. Equations (3.7) and (3.11) can be applied to V:  v  1  (3.13)  = (V*V) V* -1  v*vv  =  (3.14)  V*  From (3.1) and the p a r t i t i o n i n g (3.12) i t i s seen that  AA V = V  (3.15)  Taking the conjugate transpose of (3.15) and using (3.3) (3.13) and (3.14) y i e l d s V ^  Let r_., (j = 1, columns of R.  1  = (V*V) V* = V -1  (3.16)  1  n^) denote the l i n e a r l y independent  Consider the p a r t i t i o n i n g k+1 (3.17) k+1  where R, = r , , R 1 1 n  i t follows that  1  = R, and where the b.'s are row vectors. J  From (3.17)  (3.18)  " W+i k+ik+i +r  Applying (3.6) to  b  yields  (3.19)  = I  Applying (3.11) to 8^+1  the p a r t i t i o n e d form (3.17) y i e l d s  1 1 1  (3.20) r* k+1  k+1  P a r t i t i o n i n g of (3.20) and solving f o r b^-fi gives  \ + i  -1  <WW  =  Replacing A by  r  k i +  (  I  "  Vk+i>  (3.21)  and V by R^ i n (3.16) gives  ^W^+l  (3.22)  "  Substituting (3.18) into (3.22) and using  (3.19) y i e l d s  w =V 1  WW'  (3.23)  - \v  (3.24)  Substituting (3.23) into (3.21) gives  W = wW  1  where °k+l =  ( r  k l k+l r  +  "  r  k l\< k l r  +  )  _  1  (3.25)  +  is a scalar. The recursive algorithm (3.17) (3.23) (3.24) and (3.25). ,1 •i the matrix product R^^k.  n  f o r the computation of R  i s based on  I t i s desirable, however, to compute  (3.24) from a recursive formula.  Substituting  (3.23) and (3.24) into (3.18) y i e l d s  = Vk" v i +  ( i  - \^ k it i ) r  +  +  ( i  - v£  (3  -  26)  Since (see (3.3)) (3.27)  i t i s seen that the second matrix on the right-hand side of (3.26) i s the outer product of two vectors  Summary of the Algorithm for Finding R"*"  Initialize. \+i" b  B  k+l k+1  ( r  ][  r  \ l k l  =  r  +  (  " ^  ( I  "  ) R* -1  - 1  \J  "  I  +  = r , R* = (R*R  t+A^W  k+i k+i ~  vAi  <  Set k = 1, R  R )  V l W  = \*i v i +  (I  -A^wW 1  k+1 +1  k+1  3..  k. =  k+1  4.  If k < n  1  go to 2,  The above algorithm replaces the f i r s t Gram-Schmidt orthogonalization stage of Rust's algorithm, which, i n the interest of  accuracy,  reorthogonalizes each column a f t e r i t i s f i r s t orthogonalized.  An addi-  t i o n a l advantage of the above algorithm i s the manner i n which A i s p a r t i t i o n e d into the form (3.5).  In the Gram-Schmidt process, the de-  49  t e c t i o n of a zero vector a f t e r orthogonalization (but before normalization) indicates that the vector being orthogonalized i s l i n e a r l y dependent on the previous vectors.  This vector i s consequently moved  i n t o S and the next column vector of A i s chosen.  In the algorithm  presented here, l i n e a r dependency i s detected when CL ^  (see (3.25)) i s  K."t"-L  zero.  Since a ^ ^ i s required i n (3.24) no a d d i t i o n a l computations are +  required to test f o r l i n e a r dependency.  I f 0^+1 ^  s  z e r o  »  ^  fc  e  associated  vector i s moved into S and the next column vector of A i s chosen.  3.4  Recursive Evaluation of A  Consider the p a r t i t i o n i n g k+1 \+l. " » k+l?' ( R  S  (3.28)  \+l " k+1  where S. = s i s the f i r s t column vector of S and S = S. The generalized 1 1 n I inverse R can be obtained by the recursive algorithm given i n Section n  ?  z  (3.3).  Applying  (3.10) to A^, i t i s seen that  J  (3.29)  1  C. = U*B. J 11  (3.30)  U. = ( , u , U]L  2  , u.) = R S. .J J  (3.31)  where D. = (I + U.U*)  1  The matrix inversion lemma [9]  11  (3.32)  - \ '  \li  is  -  D  k \  +  i  (  I  "kVA'W^A  +  used to develop a recursive formula for B^.  (3.29) and using (3.33) B  1  (3.33)  Setting j = k, k + 1 i n  gives (3.34)  k i= " ViVWk+i^k ( I  +  where  e  k + 1  = .(i + u g B R u +1  k  k + 1  (3.35)  )"  The recursive algorithm for the computation of A* given based on (3.28)  (3.30)  (3.31)  is  (3.34) and (3.35).  • Summary of the Algorithm for finding A* 1.  Initialize.  Set k = 1, B . = D ' V = [I  -  ]  (1 + u*u )~ u u*]R 1  1  1  I  -1  V i • - k i\ \ it i k (I  e  R  +  3.  k = k + 1  4.  If k ^ n^, go to 2,  5.  C  n  = (R S I  2  n  )*B 2  n  +  2  A  3.5  ) B  +  1  = A nv. 1  Examples  The f i r s t two examples are taken from the l i t e r a t u r e .  This  allows a comparison to be made showing the r e l a t i v e s i m p l i c i t y of the proposed algorithm.  Example 1.  This example i s given i n [22]. I t i s required to compute  the generalized inverse of  A =  1 -1 0 0  0 1 -1 0  0 0 1 -1  -1 0 0 1  (3.36)  It i s seen that (3.36) has the form (3.5) where  0 1 -1 0  1 -1 0 0  R =  0 0 1 -1  (3.37)  s =  The algorithm of Section 3.3 i s used to f i n d R  Step 1.  Initialization.  R  l  1 -1 0 0  " l r  R* = | ( 1 , -1, 0, 0)  Step 2,  a  2  b  0  B,  R  2 2 R  =  3  = ±(1, 1, -2, 0)  |(2, -1, -1, 0)  1 3  2 -1 -1 0  -1 2 -1 0  -1 -1 2 0  0 0 0 0  Step 2. (repeated)  b„ = f (1, 1, 1, -3) 1 "4  f3 -1 -1 - l ) [2 2 - 2 2  - J  =  R  3 -1 -1 -1 2 2-2-2 1 1 1 - 3  = R =  Having found R*, the algorithm of Section (3.4) i s used to evaluate A*. Since S consists of a single column vector s^, only the i n i t i a l i z a t i o n stage i s required.  u  i  "  R  I  S  I  =  -1 -1 -1  Step 1.  •i-i  3-3-1 1 1 3 - 3 - 1 - 1 1 3-3  Step 5 C  =U  A l  1  A  Example 2.  - | ( 3, 1, -1, -3)  l Pl =  1 8  3 -3 -1 1 1 3 -3 -1 -1 1 3 -3 -3 -1 1 3  This example i s given i n [ 7 ] .  (3.38)  It i s required to compute the  generalized inverse of  1 0 1 1 0 1 - 1 0 1 1 0 1  (3.39)  53  It i s seen that (3.39) has the form (3.5) where  1 0 1  R =  0 1 1  ,  l" 0 1  1 -1 0  s =  (3.40)  The algorithm of Section (3.3) i s used to compute R*.  Step 1.  Initialization  R  l  "  , RJ =  l  r  \ (1, 0, 1)  Step 2. 0U  =  2 3  f (-1, 2, 1) B  2  R  1  =  =  f (2, -1, 1) 1 2  =  1 , 2 3 -l  -1 1 2 V  K  Having found R^, the algorithm of Section (3.4) i s used to evaluate A*,  Step 1.  Initialization  U =  D"  1  = I - u  (up  l  1  Step 2.  (  l+ u j u ^ u *  3 [0  1 IJ  u ) 2  3  v  l 2  J  B  2  = - i -v (  °  3  15 - l 5  3  )}  h  Step 5. C  2 = * 2 - T5 3 "o U  B  A  Example 3.  (  3 -1 4 3  2 " 15  0 3 5 4 -5 -1 0 3  (3.41)  This example i l l u s t r a t e s a case where A does not have the  form (3.5) i n i t i a l l y .  I t i s required to compute the generalized inverse  of A =  0 1 1 -1 1 0 0 -1 -1 0 0 0  0 0 0 1  (3.42)  Note that the t h i r d column of A i s the sum of the f i r s t two columns. Generally, the detection of l i n e a r dependency by simple observation i s not possible.  The algorithm of Section (3.3) i s applied to (3.42).  It i s seen from (3.36) that there i s no change i n the computation given i n Example 1 u n t i l k = 2.  At this stage the algorithm gives  hence the t h i r d column of A i s shifted i s chosen.  This gives  = 0,  into S and the l a s t column of A  = 1 and  b  = (0, 0, 0, 1)  3  1 2 3 1 (  3  -1 -1 0 1 -2 cr 2 - 1 - 1 0 1 1-2 0 0 0 0 3  The interchange  of the l a s t two columns of (3.42) i s accomplished  by a permutation matrix P given by  55  1 0 0 0 1 0 0 0 0 0 0 1  0 0 1 0  (3.43)  The algorithm of Section (3.4) i s now used to f i n d A*' where  A = AP, A = AP  (3.44)  The generalized inverse A''" i s then given by (see [8])  A  Step 1.  1  = PA  (3.45)  1  Initialization u^ = R s^ =  1 0 0  -1 1 0  0 -1 0  0 0 3  Step 5. U  A  1 1 B  1 3  1  =  I » °» (1  -1 1 0 0  _ 1  0 -1 0 -1  '  0 )  0 0 3 0  (3.46)  The generalized- inverse A-^ is- found- from (-3-. 45) and- is-  A  =  -1 1 0 0  0 -1 1 0  0 0 0 3  (3.47)  4.  TWO-STAGE ESTIMATION OF TIME-INVARIANT AND  TIME-VARYING PARAMETERS IN LINEAR SYSTEMS  4.1  Introduction Because of i t s importance i n system engineering,  considerable  e f f o r t has been devoted to the problem of estimating the parameters of systems whose dynamic behaviour can be described by d i f f e r e n t i a l Two  equations.  d i s t i n c t methods f o r estimating time-invariant parameters have been  brought to a high degree of refinement. i s based on a steepest-descent  The parameter tracking method  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 i s based  on optimum f i l t e r i n g .  In i t s simplest form an optimum f i l t e r minimizes  a mean-square error cost function subject to a prescribed dynamical constraint.  Detailed descriptions of optimum f i l t e r s and an  extensive  bibliography are given by Sage and Melsa [9]. The estimation of time-varying more d i f f i c u l t problem. often give acceptable  parameters i s generally a much  Parameter tracking systems and optimum f i l t e r s  estimates provided that the parameters vary  c i e n t l y slowly with respect to time.  suffi-  Lessing and Crane [10] have made  use of i n t e g r a l transforms to extend the c a p a b i l i t y of parameter tracking systems to track time-varying the i n t e g r a l transform  parameters.  P r a c t i c a l implementation of  approach requires the assumption that the parameters  can be represented by polynomials of known order.  The estimation of  time-varying parameters by optimum f i l t e r s requires e s s e n t i a l l y the same assumption. differential meters.  I t i s assumed that the parameters are solutions of  equations of known form with unknown time-invariant para-  The problem i s then reduced  to the estimation of time-invariant  parameters of a higher order augmented system.  I t i s the purpose of  this Chapter to develop a computationally simple and accurate method for estimating time-varying parameters which i s not r e s t r i c t e d by a p r i o r i assumptions about the dynamical behaviour of the parameters.  4.2  Problem Formulation  Consider the continuous-time  x  p  = Ax + Bu, p  l i n e a r system  x (t ) = x p o po  (4.1)  where x i s an n x 1 dimensional state v a r i a b l e vector and u i s an m x 1 P dimensional control vector which i s a known function of time. The system matrix A i s considered to be decomposable into a time-invariant component A^ and a time-varying component A^:  A =• A  f  + A  (-4.2)  v  The most general estimation problem associated with (4.1) i s to estimate x , A and B given an observed P  continuous-time  over an observation i n t e r v a l t  measurement vector defined  4 t <_ t ^ , and given complete s t a t i s t i c a l  information about measurement noise and system noise.  However, i f emphasis  i s placed on computationally simple methods then a less ambitious problem must be considered.  Suppose that B i s known and that the state x^ i s  58  measured during the observation i n t e r v a l .  The problem to be considered  i s the estimation of A given x^ and given no information about noise statistics.  The estimation i s performed sequentially i n two stages.  During the f i r s t stage an estimate of the time-invariant A^ i s made given an estimate of A^.  During the second stage a new estimate of A^  i s made using the estimate of A^ found during the f i r s t stage.  This  section deals with the second stage of i d e n t i f i c a t i o n . Let  z. = x - x, 1 p 1  x, = A.x. + Bu, 1 f 1  It follows from equations z., = 1  (4.3)  x.(t ) = x 1 o po  (4.4)  (4.1), (4.3) and (4.4) that  A.Z..  f 1  + A x vp  = A.z, .+ 6(x )a, t l p —  z,(t ) = 0 1 o  (4.5)  In equation (4.5), a i s a vector whose components are the time-varying elements of A . v The equation A x = 0(x )a v p p -  (4.6)  defines a matrix 0 whose elements are l i n e a r functions of x . Let P g = x  p  - x  x  h A g(t) - g ( t - T ) T  (  4  >  y  )  In (4.7), x^ i s defined by (4.4), x^ i s assumed known (through measurements)  and h. i s an approximation  to g, with the accuracy of the approximation  depending on the fixed-time increment  T.  Consider the following quadratic  cost index  J  where z^ and  =  jfl  f  [||g-^ l| Q + 2  o  1  1  llh-  (4.8)  ^||" Q ]dT 2  2  are defined by (4.5) where the unknown vector a. be r e -  placed by an estimate a_. The weighting matrices be p o s i t i v e d e f i n i t e .  and  are taken to  I t follows from (4.3)-(4.8) that l i m J = 0 proT->0  vided that the time-varying parameter vector c i i s known exactly.  Con-  sequently, (4.5) and (4.8) can be used to formulate a problem i n optimum estimation.  In the problem formulation the unknown time-varying para-  meter vector i s considered to be a control vector f o r the system described 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 cons t r a i n t given by (4.5) i s a standard problem i n optimal control. Hamiltonian  The  f o r t h i s problem i s (prime denotes transposition).  H = p ' C A ^ + 6(x ) a ) - - | ( g - z ) ' Q ( g - z ) - j ( h - A z -Q (x ) a) ' Q (h-A z.j-6 (x ) a) 1  1  1  f  ±  2  f  (4.9) The costate equation i s  p = -H  z  = -A^p + Q ( g - z ) - A^Q (h-A 1  and the gradient condition i s  1  2  fZ;L  - 6(x )a), p ( t ) = 0 p  f  (4.10)  60  H = 9' (x )[p+ Q,(h - A,z a p L r 1  G(x )a)] p —  = 0  (4.11)  Solving (4.11) for the optimum control a y i e l d s  a = (9»(x )Q9(x )) V(x )[p + Q(h 2  p  2  -  A z^ ]  (4.12)  f  Existence of an optimal control for the formulated problem can be viewed as a c o n t r o l l a b i l i t y  condition i n a function space S of  elements (g, h ) , where h i s either equal to g or an approximation to g. Taking (4.8)  to define the norm, then the system described by (4.5)  defined to be controllable i n S i f element (z^  is  a control ji exists such that the  - g, z^ - h) has minimum norm.  S a t i s f a c t i o n of this  con-  t r o l l a b i l i t y condition requires the existence of the matrix inverse  in  (4.12). Equations ( 4 . 5 ) , boundary-value problem.  (4.10) and (4.12) constitutes a two-point  Due to the l i n e a r i t y of the equations the  solution i s e a s i l y found using standard methods.  Let  L(x ) ^ (6'(x )Q6(x )) V(x ) p  2  p  p  Substituting (4.12) into (4.5)  I  C.  (4.13)  p  and (4.10) i t i s seen that  9(x )L + r  (4.14  -r' 4  where  C  ±  = (I  C  2  = -  Q  6(x )  l  LQ ) 2  A  f  + A'Q^  9(x )LQ h p  Q  1 §  2  - A^Q (I 2  0(x ) LQJh  (4.15)  61  and where I i s the unit matrix.  The general solution of (4.14) has  the  form  - *(t ,  t)  f  Q  p(t )  P(t )  c  f  +  (4.16) y  2  where  A y =  A  1  iKt,  (4.17)  T) r(x)dx  o  L2 J  fc  y  and where i|;(t, x)  f  t  = /  i s the state t r a n s i t i o n matrix for (4.14).  I t i s seen  from (4.16) that y can be found by solving (4.14) taking z ^ ( t ) = 0 = p ( t ) Q  for  the i n i t i a l values.  (4.14) taking -Kt  Q  The state t r a n s i t i o n matrix i s found by solving  , t ) = I, r = 0.  Substituting the boundary conditions into (4.16) and  parti-  tioning ^(t^,, t ) i n the form  *(t , f  t ) =  "*11  |*12  *21  j*22  Q  (4.18)  i t i s seen that  p(t ) Q  .-1 -*22  y  2  (4.19)  The i n i t i a l costate given-by (4.19) i n conjunction with (4.5), (4.10)' and t  Q  (4.12) defines the optimal estimate of a. over the observation i n t e r v a l <^ t ^ t^.  The computational  procedure w i l l 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 i n  state estimation (Sage [11]). in the problem formulations.  There i s , however, a fundamental difference In the state estimation formulation,  (4.5)  and  (4.8) are replaced by  x  = Ax + Bu, x ( t ) = x P P o po  P  ~  =  J  1  t  /  t  f  [  ll*p  - * M  1  2 q  P  +  ll*p  - % \ \  (4.20)  2  ^  <- > 4  21  The minimization of J subject to the dynamical constraint of (4.20) i s a problem i n optimum f i l t e r i n g . finding an approximation  x P  The problem can be considered as one of  to z so that the error (z - x ) has minimum P P P  norm subject to a cost penalty associated with the norm of the control effort.  The s o l u t i o n 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 R i c c a t i transformation.  E i t h e r approach results i n the Kalman f i l t e r  algorithm, f o r generating, the. optimal, estimate, x^ of the state. x the measurements z over an observation i n t e r v a l t < T < t. p o = =  p  given.  Comparing  the second terms i n the integrands of (4.8) and (4.21) i t i s seen that there i s a fundamental difference between the state and parameter e s t i mation problems.  This difference manifests i t s e l f a n a l y t i c a l l y i n the  c o n t r o l l a b i l i t y condition (the existence of L ( x ) , see (4.13)) which p  does not a r i s e i n the s o l u t i o n of the state estimation problem. 4.4  Estimation of a P e r i o d i c A . ,  v  Algorithm I I . S  Several questions a r i s e concerning algorithm I.  the estimation accuracy of  The f i r s t question concerns the e f f e c t of measurement noise  on the optimal estimate.  I t i s seen from (4.7) that h, since i t approxi-  mates g, may lose i t s s i g n i f i c a n c e i f the data representing g i s noisy. To some extent this loss of s i g n i f i c a n c e can be taken into account by  decreasing the absolute value of the elements of the weighting matrix Q (see (4.8)).  2  A second question arises concerning the s a t i s f a c t i o n of  the c o n t r o l l a b i l i t y condition.  Algorithm I must be modified i f p r a c t i c a l  answers to the above questions are to be found.  When no information i s  available about noise s t a t i s t i c s there i s no unique best method f o r f o r mulating an estimation problem, since the choice of a cost index, such as the quadratic cost index given by (4.8), i s an a r b i t r a r y one.  Con-  sequently, estimation algorithms can only be judged on the basis of trade-off between computational  complexity and estimation accuracy.  the above discussion i n mind, an algorithm i s developed  With  for estimating  time-varying p e r i o d i c parameters. Let  x  3  = z  - x  x  2  =-A x- ••+ e(-x^);a., f  (4.22)  2  x- (t ) = 0*  2  2  Q  * (4.2-3)  It i s seen from (4.5) and (4.20) that  X  Equation  3  =  A  f 3 X  +  9 ( x  2  +  X  3 -' )  X  3  ( t  o  )  =  (4.24)  0  (4.23) i s taken to represent the dynamical constraint  for the new problem formulation. the dynamics are noise-free.  I t i s seen from (4.1) and (4.23) that  A cost index of the form given by (4.8)  must be s p e c i f i e d to complete the problem formulation which i s based on the following assumptions: (2)  (1). L(x^) i s well-conditioned (see (4.13)).  an estimate I of a i s a v a i l a b l e and an estimate x^ of x^ can be found  by solving (4.24) using successive approximations.  x • x  3  0  Z  = A x f  3  + Q(x )a, 2  = A j L + 0(x.)a, J. I. 1 —  x  3(t ) =  That i s , x^ i s defined  0  0  x_(t ) = 0 I o  (4.25)  (3)  The cost index i s (see (4.8) and (4.22))  - H X  3  o  2 q  2  1 + M  H  "  X 3  "  X 2  H\]  D  (4.26)  T  Successive approximations converge i f the norm of a i s s u f f i c i e n t l y small.  This r e s t r i c t i o n i s imposed to achieve convergence i n the two-  stage estimation of A., and A . f v  (If A., i s known, then the f i r s t estimation f  stage i s unnecessary; (2.4) could then be used to define x^ and no norm constraints need be imposed on a ) . I t i s always possible to avoid  ill-  conditioning of L ( p ) by modifying the values of i t s time-varying elements. X  The modifications, however, cannot be a r b i t r a r y .  They must enter into  the problem formulations i n a mathematically consistent manner.  In the  present problem formulation i t i s seen from (4.22) and (4.26) that L(x^) i s replaced by- L ( X j - ) . .  Equations- (4. 23) , (4..25). and. (4..26). define-a--two-  point boundary-value problem f o r the optimum estimate. Let T be the length of a sampling i n t e r v a l and l e t t,. - t be a fixed estimation i n t e r v a l .  = NT  At the end t ^ of an estimation i n t e r v a l  a new measurement sequence i s taken over an observation i n t e r v a l t , < t < t + HT. f  That i s , the next estimation i n t e r v a l i s t + £T < t < o  (N + £)T (see F i g . (4.1)) .  •OBSERVATION INTERVAL FOR NEW MEASUREMENTS  0I  T  Ha  NT(UN)T  h  =H  ESTIMATION  INTERVAL  OPERATION PERIOD  flG.(4.l)  THE RECURSIVE INTERVALS.  OBSERVATION - ESTIMATION  In order to solve Eqn. to extrapolate a^(t) . i s to f i t a polynomial  (4.25) over this i n t e r v a l i t i s necessary  This can be done i n a v a r i e t y of ways. to sampled values j i ( t  o  One method  + kT) , k = 0,  N,  and  to use the polynomial approximation for extrapolation. ..Alternatively, , i f i t i s known that the length of the estimation i n t e r v a l NT i s greater than the smallest period of the p e r i o d i c components of a_, then the known p e r i o d i c i t y of a_ can be used f o r extrapolation.  The following i s a  summary of Algorithm I I : 1.  A measurement sequence g ( t  Q  + kT), h ( t  Q  a_(t + kT) over an estimation i n t e r v a l k = 0,  + kT) and an  estimate  N are assumed known.  Q  Take a new measurement sequence over an observation i n t e r v a l k = N, N + I and define a_ over this i n t e r v a l by extrapolation. 2. new  estimation i n t e r v a l for a new  (4.26)).  4.5  Solve the two-point boundary-value problem associated with the estimate  (Equations  (4.23)*, (4.25) and  The two steps of the algorithm are repeated r e c u r s i v e l y .  Estimation of a Periodic A^.  Algorithm I I I .  Algorithm II makes recursive use of measurement data and vaious estimates.  pre-  I t i s 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^ i s generated i n a smoothed form by integrating (4.25).  Conse-  quently, i t can be anticipated that algorithm II w i l l be superior to algorithm I i f measurement noise i s present.  There i s no guarantee,  however, that the c o n t r o l l a b i l i t y condition i s s a t i s f i e d . L(x^) may  be i l l - c o n d i t i o n e d .  As mentioned i n the previous section, one  method for overcoming i l l - c o n d i t i o n i n g i s to modify the elements of L by modifying  That i s ,  the problem formulation.  time-varying  Let z^ be an estimate of z^ and l e t  X  2  =  A  f 2 X  +  9  (  x  l  +  V - '  2  x  ( t  o  =  J  °  (4.27)  represent the dynamics associated with a t h i r d problem formulation.  The  cost index i s taken to have the quadratic form (see (4.8)).  J = \ / o  t f  t  [||g - *  \ \ \ + ||h - x | | Q ] d T 2  2  2  2  (4.28)  Comparing (4.27), (4.28) with (4.5), (4.8) i t i s seen that z^ should be chosen so that x  approximates z^.  2  (Ideally, i n the noise-free case  and when L(x ) i s well-conditioned, then the choice z, = z, can be made. P 1  This r e s u l t s i n x  2  1  = z^ and Algorithm III degenerates into Algorithm I ) .  I f L(x^) i s i l l - c o n d i t i o n e d Algorithm II w i l l f a i l .  This i l l - c o n d i t i o n i n g  can be overcome i n Algorithm I I I provided z^ i s chosen so that L(x^ + z^) i s well-conditioned.  A further constraint i s imposed on z^ by the require-  ment that the two-stage estimation procedure converges.  This requirement,  which was discussed i n the previous section, l e d to the use of successive approximations and the assumption that the norm of a_ must be s u f f i c i e n t l y small.  I t i s evident that z^ i s constrained by several c o n f l i c t i n g r e -  quirements.  Some sort. of. mathematically  tractable, approach must be taken,  which accounts f o r the various trade-offs associated with i l l - c o n d i t i o n i n g , dynamical approximations and two-stage convergence. Solving (4.5) by successive approximations using z^ = 0 as an i n i t i a l i z i n g function and an estimate a. gives x  2  - A x f  2  +  6( )I,  x (t ) = 0  XL  2  x„ = A.x„ + 9(x. + x ) a , 0  J  1  J  1  Z —  Q  x (t ) = 0 0  J  O  (4.29)  6;  for the next two approximating  functions.  z  1  = ax  2  Let z^ be defined by  + 3x  (4.30)  3  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 r e s u l t s i n z^ being defined by successive approximations given by (4.29).  Numerous algorithms f o r the solution of l i n e a r  algebraic equations are given i n [24]. convergence r e s u l t s when a = 2, 3 = -1.  Abramov's method for accelerated (The j u s t i f i c a t i o n f o r using  these algorithms i s e a s i l y seen by the d i s c r e t i z a t i o n and subsequent rearrangement of (4.5) into a l i n e a r algebraic equation). Algorithm I I I i s based on the solution of the two-point boundaryvalue problem defined by (4.27), (4.28), (4.29) and (4.30). parameters a. and p can. be. determined  by numerical. experimentation., or, they  can be chosen on the basis of successive approximations Abramov's method (a  4.6  The scalar  (a  = 0, 3 ='. 1) or  = 2, 3 = -1).  TwO-Stage Estimation of A^ and  A^.  During the f i r s t stage an estimate of A^ i s found given an estimate of A^.  Three d i f f e r e n t algorithms have been developed i n the  previous sections f o r the second stage during which an estimate of A^ i s found given an estimate of A^. Together  The two stages are performed sequentially.  they comprise a composite estimator f o r A.  For reasons of  s i m p l i c i t y , the discussion i s l i m i t e d to the development of a f i r s t algorithm which i s associated with second stage Algorithm I.  stage  The modifi-  cations required i f second stage Algorithm II or I I I i s used are s t r a i g h t forward.  Equation (4.4) defines the dynamics associated with the f i r s t stage of estimation.  I f x^ were available through measurements, (4.4)  could be d i s c r e t i z e d and  estimated by use of standard minimum mean-  square error estimation procedures.  However, x^ i s not d i r e c t l y available.  Furthermore, a two-stage estimation i s proposed and the i n t e r a c t i o n between stages must be accounted for i f a convergent sequence of e s t i mates i s to be obtained. Chapter 2 has suggested augmenting (4.4) by additional states and controls to reduce the e f f e c t of quantization noise which i s i n t r o duced by d i s c r e t i z a t i o n .  It i s shown that the use of augmented states  and controls results i n rapid convergence of the estimates.  The augmented  states and controls are defined by successive integrations of (4.4): y. = A y .  x  y  i  ( t  o  )  =  1  A y  i  + B u. i  fx  J  =  i I  0  i-l  d T  ft o  ,  i > 1  (4.31)  A y  >  u. = •/" u. ^dx, l o i - l t  y  l  = X  l  u, = u 1  (4.32)  For' composite estimation" of" A, successive integrations has' the b e n e f i c i a l -  e f f e c t of reducing the influence of the time-varying components of A^ on the augmented states.  The augmented system of equations (4.31) are d i s -  cretized and rearranged into a set of l i n e a r algebraic equations.  The  optimum estimate f o r A^ i f taken to be the solution with minimum meansquare error (For d e t a i l s see Chapter 2). In the case of composite estimation the above procedure must be s l i g h t l y modified since x^ i s not d i r e c t l y available.  Consequently,  x^ must be estimated from e x i s t i n g data.  From (4.3)  and (4.5)  i t is  seen that  XJL = x -  (4.33)  p  where  z-1= A fz, l+ 9(x p)a,—' z.(t 1 o) = 0  (4.34)  c  In  (4.34),  and a. are estimates of  A^  and a_ respectively.  The f i r s t  stage of estimation consists of a minimum mean-square error solution of the d i s c r e t i z e d equations  (4.31)  using  (4.33)  and  (4.34)  to define x^.  The new estimate of A^ i s used i n place of the old estimate provided that i t r e s u l t s i n improved prediction of the state of the system.  The following cost index can be used as a basis f o r making this  decision: J  In  (4.35),  change of  l  =  I t/ 'l*l' p f  x  " l X  "  S  iM Q d" 2  T  ( 4  3  i s a p o s i t i v e d e f i n i t e weighting matrix.  '  3 5 )  The rate of-  between estimation i n t e r v a l s i s useful i n making decisions  concerning the displacement length (£T) of the observation i n t e r v a l . displacement can be increased when the rate of decrease of  The  goes below  a preset threshold.  4.7  Examples  In a l l 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 i n the computer by storing the values at grid points and taking a piecewise constant approximation  70  a(t)=a(t +kT), kT < t<(k + 1)T). Q  In the case of algorithms II and I I I ,  extrapolation of _a was performed using a very simple constant approximations, _a(t^ + kT) = a_(t^) , k = 0,  , and the i n i t i a l i z i n g  functions f o r successive approximations were taken to be i d e n t i c a l l y zero.  4.7.1  Example 1. Consider a f i r s t order system x = -(1 + b sino)t)x + u  (4.36)  where the control u i s 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 a l l the examples i s taken to be 10 seconds.  The length of the estimation i n t e r v a l , as  seen from (4.37) i s NT = 2 seconds and the same has been considered for a l l the subsequent examples. of a(t) = b (4.1).  The maximum percentage error of the estimate  s i n wt for d i f f e r e n t b and w values are given i n Table  Table (4.2) shows the r e s u l t s when a(t) i s a square wave:  a(t)  f 0.5,. =\ L-0.5,  0 < L<  0..5.  0.5 4 t < 1.0  (4.38)  In the case of algorithms II and I I I , two estimation i n t e r v a l s were required to obtain the results given i n the tables. with successive approximations (a = 2, g = -1).  (a = 0, 6 = 1 )  Algorithm I I I was  and Abramov's method  tried  4.7.2  Example 2. Consider the second order system x^ = - ( 1 + b^sin io^t)x + u^ 2  x  where  2  = -(2 + t> s i n o) t)x 2  2  1  - x  2  + u  (4.39)  2  i s taken to be a step function of magnitude equal to 2 and  {  t ,  0 < t <_ 3  1 ,  t > 3  (4.40)  In addition to (4.37) the following values are used: b^ = 0.4, b  2  = 0.6, co^ = 3,  = 2.  x^(0) = 0 = x ( 0 ) , 2  Table (4.3) gives the maximum per-  centage estimation errors for the parameters a^(t) = b^ s i n co^t and a (t) = b 2  2  s i n o) t, respectively.  In the case of algorithms II and I I I ,  2  two estimation i n t e r v a l s were used.  4.7.3  Example 3. Consider the t h i r d order system x^ = (1 + b^ s i n to^t) x x  2  = -(2 + b  x  3  = -(3 + a ( t ) ) x  where the controls u^ and u  2  3  3  s i n ui^t) x  1  + u^  2  - 2x  3  - x  2  3  + u  3  4- u  2  (4.41)  are taken to be step functions of magnitudes  0.1 and 0.5, respectively and where  ft , t< 3 "  2  3  l l . « >  =  (4.42)  The parameter ^ ( t ) i s periodic and given by  1 , 0 < t < 0.25 a (t) U l , 0.25 < t <_0.5 3  (4.43)  In addition to (4.37) the following values are used: x ( 0 ) = 0, b 3  ±  = 0.3, b  = 0.5, o> = 3., OJ = 2.  2  1  and  4.7.4  Table (4.4) gives the  2  maximum percentage estimation e r r o r s .  x^(0) = ^ ( 0 ) =  In the case of algorithms I I  I I I , four estimation i n t e r v a l s were used.  Example 4.  The next two- examples i l l u s t r a t e the composite two-stage estimator.  In both examples the f i r s t - s t a g e of estimation was i n i t i a l i z e d  by taking a. = 0 as the i n i t i a l estimate f o r c i and only one set of augmented states and controls was used ( i = 1, 2, see (4.31) and (4.32)). Consider the second order system  x^ = (1 + b^ s i n to^t)x + u^ 2  x  = (b-- s i n o) t.)x - + 2 2 a  2  2  where the controls u^ and u  2  2  the following values are used: 2  3  = 1.  +  U  2  (4.-44)-  are taken to be step functions of magnitude  equal to 2. and 3., respectively.  co^ = 1, c o = 2, a  x  1  In (4.37), £ i s taken equal to two and  x^(0) =  =  ^» ^1  =  0.3, t> = -0.5, 2  After two estimation i n t e r v a l s the maximum per-  centage estimation error f o r a^ was 1.2% and a f t e r ten estimation i n t e r v a l s , the maximum percentage estimation errors f o r a  2  and a  3  were 1.8%  72  and 1.0%, r e s p e c t i v e l y .  4.7.5  Example 5.  Consider the second order system s i n ui^t)x^ +  x^ = (1 + X  2  =  a  2 l X  +  a  3 2 X  +  U  2  (4-45)  where b^ = 0.4, a^ = -2, a^ = - 1 , co^ = 2 and where a l l remaining data i s the same as i n Example 4.  A f t e r the second estimation i n t e r v a l , the  maximum percentage estimation errors f o r a^, a^,  and a^ were 0.9%, 0.4%  and 0.3%, r e s p e c t i v e l y .  4.7.6  Example 6.  In the previous examples the measurements were noise-free. Example 1 has been t r i e d with zero mean gaussian measurement noise with unity standard deviation and the random noise generator has been randomly initialized.  A peak r.m.s. s i g n a l to noise r a t i o of 10 was taken.  following values were chosen:  The  b = 0.4, co = 1, a = 2, (3 = -1. The maximum  percentage estimation errors f o r algorithms. I.,.. I I and. III. were 6%., 3.2% and 2%, respectively.  In the case of algorithms II and I I I , seven and  four estimation i n t e r v a l s were used, respectively. The examples show that the proposed  algorithms and the two-  stage estimation method can give accurate estimates of l i n e a r system parameters.  I f the measurements are noise free, Algorithm I gives the  most accurate estimates.  In the case of measurement noise, Algorithms  II and I I I give better estimates.  In the above examples, no problem  11  was encountered with i l l - c o n d i t i o n i n g i n any of the above cases.  75  b  1  W  /  0.5  1  1.0  0.6  %  3-0  10  %  2.0  0.7%  3.0  0.9  3.0  0.93%  0.75%  15  %  JT 777  °(=2,  3=-1  777  TABLE  TABLE  (4.3)  0.85% %  TABLE  1.0% 0.95%  (4.1)  TABLE  (4.4)  (4.2)  5.  OPTIMUM CONSTANT FEEDBACK GAINS FOR LINEAR SYSTEMS AND  THE PROBLEM OF  TRAJECTORY SENSITIVITY REDUCTION  5.1  Introduction  The optimum s o l u t i o n of the l i n e a r state regulator problem with quadratic cost has assumed a central p o s i t i o n i n the development of modern control theory.  The optimum gain, however, i s i n general  time-varying and this has severely r e s t r i c t e d p r a c t i c a l applications. P r a c t i c a l considerations generally specify constant or piece-wise stant feedback gains.  con-  Considerable e f f o r t has gone into the development  of design procedures f o r constant feedback gains and i n the design of s p e c i f i c optimal c o n t r o l l e r s [25-29]. There are, however, computational  as w e l l as p r a c t i c a l d i f f i -  c u l t i e s associated with optimum constant gains.  Specifying a constant  gain r e s u l t s i n a constrained optimization problem and the s o l u t i o n of the associated two point boundary value problem i s generally more d i f f i cult.  Furthermore, the optimum constant gain depends on the i n i t i a l  conditions.  One method that has been proposed to overcome these d i f f i -  c u l t i e s i s to introduce the trace of a cost matrix which then allows  ini-  t i a l conditions to be eliminated i n the problem formulation by a r b i t r a r i l y assuming a uniform d i s t r i b u t i o n of i n i t i a l states over a sphere [25-27]. From a design standpoint, however, i t i s important  to i n v e s t i -  gate the e f f e c t of i n i t i a l conditions on optimum feedback gains before computing average suboptimum gains.  Furthermore, i n many p r a c t i c a l s i t u -  ations, state t r a j e c t o r i e s are controlled, and the effect of i n i t i a l  conitions on optimum gains can become important.  This occurs, f o r example,  i n c e r t a i n schemes f o r combined estimation and c o n t r o l .  During a f i n i t e  observation i n t e r v a l , measurements of the state are taken and used to determine parameter estimates of a l i n e a r model.  The l i n e a r model i s  used to compute a control strategy which i s applied i n a subsequent control i n t e r v a l [11], [ 5 ] . systems.  This scheme seems p a r t i c u l a r l y a t t r a c t i v e f o r adaptive  I t s p r a c t i c a l implementation, however, requires an on-line  computational  algorithm f o r determining  optimal or suboptimal  Trajectory s e n s i t i v i t y minimization i s an important problem that has received considerable attention [30-37],  gains. practical  A widely  adopted method i s based on augmenting the state vector with the s e n s i t i v i t y vector and formulating a l i n e a r state regulator problem f o r the augmented state.  However, i f time-varying gains are permitted, the  problem, i s i l l - p o s e d [35, 38]. A well-posed, traj.ectory. s e n s i t i v i t y minimization problem requires that the gain be constrained to be time-invariant. In general, the s e n s i t i v i t y vector i s affected by i n i t i a l conditions.  Consequently, i f feedback i s to minimize trajectory s e n s i -  t i v i t y , the e f f e c t of i n i t i a l conditions on the optimum constant must be evaluated.  I t may be possible to choose a suboptimal  gain  average gain.  A u n i f i e d treatment of several problem formulations dealing with time-varying gains and constant gains i s presented.  Extensive" use  i s made of the known r e s u l t s f o r the l i n e a r state regulator problem to develop a simple algorithm which results i n rapid convergence.  An  approximate method i s developed to replace the nonlinear two point boundary value problem by an algebraic optimization problem.  This method  appears to be suitable f o r on-line computation of suboptimum constant feedback gains.  5.2  Problem Formulation  Consider a l i n e a r time-varying  x = Ax + Bu,  system  x(0) = X  (5.1)  q  where the state x and the control u are n x 1 and m x 1 vectors, respectively.  The optimal control which minimizes the quadratic  cost  index (prime denotes transposition) 1 f = ± / (x'Q-x + u'Ru)dt, o 2 o 1 t  J  where  (5.2)  i s p o s i t i v e semidefinite and R p o s i t i v e d e f i n i t e , i s given i n  feedback form by u = -R  The time-varying  gain matrix  •''B'K.jX.  (5.3)  i s found by solving a matrix R i c c a t i  d i f f e r e n t i a l 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 i n the matrix A i s considered to have a s i g n i f i c a n t e f f e c t on s e n s i t i v i t y .  Let z =  represent  the  dq  closed-loop  t r a j e c t o r y sensitivity vector.  A feedback structure of the  form u = -R B'(K x + K z)  (5.4)  -1  2  i s postulated, where  = K| and  = K^.  This s p e c i f i c structure i s  postulated because i t resembles the feedback structure of the  optimal  control r e s u l t i n g from solving the state l i n e a r regulator problem (see (5.3)).  Substituting (5.4) into (5.1) and neglecting the seconds '  order s e n s i t i v i t y function —j  9q  yields  z = 1^ x + (A - F K j z , 3q J.  z(0) = 0  (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 d i f f e r e n t i a l  y = (A - FK)y = Cy,  y(0) = y  equation  (5.6)  r  where  F  0  A =  -  A  K  , K = 3q  Equation  1  0  J  (5.6) represents  0  l  K  x  2 ' o y  L 2 K  K  (5.7)  =  3  (5.1) and (5.5) when  = K^.  This constraint  i s relaxed i n (5.6) i n order to obtain a u n i f i e d treatment of the l i n e a r and nonlinear  formulations.  The design objective i s to choose gain matrices  and  so  that (5.2) i s minimized subject to the a d d i t i o n a l constraint that the feedback (5.4) reduces, s e n s i t i v i t y vector z.  i n some sense, the closed-loop t r a j e c t o r y A mathematically  tractable approach that has  been widely adopted i n s e n s i t i v i t y studies i s to neglect second-order s e n s i t i v i t y terms so that (5.6) i s v a l i d and to introduce a cost index J(K  r  K, t ) = j /  f  f  (5.8)  y'(Q + KFK)y dt  where Q i s an augmented p o s i t i v e semidefinite weighting matrix  [32-38],  The optimum gain matrix minimizes (5.8) subject to the d i f f e r e n t i a l constraint (5.6).  I t w i l l be seen that the s e n s i t i v i t y problem formu-  lated by (5.6) and (5.8) i s well-posed  only i f the gain matrices are  constant. The argument l i s t i n the cost index (5.8) i s used to d i s t i n g u i s h  between d i f f e r e n t problem formulations. where  Let I(K^, K, t^) and I(K^, K, t ^ ) ,  i s a r b i t r a r y , represent the l i n e a r and nonlinear formulations,  respectively ( I f (5.6) i s rewritten i n the form y = Ay + Fu, then for the l i n e a r formulation the matrices A and F are completely known and independent of K^.  This r e s u l t s i n l i n e a r dynamic constraints.  For the  nonlinear formulation the parametric matrix A depends on the unknown gain matrix K^.  This r e s u l t s i n nonlinear dynamic c o n s t r a i n t s ) . The  Hamiltonian  for e i t h e r formulation i s  H = p'Cy  - |- y'(Q + KFK)y  (5.9)  where the costate vector p i s defined by p = -C»p + (Q + KFK)y,  Gradient matrices  p(t ) = p f  f  = 0  (5.10)  (see Appendix III) are used to derive the gradient  condition for a minimum.  In order to use gradient matrices a l l matrix  elements must be independent and cannot be constrained by the symmetry condition K =  K'.  This i s e a s i l y accomplished by taking K = K + K', where K i s a matrix with a r b i t r a r y elements.  In order to introduce a compact notation  which proves useful i n subsequent discussions, i t i s assumed that F i s time-invariant.  The gradient condition f o r the l i n e a r formulation  I(K^, K, t^) can then be expressed  8J(K  K,t ) i — — = 3K  i n the form (see Appendix V)  t / -_-5-dt=$  KF + FK$  J  o  3K  yy  +F$ yy  F = 0  +<S>  py  yp (5.11)  where  $  yv  k  o  y Jv  »  dt  (5.12)  81  (Equation (5.11) remains v a l i d i n the general case of time-varying F provided that KF and F are moved under the respective i n t e g r a l signs. Keeping this i n mind, i t i s seen that the form (5.11) i s a c t u a l l y not a r e s t r i c t i o n but rather a notational convenience). The gradient condition for the nonlinear formulation I(K^, K, t^) d i f f e r s from (5.11) i n an a d d i t i o n a l term which arises from the dependency of A on K^.  Let p' = (p|, p p  be a p a r t i t i o n i n g of the  costate vector which corresponds to the p a r t i t i o n i n g of y and l e t K^ = L + L', where L i s a r b i t r a r y .  The a d d i t i o n a l term which i s associated  with K^ i s  a  t  ~ f 3L o  f  Tr(-FK p!)dt 1 2 lZ r  (5.13)  = F$ + $ F p z zp„  Introducing the augmented matrix  F$ E =  + $  P z  Z p  2  2  F ' 0 !  0  1  (5.14)  0  allows the gradient condition for the nonlinear formulation to be expressed i n the form  8K  "  = E + $ KF + FKo + F'f + $ F = 0 yy yy py yp  Minimizing J(K^, K, t^) i s a constrained optimization problem.  (5.15)  Conse-  quently J(K*, K*, t ) = J ( K f  1 }  K, t ) f  where K* i s the optimum gain matrix f o r the nonlinear formulation.  (5.16)  If  K^ = K*, i t i s seen from (5.16) that the minimum for the l i n e a r formulation  82  i s achieved when K = K*.  The condition E* = 0, which i s s a t i s f i e d when F$  p z  + $  2  zp  F = 0  (5.17)  2  i s an optimality condition.  I t i s i n t e r e s t i n g to note that i f time-  varying gains are permitted,  the optimality condition (5.17) becomes a  s i n g u l a r i t y condition which i s impossible  to s a t i s f y [32, 37].  indicates that the s e n s i t i v i t y problem formulated  with  This  time-varying  gains i s i l l - p o s e d and cannot, i n general, meet the design objectives. Equations (5.11) and introducing a matrix S defined  FS$  + $ yy  (5.15) can be s i m p l i f i e d and u n i f i e d by by  SF = -F$ yy  - $ py  F - E  (5.18)  yp  (Several methods are available for solving an equation of the form (5.18) for S.  See Appendix I I ) .  Substituting (5.18) into (5.15) y i e l d s  9J(K K,t ) \ = $ (K - S)F + F(K - S)$ =0 9K yy yy  (5.19)  v  Equations (5.6), (5.10) and  (5.19) define a nonlinear  two-  point boundary value problem which must be solved i n order to determine the optimum gain matrix.  5.3  U n i f i e d Treatment of D i f f e r e n t Formulations  Because i t i s a constrained optimization problem, the computation of optimum constant  gain matrices  computation of unconstrained  i s generally more d i f f i c u l t than the  time-varying  gains.  It i s only i n the  s p e c i a l case of time-invariant systems and i n f i n i t e t  f  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 d i f f i c u l t i e s . evidently desirable to r e t a i n the comparative computational  It i s simplicity  of the unconstrained time-varying gain case, as f a r as this i s possible. This requires developing a u n i f i e d treatment  f o r the various formulations  and then making maximum use of the known results f o r the l i n e a r state regulator problem with time-varying gain and quadratic cost. Substituting the R i c c a t i transformation p = -Sy  (5.20)  into (5.10) and using (5.6) y i e l d s ' 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 l i n e a r formulation I(K^, K, t ^ ) , E = 0. Substituting  (5.20) into (5.18) y i e l d s t /  f  (F(S - S)yy' + yy'(S - S)F)dt = 0  (5.22)  Equation (5.22) i s s a t i s f i e d when S = S and (5.19) i s s a t i s f i e d when K - S =• S.  Substituting this- constraint into (5.21)- gives- the- standard  matrix R i c c a t i d i f f e r e n t i a l equation f o r the optimum time-varying gain. In the case of a constant gain matrix the above procedure must be modified.  H e u r i s t i c considerations indicate that some sort of time  average of S should r e s u l t i n 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 i s given by (5.22).  Replacing  Sy by p, (5.22) can be solved f o r S by use of Kronecker products (See Appendix I I I ) . The vector representation (S) s a t i s f i e s the l i n e a r  of the matrix S  equation  (5.23)  A method of successive approximations could be considered for solving the two-point boundary value problem. K = K^,  Using a nominal gain  (5.6) i s integrated i n the forward d i r e c t i o n and  backward d i r e c t i o n .  Equation  (5.16) i n the  (5.23) can then be solved for S and (5.19)  used to determine the gain increment.  The updated gain i s (5.24)  and the procedure i s repeated i t e r a t i v e l y .  While successive  approximations'  seems straightforward, they generally f a i l to converge i n the case of nonlinear two-point boundary value problems.  The s p e c i a l structure of  the formulated problem, however, allows some insight to be concerning  5.4  the convergence of successive  obtained  approximations.  Gradient Descent and Successive Approximations  Let (6K)  v  and  3J (—) be the vector representations for a matrix 3K v  gain increment and the gradient matrix, respectively.  The  incremental  cost associated with 6K i s given by  (5.25)  where  (5.26) i s p o s i t i v e semidefinite (see Appendix I I I ) .  Consequently, i f  6K = w(S - K), 0 < w = 1  (5.27)  then 5J = -w(S - K)'W(S - K)  (5.28)  = 0  The gain i s updated by  \+l " \  +W (  (5.29)  \"V  where the subscript k i s used to indicate an i t e r a t i o n stage.  For the  problem formulation under discussion, successive approximations (5.24) i s a s p e c i a l case (w = 1) of. gradient descent (5.29) with a weighting matrix (5.26).  Consequently, the proposed successive approximations  method i s equivalent to a deflected gradient method.  In the case of  quadratic cost indices i t is.known that some deflected gradient methods (such as conjugate gradients) r e s u l t i n s i g n i f i c a n t l y improved rates.  convergence  Having established the property that (5.24) i s a deflected  gradient descent, i t i s of i n t e r e s t to see whether i t has the properties of improved convergence, such as quadratic convergence.  I t i s known that  successive approximations r e s u l t s i n quadratic convergence i n the case of unconstrained time-varying gains.  Quadratic convergence can be proven  by considering (5.21) f o r two successive i t e r a t i o n s .  k+1  Introducing  S, , the l i n e a r matrix d i f f e r e n t i a l equation f o r AS, can be  solved using standard procedures, giving  -AS (0) o /  f c  k  *k+l >°  C  ( T  f  f  K  +  ±  (  T  ,  0  \+l ~  ) C (  )  - V ^ C t ,  ^ -  W  ~  +  (  \ " \+l  ) F  0)dx +  S  " \+l " (  k W '° ) ]  ( T  ) d T  (5.30) ' In (5.30) ^ ( t , x) i s the s t a t e - t r a n s i t i o n matrix of (5.6) at the k-th i t e r a t i o n stage.  I t i s e a s i l y shown that the cost index at stage k i s  given by J , = — y'S.y . k 2 o k o  The incremental cost A J , = J , , , - J , i s found k k+1 k  from (5.30) and can be expressed i n the form  -AJ  k  - T r t l y y C k K K ^ - \)H\  ~ VI  +±  + T r K ^ ) ^ - K^) ]  (5.31)  where  - C Wk+i \+i - V f  [  (  F + F(  \+i - V W k + i ^  (5  -  32)  and where the notation $ (k) i s used to indicate that (5.12) i s evaluated yy at stage k.  The gradient condition (5.19) can be obtained from (5.31)  by noting that Tr[(A + A')B] = 0 f o r a r b i t r a r y symmetric matrix B only when A + A' = 0 .  Equating (5.32) to zero y i e l d s (5.19).  The remarkable properties of the case with unconstrained timevarying gain are seen from (5.32).  By setting K ^ = S k+  k >  (5.32) vanishes  and the decrease i n cost (5.31) i s quadratic [38]. This i s accomplished by solving (5.21) .using successive approximations.  The same quadratic  convergence would r e s u l t i n the case of constant gain i f chosen to make (5.32) vanish. However, y^-^ depends on K^-^* condition ( ^ ) ^  =  could be  In p r i n c i p l e this may be possible. Consequently, any attempt to s a t i s f y the  0 would r e s u l t i n poor  convergence.  It i s reasonable to look f o r an approximation by replacing  87  y ^ ^ i n (5.32) by y^.  Let ( ^ r ) ^ represent the matrix a r i s i n g from this  +  substitution.  I t i s seen that  <M>k • V  k ) (  V-l  "  V  F+  f  (  \  +  l  " V yy $  where S^ i s defined by (5.22) evaluated at stage k. (5.33) that (^O^ = 0 i f (set w = 1 i n (5.29)).  ( k >  ( 5 i 3 3 )  I t follows from  i s determined by successive approximations  Since ( 7 ^ ) .  i s an approximation to ( 7 ^ ) ,  i t is  reasonable to anticipate that successive approximations r e s u l t s i n rapid convergence.  A rigorous proof of quadratic convergence would require  shoxtfing that the second term of (5.31) i s either never negative or, i f i t i s negative for a certain gain sequence, that i t s absolute value i s less than a fixed f r a c t i o n of the f i r s t term, which for Kj -^ f  1  c+  positive.  i s always  From a p r a c t i c a l point of view such a proof i s not e s s e n t i a l .  If the choice w = 1 should happen to r e s u l t i n divergence, then a choice w < 1 can be made which w i l l always give a convergent sequence (see (5.28)). The proposed computational algorithm i s the following. a nominal K = K^.  Integrate (5.6) i n the forward d i r e c t i o n and (5.10) i n  the backward d i r e c t i o n . use of (5.29).  Choose  Solve (5.18) for S^ and update the gain by  The rate of convergence i s controlled by w.  Computationally,  the proposed method i s s i m i l a r to the conventional steepest descent procedure.  The. essential, difference i s i n the use. of a variable weighting  matrix (5.26) to give a deflected gradient.  The associated change i n the  gain matrix i s actually given by the simple successive approximations formula (5.29). Limited computational experience to date indicates that this variable gradient d e f l e c t i o n allows a large step (w = 1) to be taken r e s u l t i n g i n rapid  convergence.  Note that the l i n e a r and nonlinear formulations of the s e n s i t i v i t y problem are treated i n a u n i f i e d manner.  This arises from the  fact that i n both formulations a nominal gain K i s chosen to perform the computations. The above algorithm i s s u i t a b l e for o f f - l i n e design studies. Since the optimum constant gain depends on i n i t i a l conditions, there i s a problem of p r a c t i c a l on-line implementation. i n i t i a l conditions and deriving suboptimal given i n . [ 3 9 ] .  A method for eliminating  feedback control laws i s  A method which i s s u i t a b l e f o r combined on-line e s t i -  mation and control i s discussed i n Section (5.6).  5.5  Optimum Constant Gain:  Trajectory S e n s i t i v i t y Reduction  Because of i t s engineering importance, open-loop and  closed-loop  system s e n s i t i v i t i e s have been extensively investigated (see Reference [32] for a bibliography).  However, for the closed-loop  s e n s i t i v i t y formulation given by (5.6) and unconstrained  (5.8), a paradox arises i f  time-varying gains are permitted.  This paradox has been  discussed i n the l i t e r a t u r e but i t s e f f e c t on s e n s i t i v i t y reduction does not seem to have been f u l l y appreciated [33, 37]. example can be- used to i l l u s t r a t e - the paradox.  x = -ax + u,  A simple  Consider-  x(0) = x  (5.34) o  z = -x - az - aKjZ,  z(0) = 0  u = -Hi K - K z = 0  (5.36)  2  J = -| /°°(x + Q z 2  2  2  (5.35)  + Ru )dt 2  (5.37)  and suppose that  and R are both chosen large compared to unity.  This  would seem a reasonable choice to make i n order to keep both s e n s i t i v i t y and control e f f o r t small.  The optimum time-varying gains, however,  s a t i s f y the constraint (5.36) and are i n f i n i t e . taking K  -> + » i  n  (5.35).  Equations 1 aK K  2  = K  This can be seen by  (5.34)-(5.37) then y i e l d .  x = 0  (5.38)  (5.39)  2  (5.40)  The optimum ( i n f i n i t e ) gain can r e s u l t i n a very s i g n i f i c a n t smaller cost than that associated with other suboptimal gains (see [37]).  (but p h y s i c a l l y realizable)  I t i s evident, however, that minimization of (5.37)  has not succeeded i n the design objective, which i s to decrease system sensitivity. open-loop.  The system (5.34) with the optimum control (5.36) operates The same d i f f i c u l t y i s experienced with a general n-th order  system (see (5.1), (5.4) and (5.5)).  The s e n s i t i v i t y problem defined by  (5.6) and (5.8) i s consequently i l l - p o s e d i f unconstrained time-varying gains are permitted.  From both a computational and an engineering point  of view the simplest constraint i s that of time-invariant gains.  It i s  i n t e r e s t i n g to note that this choice was made by Kreindler i n h i s design studies which dealt exclusively with l i n e a r time-invariant systems and t ^ = + <» [3].  I t seems natural to choose constant gains i n this case.  However, i f t ^ i s f i n i t e and/or the system i s time-varying, then the choice of constant gains i s not obvious.  90  5.6  On-Line Computation of Gains  From an engineering point of view a very a t t r a c t i v e f e a s i b l e on-line adaptive c o n t r o l scheme i s to i d e n t i f y a process as a l i n e a r time-invariant system of the form (5.1) over a f i n i t e observation i n t e r v a l and then apply a c o n t r o l of the form (5.3) over a subsequent f i n i t e control i n t e r v a l .  In such a combined estimation and control scheme i t  i s often important  to reduce the system s e n s i t i v i t y to incorrect parameter  estimates.  I t i s , however, not generally computationally f e a s i b l e to  solve a nonlinear two-point boundary value problem on-line. Some form of suboptimal control must be introduced.  It i s  known that the steady-state solution of the matrix R i c c a t i equation (take K = S i n (5.21)) often results i n a good suboptimal control.  With  this i n mind l e t P = -My  where M and N are constant matrices.  - Ny + e  (5.41)  Substituting (5.41) into (5.10)  and choosing M to be defined by MC + C M  = -Q - K F K ,  (5.42)  which i s the steady - state form of (5.21) gives  e + Ce  = -(NC  + CN)y,  e  f  = (M + N)y  (5.43)  f  A  For reasons of n o t a t i o n a l convenience the l i n e a r formulation I(K^, K, i s discussed.  t^)  Substituting (5.41) into (5.11) y i e l d s  F(S - M - N)$  yy  + $  yy  (S - M - N)F = -F$  ey  - $  ye  F  (5.44)  The decomposition (5.41) does not uniquely define N and e.  Consider the  p o s s i b i l i t y of defining N i n terms of observational data so that the r i g h t hand side of (5.44) has a n e g l i g i b l e e f f e c t on S.  The updated gain,  as given by (5.24), i s then S = K  2  = M + N  (5.45)  It i s seen from (5.43) that e i s "small" when N and e^ are "small". Besides the obvious trade-offs involved i n keeping N and e^  simultaneously  "small", there i s also the requirement f o r on-line computation using observational data.  The following cost index accounts f o r the various  factors which enter into an "optimal" decomposition (5.41):  Je = ^ / 2 o In (5.46) W^ and  f  (e + C'e)*W. (e + C'e)dt + ^ e/W.e. 1 2 r 2 f  are p o s i t i v e d e f i n i t e weighting matrices.  (5.46) Sustituting  (5.43) into (5.46) i t i s seen that Je i s an algebraic function of N, $  and y^y ^' Setting N = L + L', where L i s a r b i t r a r y , allows  matrices  1  to be used to determine the minimum. 1^oL  gradient  This y i e l d s  = TC + CT + W (M + N)i|) + iK(M + N)W_ = 0 2 t r 2 1  0  (5.47)  where T = W S$ + $ SW, 1 yy yy 1 1  S = NC + C'N <J> = y y ^ f  f  Equation (5.47) can be solved f o r N by use of Kronecker products. seen that  (5.48) It is  (S) (T)  v v  = G. (N) 1 v = G (S) 2  (TC* + CT)  v  [ G  1 2 3 G  G  +  4  G  ] ( N )  v  v  = G.(T) 3 v ~ 4  =  G  ( M )  v  ( 5  *  4 9 )  where G  1  = C  2 G  4  ,  © I + I © C ' =G  yy  1  = <|» ® W f  2  1  + W O ^ 2  3  yy (5.50)  f  E f f i c i e n t algorithms are available f o r evaluating (5.50) [40]. The solution of (5.49) f o r N involves one matrix inversion.  The system  parameters... in. the. C. matrix- are- determined, by the i d e n t i f i e r from, ob.servational data.  The matrices $ and i i , can be found by a forward yy f  integration of (5.6).  A l t e r n a t i v e l y , i t may be possible to use measured  states or state estimates to evaluate these matrices.  The gain compu-  tation (5.45) i s 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 i s computationally equivalent to solving these equations by successive approximations. optimum gain.  The f i r s t i t e r a t i o n must result i n an improved subI t i s desirable, but no e s s e n t i a l , that successive i t e r a -  tions r e s u l t i n a convergent Since  and W  2  suboptimum sequence.  are a r b i t r a r y p o s i t i v e d e f i n i t e weighting  matrices, i t i s not possible to reach general conclusions concerning the convergence of successive approximations.  I t i s easy to show, how-  ever, that given a nominal K, weighting matrices e x i s t which r e s u l t Suppose that the optimum gain K* i s determined  i n one-step convergence. by o f f - l i n e computations.  Equation  (5.42) can be solved for M.  (see (5.45)) N = K* - M and take W^ = I.  Equation  Let  (5.47) i s a Liapunov  type equation which can be solved for a p o s i t i v e d e f i n i t e l ^ .  With the  weighting matrices so defined there i s one-step convergence to the o  optimum gain.  It should therefore be possible to choose weighting,  matrices so that the f i r s t i t e r a t i o n r e s u l t s i n an improved suboptimum gain for small v a r i a t i o n s about nominal parameter values.  I f the system  i s time-invariant and i f N i s small i n the sense that terms involving products of the elements of N can be neglected, then i t i s seen by subs t i t u t i n g (5.45) into (5.42) that M i s 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 n o n - t r i v i a l i l l u s t r a t i o n of the successive approximations as w e l l as the approximation method discussed i n Section  (5.6)  t ^ should be chosen smaller than the s e t t l i n g time of the open-loop response for the nominal gain, and the nominal gain should d i f f e r s i g n i f i c a n t l y from the optimum gain.  These requirements are s a t i s f i e d i n the  case of the second order system  C  x  l 2  =  " 2 qX  = x  +  U  '  x  ,  x  l(°)  =  (5.51)  1 , 0  x (0) = 2  -1.0  t J  o  = \ f I o  f  (x  2  1  + 3x  2  + u )dt 2  2.  94  by taking q = 1, t f = 2.0, and  = (-10, 10). (The system i s an o s c i l l a t o r  for u = 0 and i s unstable f o r the chosen nominal gain). Successive approximations ((5.6), (5.10), (5.23) and (5.24)) required eight i t e r a t i o n s to converge.  No attempt was made to determine  weighting matrices which would r e s u l t i n rapid convergence for the on-line approximation method ((5.42), (5.49) and (5.45)). = W  2  = I was made.  An a r b i t r a r y choice  Six i t e r a t i o n s were required.  The r e s u l t s are  the following: 1.  Successive approximations K^** = (-1.06, - .215),  2.  J ** = 1.518  On-line approximation method. K ** = (-1.713, - .952), J * * = 1.565 q  (5.52)  The minimum cost for unconstrained time-varying gain i s J * = 1.503. It i s seen that both methods give very good suboptimal constant gains. The increases are only 1% and 4% for case (1) and (2), respectively. To i l l u s t r a t e the s e n s i t i v i t y problem consider (5.51) to be augmented by (5.5) and take J = 5 f o  tf  ( z + z )dt + J 1 2 o 2  (5.53)  2  A comparison i s made with Lee's et a l . method which uses a time-varying gain and the choice t ^ = 10 i s made.  The following results  are obtained  (a)  Lee's method [35] (augmented  state feedback; t = t )  u = -2x - 1.5x + .25z - .5z , J * = 1.94 ±  2  1  2  (5.54)  (b)  Constant gain (augmented state feedback) u = -6.8x  (c)  i  - 6.6x + 4.66 2  Z;L  + 20z  2>  J * = 1.70  (5.55)  Constant gain (Non-augmented state feedback) u = -1.812X  - 1.22x , J = 2.04  (5.56)  2  Figures (5.1)-(5.3) i l l u s t r a t e the improvements that result byuse of constant gain and augmented state feedback.  The magnitude of the  control i s smaller, the cost index i s smaller and there i s a s i g n i f i c a n t reduction i n trajectory s e n s i t i v i t y .  That the design objective of t r a -  jectory s e n s i t i v i t y reduction i s achieved can be seen by comparing the s e n s i t i v i t y functions of case (b) and case (c) (case (c) i s based on (5.51) taking t  f  = 10).  The optimality condition (5.17) was .0004).  checked and gave (-.0002,  I t should be noted that Lee's method uses a time-varying  gain  which i s a suboptimum solution of an i l l - p o s e d s e n s i t i v i t y problem. (The gain values given i n (5.54) are for t = t ).  It i s therefore not  surprising that an optimum constant gain can give s i g n i f i c a n t l y better results.  Note that increasing t ^ from 2 to 10 and using s e n s i t i v i t y  augmentation does not result i n an appreciable increase i n the cost index (compare (5.55) and  (5.52).  The convergence of successive approximation (5.29) was gated by  choosing K  ±  as the i n i t i a l gain. optimum  investi-  = (-1.5, -1.0), K  2  = (-0.1, -0.5)  This choice i s s i g n i f i c a n t l y different  gain (see (5.55)).  (5.57) from the  Only f i v e i t e r a t i o n s were required to converge  to the optimum. The following example i s used to i l l u s t r a t e the o n - l i n e approximation method when i t i s applied to a r e a l i s t i c system of moderate complexity. Figure (5.4) i l l u s t r a t e s a block diagram for the feedback control of a power generator connected to an i n f i n i t e bus [41]. The matrices i n (5.1) are  0.0 -0.676  0.0  0.0  0.0  25.0  0.0  0.0  1.0 -.25  A=  , 0.0  0.0  0.0  -0.06  (5.58)  B=  2.0  0.0  0.0 -2.0  2.0  -2.0  The components of the state vector represent the following v a r i a b l e s ; x^ = power angle deviation, x  2  = angular frequency deviation,  x^ = mechanical power deviation, x^ = governor p o s i t i o n v a r i a t i o n , u = governor p o s i t i o n c o n t r o l . 3.0,  The i n i t i a l state i s taken to be x^ = ( 0 . ,  2 . 0 , 0.0) and i n (5.2) and (5.46) the choice Q = W = W = I and ±  R = 1 i s made.  2  The choice t^ = 1 sec. i s made to accentuate convergence  problems, i n case they should e x i s t . The minimum cost using the optimum time-varying gain S(t) (5.21)) i s J * = 33.215.  (see  Successive approximations for the optimum con-  stant gain i s i n i t i a l i z e d 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  (b)  u* = -.069x  ±  (5.59)  - 1.22x -8.44x -3.35x 2  3  4  -1.47x -10.67x -5.47x 2  3  4  (5.60)  The suboptimal choice K = M i n case (a) amounts to taking the solution of the i n f i n i t e time case f o r the constant gain.  I t i s known  that (5.42) converges quadratically to the optimum gain f o r the i n f i n i t e time case.  The suboptimum choice K = M + N i n case (b) i s based on the  approximation technique discussed i n Section (5.6).  It i s seen to have  a quadratic-type convergence and r e s u l t s i n a smaller cost index than f o r case (a). In fact, the f i r s t i t e r a t i o n results i n a s i g n i f i c a n t decrease. A one-step improvement of this kind could be of p r a c t i c a l s i g n i f i c a n c e i n an on-line a p p l i c a t i o n .  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  0.0  to co-0.08 ^~  o UJ  £  -0.24  FIG. (5.2) EXAMPLE 1. TRAJECTORY SENSITIVITY FUNCTION 3 FOR(a) LEE'S METHOD- (b) CONSTANT GAIN AND AUGMENTED STATES, (c) CONSTANT GAIN AND NON - AUGMENTED STATES. 7  VO VO  FIG. (5.3) EXAMPLE I TRAJECTORY SENSITIVITY FUNCTION Z (a) LEE'S METHOD, (b) CONSTANT GAIN AND AUGMENTED STATES, (c) CONSTANT GAIN AND NON- AUGMENTED STATES. 2  F  0  R  :  SPEED GOVERNOR  FIG. (5.4)  TURBINE  EXAMPLE 2. BLOCK DIAGRAM VOLTAGE CONTROL SYSTEM.  GENERATOR  OF A SYNCHRONOUS  GENERATOR  6.  6.1  OPTIMAL ADAPTIVE CONTROL  Introduction  In this chapter we w i l l discuss one possible strategy for optimal adaptive control which i s p a r t i c u l a r l y a t t r a c t i v e point of view.  from an engineering  The strategy i s to i d e n t i f y system parameters at the end  of an observation i n t e r v a l and then to use the parameters to derive an "optimal" control for a subsequent control i n t e r v a l .  This chapter develops  one version of this strategy which i s based on the i d e n t i f i e r s and optimum controls developed i n Chapters 2, and  6.2  5.  Optimum Control During the F i r s t and Subsequent Optimization Intervals  The p a r t i c u l a r  optimization technique  to be used during the  f i r s t optimization i n t e r v a l i s of importance since the i n i t i a l  information  about the system i s usually i n s u f f i c i e n t to derive a good c o n t r o l . parameter estimates are poor then the optimum constant gains may s u f f i c i e n t l y close to the correct values.  If  not be  The problem of convergence can  become serious during the f i r s t optimization i n t e r v a l .  Hence, i t i s 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. operating at time t (6.1-a,  The system i s considered to s t a r t  and to terminate i t s operation at time t ^ .  d) (T) i s the length of the sampling i n t e r v a l .  The  In Figures length  of each observation i n t e r v a l i s (NT) seconds and the s t a r t of the i - t h  INITIALIZATION INTERVAL OPERATION  R  PERIOD  -sa-f«a  J  t -NT FIG. (6.1-a) 0  L  to THE INITIALIZATION INTERVAL AND THE OPERATION PERIOD. rT rT  t -NT 0  FIG. (6.1-b)  FlG.(6.1-c)  t tj° 0  NT  t°+NT  THE SUCCESSIVE INTERVALS.  OBSERVATION  THE SUCCESSIVE OPTIMIZATION INTERVALS (l = r).  IT  FIG. (6.1-d)  t  ,  T  c  THE SUCCESSIVE OPTIMIZATION INTERVALS (l = 2r).  observation i n t e r v a l i s denoted by the instant ( t ^ ) .  The displacement  between two successive observation i n t e r v a l s i s (rT) seconds (see F i g . (6.1-b)) where r i s an integer number.  The length of each optimization  i n t e r v a l i s (T ) seconds and the s t a r t of the i - t h optimization i n t e r v a l i s denoted by the instant (t^) (see Fig's (6.1-c, d ) ) . The displacement between two successive optimization i n t e r v a l s i s (&T) seconds. meter I i s a multiple of r , i . e . , I = mr.  The para-,  The value of m need not be  the  same f o r every two successive optimization i n t e r v a l s .  The choice of  the  parameter (m) i s based on the rate of v a r i a t i o n of certain perfor-  mance indices to be subsequently discussed (see subsection(6.3.4)). During the i n i t i a l i z a t i o n i n t e r v a l (t -NT < t < t ) s p e c i a l input test Q  signals can be used to drive the system.  Measurements of the input and  output at successive sampling instants are used to determine i n i t i a l estimates of the unknown parameters of the system.  The choice of the  i n i t i a l conditions at the s t a r t of the i n i t i a l i z a t i o n i n t e r v a l (t = t -NT) o i s completely a r b i t r a r y .  At the end of each observation i n t e r v a l new  estimates of the plant parameters are computed.  At the beginning of each  optimization i n t e r v a l new estimates of the "optimum" constant gains are made using the l a t e s t estimates of the plant unknown parameters.  These  new gains are used f o r the coming control i n t e r v a l of IT seconds.  6.3  Examples In this section a l i n e a r , continuous system with time-invariant  parameters i s considered.  Some of the system parameters are unknown  and an optimal adaptive control i s computed f o r the system.  The recursive  observation-optimization strategy described i n the previous section i s followed.  The i d e n t i f i c a t i o n of the unknown parameters i s performed using  the d i g i t a l i d e n t i f i e r developed i n 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 i n i t i a l i z a t i o n i n t e r v a l are taken to be step functions with unity magnitude, i n both examples. at  The i n i t i a l  conditions  the s t a r t of the i n i t i a l i z a t i o n i n t e r v a l are taken to be i d e n t i c a l l y  zero i n both Example 1 and  2.  It should be noted that use of a d i g i t a l i d e n t i f i e r r e s u l t s i n quantization noise.  Consequently, even though the examples treat deter-  m i n i s t i c cases, the system models are perturbed by quantization noise. This noise also a f f e c t s the value of the optimum constant which are computed using the i d e n t i f i e d parameters. noted that the optimum constant  It should be further  gains depend on the system state at the  s t a r t of an optimization i n t e r v a l . of" the proposed" adaptive"strategy  6.3.2  feedback gains  The examples i l l u s t r a t e convergence i n the presence of quantization noise.  Example 1 The system considered  i n this example i s represented  by  the  dynamic model l  X  x  2  =  a  l  x  2  +  U  =- a^x^  (6-. 1)  where the nominal values for a^ and a^ are -1.0 The i n i t i a l conditions are x^CO)  = ^(0)  of the recursive observation-optimization  = 1.0.  and 1.0,  respectively.  The e s s e n t i a l parameters  strategy (see Figs. (6.1-a,  d)) are taken to be: T = 0.2, T  t , = 10., f = NT = 2.0  t  = 0.,  r = I = 2, N = 10  o (6.2)  It i s required to compute the optimum constant feedback gains and  such that the performance index  2+t? J. = \ f l  (x + x + u )dt 1 L  1  2  • C  L  2  (6.3)  2  I i s minimized using the new estimates of a^ and a^ which are obtained at the end of the i - th observation i n t e r v a l .  The optimum constant gains  are computed using the successive approximations method developed i n Section (5.2) with w = 1.  The optimum solution  i s assumed to be obtained  when the inequality  ' i  " i  J  J  _  1  |  =  1 0  "  '  5  i s s a t i s f i e d , where J . i s the value of J . at the end of the k-th l  (  6  ,  4  )  iteration.  l  The time t ^ (see Fig.'s (6.1-c, d)) denotes the beginning of the i - t h , optimization i n t e r v a l .  During the f i r s t optimization i n t e r v a l the i n i t i a l  guess f o r the feedback gains  and  i s taken to be  = Y.^  =  -10.  For the subsequent optimization i n t e r v a l s , the i n i t i a l guess f o r the i - t h optimization i n t e r v a l i s taken to be the optimum values obtained f o r the ( i - l ) - t h optimization i n t e r v a l . after s i x i t e r a t i o n s  The stopping rule  (6.4) was s a t i s f i e d  f o r the f i r s t optimization i n t e r v a l .  No more than  two i t e r a t i o n s were required to s a t i s f y the stopping rule given by (6^4) for the subsequent optimization i n t e r v a l s .  The performance index of the  optimal adaptive system i s taken to be J* = •  a  25 E AJ* i=l  (6.5)  1  whe re the optimum value of J ^ (see (6.3)) i s denoted by J * and AJ* i s  i  i  the value of J * f o r an integration i n t e r v a l t ^ < t < t ^ + 0.4.  Since  the optimum constant gains computed at the s t a r t of the i - t h optimization i n t e r v a l are only used f o r a period of (JIT) seconds s t a r t i n g from the c instant t . , then J * as given by (6.5) i s the t o t a l operation cost. X  For  cL  the sake of comparison the optimum value of the performance index  J  o  = \ f I o  (x . + xl + u ) d t 1 l 2  (6.6)  2  i s computed using (6.1) with a^ and a^ being equal to t h e i r respective nominal values and the control feedback gains are taken to be the optimal time-varying gains.  The optimum value of J  q  i s found to be 2.046 while  the value of J * (see (6.5)) i s equal to 2.202, i . e . , the increase i n the performance index value does not exceed 8%.  Figures (6.2) and (6.3) show  that the i d e n t i f i c a t i o n errors do not exceed 5%.  6.3.3  Example 2 The system considered i n this example i s represented by the  dynamical model x^ = a^c^ + 2u x  2  = a x  3  + 2u  x  3  = a x  2  - u  2  3  where the nominal values f o r a^, a , and a 2  respectively.  3  (6.7)  are -1.0, 2.0, and -2.0,  The i n i t i a l conditions are x^(0) = x ( 0 ) = ^ ( 0 ) = 1.0. x  2  The e s s e n t i a l parameters of the recursive observation-optimization strategy are as given by (6.2). gainsK , K  2  and K  3  It i s required to compute the optimum constant  such that the performance index  108  2+t° J.=|J  1  (x  2 +  x  2 +  4  u )dt  (6.8)  2  +  i  i s minimized  during the i - t h optimization i n t e r v a l .  The optimum solution  i s assumed to be obtained when the inequality (6.4) i s s a t i s f i e d . In  this example, using the successive approximations method  with w = 1 (see equation (5.27)) during the f i r s t optimization i n t e r v a l r e s u l t s i n a diverging sequence of i t e r a t i o n s .  Thus w i s taken to be  and the i n i t i a l guess f o r the feedback gains are taken to be = -2, and  = -3.  rule (see (6.4)).  0.1  = -1,  Twenty i t e r a t i o n s are required to s a t i s f y the stopping  For the subsequent optimization i n t e r v a l s , the i n i t i a l  guess f o r the i - t h optimization i n t e r v a l i s taken to be the optimum values obtained f o r the ( i - l ) - t h  optimization i n t e r v a l .  No more than  two iterations, are required to satisfy. (6.4) for the subsequent optimization i n t e r v a l s and w i s taken equal to unity.  For the sake of comparison  the optimum value of the performance index  J  1 10 2 2 2 2 = i / (xf + x^ + x^ + u )dt o 2 o 1 2 3  i s computed using (6.7) with a^, a^,  (see (6.5)) was  The  found to be 0.655 while the value of J * a  equal to 0.705, that i s , the increase i n the performance  index value does not exceed 7.5%.  Figures (6.6), (6.7), and (6.8) show  the i d e n t i f i c a t i o n results for the unknown parameters a^, a^, and a^, respectively.  «  and a^ being equal to their respective  nominal values and the control feedback gains are time-varying. optimum value, J * . of J was ^ ' o o  (6.9)  1 U  The maximum i d e n t i f i c a t i o n errors do not exceed  5%.  The rapid rate of i d e n t i f i c a t i o n i s evident i n both this example and the previous one since a l l the unknown parameters have been i d e n t i f i e d  109  to a good degree of accuracy a f t e r only one observation i n t e r v a l . In Figures (6.4),  and (6.5),  the piece-wise constant curves  which are drawn as s o l i d l i n e s represent the optimum constant gain K* and K*, r e s p e c t i v e l y .  The parameters a^, and  (see (6.1)) are updated  at the end of each observation i n t e r v a l and a new estimate of the feedback gains i s made at the s t a r t of the following optimization i n t e r v a l . Since the optimum values of the constant feedback gains depend on the plant state at the s t a r t of the optimization i n t e r v a l ,  these optimum  values are not constant during the whole operation period (see Figure (6.1-a)).  The curves drawn as dashed l i n e s i n Figures (6.4),  represent the optimum time-varying feedback gains K^, and  and (6.5) 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 t h e i r respective nominal values. Figures (6.9), K*, K* and K*,  (6.10), and (6.11) show the optimum feedback gains  respectively.  the optimum constant gains.  The curves drawn as s o l i d l i n e s represent  The parameters a^, a^, and a^ (see (6.7)) are  updated at the end of each observation i n t e r v a l and a new estimate of the feedback gains i s made at the s t a r t of the following optimization i n t e r v a l . The optimum constant gains are not constant along the whole operation period but unlike the s i t u a t i o n i n Example 1 the feedback gains for Example 2 (see Figures (6.9),  (6.10),  (6.11)) change less  The curves draxm as dashed l i n e s i n Figures (6.9),  frequently.  (6.10), and (6.11)  represent the optimum time-varying feedback gains K^, K^, and These gains minimize the performance index given by equation (6.9)  respectively. subject  to the dynamical constraints given by (6.7) with a^, a^, and a^ being equal to t h e i r respective nominal values.  110  The recursive observation-optimization Section (6.2)  strategy described i n  i s characterized by i t s f l e x i b i l i t y .  due to the fact that i t  This f l e x i b i l i t y  i s possible to adjust the displacement distance  ( £ T ) between the successive optimization intervals c, d))  is  (see Figures  (6.1-b,  according to the behaviour of the system under consideration.  The  shortest time i n t e r v a l during which the gains remain constant i s equal to (2rT)  i n Example 1 (see Figures  Example 2 (see Figures ( 6 . 9 ) ,  (6.4),  (6.10),  4r for Examples 1 and 2, respectively, with less computational e f f o r t .  (6.5)) and equal to (4rT) (6.11)).  in  Hence, using I = 2r and  results i n a s a t i s f a c t o r y performance  The value of % used i n both Example 1  and 2 was % = r. In the following discussion two performance indices w i l l be discussed. not i t  These performance indices can be used to decide whether or  i s necessary to update the i d e n t i f i e d parameters and/or the con-  stant feedback gains at the end of each observation  6.3.4  interval.  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 i s a p o s i t i v e d e f i n i t e weighting matrix,  x  p  i s the measured plant state during the l a t e s t observation i n t e r v a l , and x  m  i s the corresponding model state during the same i n t e r v a l using the  l a t e s t i d e n t i f i e d 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 * ( t ..)x (t? C  j  In Eqn.  =  P  1  P  P  1  (6.11), x (t^) I  s  t n e  i - t h optimization i n t e r v a l (the present  X  -  X  P  ^  )  (6.11)  plant state at the s t a r t of the i n t e r v a l ) , and  x  (t^- -j) i s the  plant state at the s t a r t of the ( i - l ) - t h optimization i n t e r v a l . The assignment of values for the upper bounds of  and  J^y  respectively, would depend upon the degree of f a m i l i a r i t y with the p a r t i cular system under i n v e s t i g a t i o n (as can be seen from Figures (6.4), (6.5), (6.9), (6.10), and  (6.11)).  choice for the system considered  Taking IT = 2rT would be a reasonable  i n Example 1.  £T = 4rT would be a reasonable choice. be performed whenever  or  For Example 2 taking  Parameter or gain updating would  exceeded t h e i r respective bounds.  0-4. 0.2 0-0. 1  2  3  4  5  6  7  8  9  W t  0.2 'I 0.4 . -0.5. -0.8 . -10. FIG. (6.2) EXAMPLE 1. IDENTIFICATION OF PARAMETER QUANTIZATION NOISE.  o, WITH  1.0. 0.8. 0.6 . 0-4V 2 0.2 0.0 10 0.2 0.4 FIG.(6.3) EXAMPLE I. IDENTIFICATION OF PARAMETER QUANTIZATION NOISE.  a WITH 2  t  OPTIMAL SOLUTION ADAPTIVE SOLUTION  0.2  0.0 -0.2\  i  i  1 1 1  1  2  3  4  5  6  — 1— I 1 7  8  r9  -0-4 «1  -0.6 -0-0 -1.0 -1.2 -1.4  FIG. (6.4) EXAMPLE 1. THE OPTIMUM PIECEWISE CONSTANT FEEDBACK GAIN  -0.61  FIG.(5.5) EXAMPLE J. THE OPTIMUM GAIN K . 2  PIECEWISE CONSTANT FEEDBACK  -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 08 0-4  10 -0-4 -0-8  FI61S-7) EXAMPLE 2. IDENTIFICATION QUANTIZATION NOISE.  OF PARAMETER  aj WITH  as. 0-4.  ao. -0.4  "3  -as. -1-2 .  -IS -2X1.  FIGIS-S)  EXAMPLE 2. IDENTIFICATION QUANTIZATION NOISE.  OF PARAMETER  o, WITH  I  OPTIMAL ADAPTIVE  SOLTN . SOLTN  02  0.0  -0.1  -0.4  FIG. (5.9) EXAMPLE CAIN K  2. THE OPTIMUM  PIECEWISE  CONSTANT  FEEDBACK  t  OPTIMAL ADAPTIVE  0.4  SOLTN — SOLTN  —  —  0.2 0.0  I  2  3  4  S  S  7  S  9  -0.2 -0-4 -OS -08  -t-0 FIG.t6.IO) EXAMPLE GAIN Kj  2. THE OPTIMUM  PIECEWISE  CONSTANT  FEEDBACK  1.0. OS OS 0.4 3 0.2. OO -0.2  10 t OPTIMAL ADAPTIVE  SOLTN SOLTN  -  -  -0.4  FIG.(S-U) EXAMPLE 2. THE OPTIMUM GAIN Kj.  PIECEWISE  CONSTANT  FEEDBACK  7.  CONCLUSIONS  In Chapter 2 the generalized equation error has been shown to be applicable to a mean-square method of rapid d i g i t a l estimation of l i n e a r system parameters.  Due to the imposed structure of the estimator  the manipulation of high order matrices i s avoided. by the estimator have a block diagonal structure.  The matrices used 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 i l l u s t r a t e the  effectiveness of the estimator f o r a v a r i e t y of cases dealing with quantization noise as w e l l as measurements noise. In Chapter 3 a simple algorithm f o r computing the generalized inverse of a matrix has been developed.  The advantage of this algorithm  over other methods i s that i t eliminates the need for Gram-Schmidt orthogonalization and i t s associated ( i n the interest of accuracy) a l i z a t i o n as new vectors are introduced.  reorthogon-  Inherent i n the method i s a  very simple check on the l i n e a r dependence of the d i f f e r e n t matrix columns. In Chapter 4 a two-stage method f o r estimating, time-invariant and time-varying parameters i n l i n e a r systems has been developed.  During  the second stage, the time-varying parameters are considered as unknown control inputs to a l i n e a r subsystem of known dynamics.  This method f o r  estimating time-varying parameters i s computationally simple i n that the estimates are obtained by the solution of a l i n e a r two-point boundaryvalue problem.  Standard methods of solution r e s u l t i n one-step  convergence  A method of stage and control augmentation i n conjunction with d i s c r e t i -  zation and mean-square error minimization has been shown to be e f f e c t i v e for  the f i r s t stage of estimation.  This approach reduces undesirable  coupling between the two-stages which could r e s u l t i n convergence problems. Numerous examples i l l u s t r a t e the effectiveness of the method. In Chapter 5 method of successive approximations f o r solving the  two point boundary-value problem f o r optimum constant gain matrices  has been developed. The method i s shown to be computationally equivalent to a def l e c t e d gradient method.  Convergence can always be achieved by choice  of a scalar step-size parameter.  Rapid convergence i s important i n o f f -  l i n e design studies when the e f f e c t of various weighting matrices on optimal control performance i s investigated.  A close relationship be-  tween successive approximations and the solution technique for the l i n e a r time invariant case with i n f i n i t e t ^ i s exploited so that i n many cases rapid convergence i s achieved. The problem of trajectory s e n s i t i v i t y reduction by augmented state feedback i s shown to be well-posed i f gain matrices are constrained to be constant.  The problem of on-line implementation of a constant  gain matrix whose elements depend on i n i t i a l conditions i s discussed. The simplest approach i s to take a weighted average of the computed optimum gain matrices. feedback control law. back control  This r e s u l t s i n an average suboptimal l i n e a r A second approach i s to introduce nonlinear feed-  [39].  An on-line approximate method i s developed which appears suitable for systems whose parameters must be i d e n t i f i e d .  The two point  boundary-value problem i s replaced by an algebraic problem, the solution of which gives a suboptimal constant gain.  I t does not appear f e a s i b l e  to give a formal proof of convergence or a proof of improved  response.  By proper choice of weighting matrices, however, i t does seem possible to achieve an improvement with the f i r s t i t e r a t i o n .  Examples and para-  meter values were chosen to d e l i b e r a t e l y accentuate situations where convergence problems could a r i s e . countered and rapid convergence was  No d i f f i c u l t i e s , however, were enachieved i n a l l cases investigated.  A natural development of the present research i s to augment the s e n s i t i v i t y vector with the s e n s i t i v i t y functions with respect to the i n i t i a l conditions.  This w i l l reduce the s e n s i t i v i t y of the optimum  constant gains with respect to the state of the plant at the s t a r t of the optimization i n t e r v a l .  Accordingly, the optimum constant gains w i l l  be updated less frequently. Further research i s needed concerning the "best" choice of the constrained structure of the feedback gains.  One p o s s i b i l i t y has  been investigated i n this thesis where the gains are taken to be invariant. Let  time-  A more general choice i s to use s p e c i f i c time-varying gains.  the feedback gains be constrained by the s p e c i f i c structure K(t) =  FG(t), where the time-varying matrix G(t) i s to be s p e c i f i e d before s t a r t i n g the optimization process.  The time-invariant gain matrix F i s  to be computed according to the algorithm developed problem i s how  i n Chapter 5.  to choose G(t) systematically such that the r e s u l t i n g  time-varying gain matrix K(t) i s "optimum".  A s p e c i f i c time-varying  gain could result i n better r e s u l t s concerning s e n s i t i v i t y function reduction.  The  APPENDIX I INITIAL CONDITIONS LINEAR LEAST-SQUARES ESTIMATOR  The values of the optimum constant feedback gains as computed by the algorithm developed  i n Chapter 5 depend on the plant state at the  s t a r t of the optimization i n t e r v a l .  Accordingly, a good estimate of the  plant i n i t i a l conditions i s required.  In a noise free s i t u a t i o n , i t may  be possible to obtain the i n i t i a l conditions by measuring the state of the plant.  However, i f the measurements of the plant state are conta-  minated by noise then an i n i t i a l conditions estimator i s required. possible choice of such an estimator i s the l i n e a r least-squares  One  one.  This estimator does not require a p r i o r i knowledge about the noise s t a t i s t i c s and i s characterized by  i t s simplicity.  Considering-the- i-ttv optimization i n t e r v a l (see Figures- 6vl-a, b, c, d), the instant (t^) indicates the s t a r t of this i n t e r v a l as well as the completion of the i - t h observation i n t e r v a l . Let  the plant be represented  c c (for t ^ - NT < t < t^) by  x = Ax + B ±  where A^, and  (1.1)  u  ±  are the i d e n t i f i e d values of the plant parameteric  matrices A, and B, r e s p e c t i v e l y . These values are obtained according to the measurements taken during the i ~ t h observation i n t e r v a l . term on the right hand side of equation (1.1) B.u = -F.K 1 1 i  jX  can be expressed  The second as (1.2)  where the d e f i n i t i o n of the matrix F. i s as given i n Chapter 5 (F = BR "*"B')  and  represent the piecewise constant optimum feedback gain matrix  being used during the i - t h observation i n t e r v a l . Substituting (1.2) into (1.1) y i e l d s  x = (A. - F.K.)x i 1 1  (1.3)  The t r a n s i t i o n matrix i^(t^, x) of the dynamic system represented by  (1.3)  can be obtained by solving the following matrix d i f f e r e n t i a l equation  <L = -(A. - F .K. ) ' ii; r \ l l  i K t ? , t?) = I T  1  T  x  (1.4)  1  Let X j , represent the noisy measurements of the plant-state at  /s  C  t = t ^ - KT, and l e t x^ be an estimate of x^ defined by  *K  *K  =  x ( t  (1.5)  i>  where x(t^) i s the estimate of x(t^) which minimizes  the quadratic error  function e =  (x  -r  2  a  - x ) ' (x a a  - x ) a  (1.6)  where *K-*<^  The state vector x  ' i "  1  ^  X  a =  (1.7)  i s given by a  b  J  x  a  = ij, x(t?) a l  (1.8)  where (1.9)  N  The value of x(t^) which minimizes the error function e given by  (1.6)  i s known to be x(t ) = C  1  W a faJ ' V xa ,a  (1.10)  Equation (I.10) gives the best estimate of the plant state at t = t ^ which i s the i n i t i a l condition for the i - t h optimization (the present  interval).  to compute new  interval  These i n i t i a l conditions are required i n order  values for the optimum constant  feedback gains.  APPENDIX I I ANALYTICAL SOLUTION FOR THE LIAPUNOV MATRIX EQUATION [42]  This algorithm i s characterized by i t s s i m p l i c i t y since no matrix inversion i s needed and the solution i s obtained i n a closed form and not i n a numerical or i t e r a t i v e way. Consider the matrix equation S'X + XS = -Q where S, X and Q are a l l of dimension  (II.1)  (n x n)  Let A S, = k  k-i Z X.S i=0 k  (II.2)  1  where the parameters X ^ ( i = 0,  n-l) are the c o e f f i c i e n t s of the  c h r a c t e r i s t i c equation of the matrix S.  These parameters can be computed  e a s i l y using Faddeev's method.  X. = - ~ Tr[SS. . ] , S. = SS. . + X.I, i = 1, 1 x i - l' l i - l x ' ' X  Q  = 1  S  Q  = I,  S  '  n  = 0  (II.3)  The l a s t r e l a t i o n i s used as a check. We may write down the solution to equation (II.1) as  X =  1  2 n  "  2  j j - i E (-l) R(S, -S) k=0 n  K  k-l (S, -S) E s ; ( - l ) *S £=0 k  A  k  Q  1  %  k  %  (II.4)  123  where: R(S, -S) = ( - l ) 2 n  H (S) = n  n  T°- H (S) n  A,1 A o A3  0  A2  A  (II.5)  2  Ao  Ai 1  0  0 (II.6)  n  The determinant H (S) i s known as the Hurwiz determinant,  r  A (S, - S )  k odd  o  K+1  _ l» 2  where G^,  H (S) —7 G, (S) An k  n  f o r k even, i s the (- + l ) t h .  k even  (II.7)  Cofactor of the f i r s t column  of H (S). n Hence equation (II,. 4) can be written i n a s i m p l i f i e d form as  1  X = 2*^  n  n  (S)  _  1  E i=0  2 1  G (S)£,=0 E 2±  I 2i-£ (-iyS\QS  (II.8)  n4  From (II.8) i t i s clear that the condition for a unique solution i s H (S) ^ 0 which- i s equivalent, to. the condition, that, the dynamical system. (X = SX) i s asymptotically stable. For the n o n - t r i v i a l case where the matrix Q i s symmetrical, the term 2 1  E (-D A=O  I  S;QS  (II.9)  124 w i l l be symmetrical and since G„.(S) and H (S) are scalers, X i s also 2 i n J  symmetrical and we need only to solve n(n + l)/2 entries of X instead  2 of n . Thus f o r this case we have  X  U 1 1  where  =  -  r  T -  i  V G (S)Z  2 H (S) 1=0 N  L  <-!>*<SjQS  Si=0  %  i s the upper t r i a n g l e of X and ( pQ 2i-X?]J s1  S  l  i  s  ) X  %  t r i e  '  (11.10)  U  u  PP  e r  triangle  f o r (SlQS„. „•) f o r the case i = £, while f o r i ^ I i t i s a triangular matrix with i t s main diagonal equal to twice the main diagonal of (SlQS  ) and every element ( i , j ) i s equal to the sum of the corresponding  element ( i , j ) i n the upper t r i a n g l e of ( £ Q S  element i n the lower t r i a n g l e .  s 2 i  _ ) P £  l u s  the ( i + 1, j - 1)  125  APPENDIX I I I  The trace properties and gradient matrices used i n deriving the gradient conditions are the following [43],  Tr[AB] = Tr[BA],  Tr[AB ] = Tr[BA ] f  1  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 u t i l i z e d i n the algebraic manipulation of matrix equations [44], (AB)  v  = ( 1 0 A) ( B ) = (B' © I) ( A ) v  v  (A ® B) (C © D) = (AC) © (BD) (A + B) ® (C + D) = A ® C + A ® D + B ® C + B ® D  If Aj,, i = 1, .. ., n, and u , j = 1, of A and B, respectively, then  n, are the eigenvalues  u^. are the eigenvalues of A @ B.  Because of the p a r t i t i o n i n g (5.7) of F, can be found from the eigenvalues of F. equations:  (III. 2)  the eigenvalues of F  Consider the following eigenvalue  Fz. = X.z. x  $ V. yy J  x x =  (111,3)  U.V. i i  It i s seen from the defining equations f o r F and $  that (see (5.5)  and  (5.12)) (B'z.)'R  / o  1  (B'z.) = X z ' z. = 0 i  (III.4)  i  2 (v.'y) dt = u.v.'v. = 0  1  111  The matrix W i s therefore p o s i t i v e semi-definite.  Since ("j^)  =  W(S  - K)^  i t follows that the equality sign i n (5.28) holds only when the gradient condition i s s a t i s f i e d .  127  APPENDIX IV The confidence i n t e r v a l length deduced i n Subsection(2.7.3) (see  (2.101)) i s based on the assumptions that the random vector (c - c)  i s normally d i s t r i b u t e d and that the random vector V i s s t a t i s t i c a l l y r.  independent of (c - c) and i s normally d i s t r i b u t e d .  Also the random  2 variable V'V i s assumed to be distributed as % . The above mentioned assumptions can only be j u s t i f i e d i f the matrix 6^ (see (2.11)) i s deterministic.  In the following we w i l l prove  the above assumptions f o r the case where 6^ i s a deterministic  quantity:  Using (2.12) and (2.38) y i e l d s c - c = P 0'Z - c N NN M  M  = N N N P  6  ( 6  C +  N  )  "  C  = P 6^N  (IV.1)  N  where, f o r convenience, the random vector N i s taken to be normally d i s t r i b u t e d with E[N] = 0 and E[N N'] = 1 .  I t follows d i r e c t l y from  (IV.1) that (c - c) i s a normally distributed random vector. Using (2.63) and(EV.l) y i e l d s .  E[(c - c ) V ] = E[P 0^NN'(6 P 6^ - I ) ] N  =  P  N  N N N N N ~  =0  6  (9  P  e  N  I  }  (IV.2)  Thus the random vectors (c - c) and V are s t a t i s t i c a l l y independent.  The vector V i s normally distributed as can d i r e c t l y be seen  from (2.63) where  128  V  =  (0  N N N " P  I  )  S  =  (W  N  " ^  ( 2  '  6 3 )  Since 8JJ and P ^ are deterministic matrices then V i s normally distributed as N . Using (2.63), (2.64), and (2.65) y i e l d s  W  - D*(W - I)N  = N'(W  N  N  = N'(W  - I) N 2  N  = N ( I - W )N  (IV.3)  r  N  Considering the orthogonal matrix <J>^ (see (2.66)), equation (IV.3) can be written as V'V = N$ $'(I - W )$ $'N V N V N N V  V  V  N'$ (I - D ) ^ N N  N  = EyJ  (IV.4)  where y. i s a l i n e a r combination of a l l the elements of the random vector (the matrix (I - D^) i s a diagonal one with the main diagonal elements being equal to 1 or 0 (see (2.69))).  Thus the random vector y^  i s d i s t r i b u t e d as the vector $' N. The vector $'N i s distributed as N N the random vector N which i s a d i r e c t result of Fisher's Theorem [17]. Accordingly, the random vector 2 2 Ey^ i s d i s t r i b u t e d as % . 0^.  i s normally distributed l i k e N and  The above proofs are based on a deterministic  I f 0 ^ i s a random v a r i a b l e , i t i s no longer possible to prove that  the random vectors (c - c) and V are normally d i s t r i b u t e d along with the assumption  that V V i s d i s t r i b u t e d i s % . r  There appears to be no mathe-  matically tractable approach to derive confidence i n t e r v a l s i n the case  129  where 0^ contains quantization noise. Equation  To r e t a i n the basic s i m p l i c i t y of  (2.101) the following procedure i s suggested.  In equations  (IV.1 - 4 ) , i t i s possible to replace Z ( k ) by Z ( k ) where ^ ( k ) i s N  obtained by use of the best estimate (see  N  N  c at stage k, that i s by use of  (2.10)).  z  N  =? ;  dv.5)  N  The "smoothed" values 0^ of 0^ are used i n equations (IV.1-4).  Since  the randomness i n 6^ tends to be eliminated, i t seems j u s t i f i a b l e 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 ( I I I . 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'] - |- T r [ y y ' K ' F K ' ]  (V.2)  The following gradient matrix operations are required  g | Tr[FKyp'] = F(py') ^  Tr[FK'yp'] = (yp')F  " \ ~k  T r  tyy'  K f K  ] = " f ( y y ' ) K ' F - \ FK' (yy')  - \ g | Tr[yy'KFK'] = - ( y y ' ) K F -  T r  [ F K y y ' K ' ] = -FK(yy')  " \ a! [ y y ' ' T r  K  f K  ']  (  = - \ ^ ' ( y y ' ) - f(yy')K'F  The above operations allow ~  to be written as  v  .  3  )  13  3H  —  1  1  = - (py») - (yp')p - | ( « ) ' F - -| FK»(yy*) - (yy')KF F  yy  K  - 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)  where K = K + K  ?  From (V.4) and f o r time-invariant F and K equation (V.4) gives  t  a  H  _  t  £ - 9K = F / 0  t t t py'dt + / yp'dt F + / yy'dt K F + F K / yy'dt = 0 0 0 o (V.5)  Using the d e f i n i t i o n A  & = / yv  t J  yv'dt  O  Equation (V.5) can be written as t  gxr  _  _  _  / - -r^r = F<£> + $ F + $ KF + FK $ =0 9K py yp yy yy  q  This completes  the proof of equation (5.11).  (V.6)  132  REFERENCES  1.  P. V. Kokotovic', J . B. Cruz, J r . , J . E. H e l l e r and P. Sannuti, "Synthesis of Optimally Sensitive Systems", Proc. I.E.E.E., Vol. 56, No. 8, August, 1968.  2.  Lion, P., "Rapid I d e n t i f i c a t i o n of Linear and Nonlinear Systems", Joint Automatic Control Conference, August, 1966, pp. 605-615.  3.  Pazde ra, J . , and Pottinger, H., "Linear System I d e n t i f i c a t i o n Via Liapunov Design Techniques", Joint Automatic Control Conference, August,1969, pp. 795-801.  4.  Rose, R., and Lona, G., "An Approximate Steepest Descent Method f o r Parameter I d e n t i f i c a t i o n " , Joint Automatic Control Conference, August, 1969, pp. 483-487.  5.  K i s h i , F. H., "On-Line Computer Control Techniques and Their A p p l i cation to Re-entry Aerospace Vehicle Control", C. T. Leondes, ed. Advances i n Control Systems Theory arid Applications, V o l . 1, Academic Press, New York, 1964.  6.  C. P r i c e , "The Matrix Pseudo-Inverse and Minimal Variance 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 I d e n t i f i c a t i o n , Academic Press, 1971.  Estimates','  10.  H. C. Lessing and D. F. Crane, "The Use of Integral Transforms i n the Estimation of Time Variable Parameters", NASA Techn. Note, NASA TN D-5194, May 1969.  11.  Sage, A. P., Optimum Systems Control, P r e n t i c e - H a l l , Inc., 1968, pp. 507.  12.  Mayne, D. 0., "Optimal Non-Stationary Estimation of the Parameters of a Linear System with Gaussian Inputs", Journal of E l e c t r o n i c s and Control, V o l . 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 , " I d e n t i f i c a t i o n 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 i n the Paper-Making i n System Control, 1968, pp. 113-192.  17.  Y. V. Linnik, Method of Least Squares and the P r i n c i p l e s 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 E s t i mation i n a Dynamic System", JACC, 1971, pp. 820-829.  20.  D. Luenberger, Optimization by Vector Space Methods 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 General i z e d 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 Timevarying Systems", I.E.E.E. Trans., Automatic Control, V o l . AC-13, pp. 150-159, A p r i l 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, V o l . 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, V o l . 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, V o l . 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 S e n s i t i v i t y f o r a Large E l e c t r i c Booster", Proc. of the 5th Annual A l l e r t o n Conference, pp. 142-148, 1967.  Industry", Case Studies  9  John Wiley  134  31.  E. Kreindler, "On Minimization of the Trajectory S e n s i t i v i t y " , Intern. J . Control, Vol. 8, No. 1, pp. 89-96, 1968.  32.  E. Kreindler, "Synthesis of F l i g h t Control Systems Subject to Vehicle Parameters Variations", Grumman A i r c r a f t Engg. Corp., Bethpage, N. Y., Techn. Rept. AFFDL-TR-66-209, A p r i l 1967.  33.  E. Kreindler, "Formulation of the Minimum Trajectory S e n s i t i v i t y Problems", I.E.E.E. Trans., Automatic Control, Vol. AC-14, pp. 206-207, A p r i l 1969.  34.  J . F. Cassidy J r . and I. Lee, "On the Optimal Feedback Control of a Large Launch Vehicle to Reduce Trajectory S e n s i t i v i t y " , 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 J r . , "A Note on Trajectory S e n s i t i v i t y of Optimal Control Systems", I.E.E.E. Trans., Automatic Control, Vol. AC-13, pp. 111-112, February 1968.  37.  A. J . Bradt, " S e n s i t i v i t y Functions i n the Design of Optimal C o n t r o l l e r s " , I.E.E.E. Trans., Automatic Control, V o l . AC-13, pp. 110-111, February 1968.  38.  D. L. Kleinman, "On an I t e r a t i v e Technique for R i c c a t i 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 Feedback 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, V o l . AC-13, pp. 122-123, February 1968.  41.  I i . 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. L i n c o l i n 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., V o l . 17, No. 3, pp. 603-606, May 1969.  

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 11 25
United States 11 0
France 4 0
Sweden 3 0
Mexico 2 0
Russia 2 0
Canada 1 0
Brazil 1 0
Czech Republic 1 0
City Views Downloads
Unknown 10 10
Ashburn 5 0
Beijing 5 0
Harbin 3 0
Shenzhen 3 25
Stockholm 3 0
Mountain View 3 0
Salvador 1 0
London 1 0
Sunnyvale 1 0
Mladá Boleslav 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

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