IDENTIFICATION OF PARAMETERS IN DISTRIBUTED PARAMETER SYSTEMS by NORMAN STEPHEN BORDAN B. Eng., S i r George Williams U n i v e r s i t y , 1968 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE We accept t h i s t h e s i s as conforming to the required standard Research Supervisor Members of Committee Head of Department Members of the Department of E l e c t r i c a l Engineering THE UNIVERSITY OF BRITISH COLUMBIA May, 1971 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. 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 c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head o f my Department or by h i s r e p r e s e n t a t i v e s . I t i s understood t h a t c o p y i n g or p u b l i c a t i o n of t h i s 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 allowed without my w r i t t e n p e r m i s s i o n . _ E l e c t r i c a l Engineering Department o f _ 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 25th May 1971 /ABSTRACT / / This thesis deals with the i d e n t i f i c a t i o n of parameters i n d i s t r i b u t e d parameter systems. Two s e n s i t i v i t y methods, namely; Meissinger's method and the method of s t r u c t u r a l s e n s i t i v i t y are extended to obtain the s e n s i t i v i t y c o e f f i -cients of d i s c r e t i z e d d i s t r i b u t e d parameter systems. The method of Bingulac and Kokotovic i s extended to i d e n t i f y parameters i n the one and the two dimen-s i o n a l p a r a b o l i c d i f f e r e n t i a l equations. CSMP (continuous system modeling programme) i s used throughout to simulate the systems. Results f o r both s e n s i t i v i t y schemes are obtained, and i t i s found that although s t r u c t u r a l s e n s i t i v i t y i s advantageous f or parameter i d e n t i f i c a t i o n i n ordinary d i f f e r e n t i a l equations, t h i s i s not the case f or p a r t i a l d i f f e r e n t i a l equations. i i / TABLE OF CONTENTS Page. LIST OF ILLUSTRATIONS ^ - v ACKNOWLEDGEMENT / > v i i 1. INTRODUCTION 1.1 I d e n t i f i c a t i o n 1 1.2 S e n s i t i v i t y Analysis 4 1.3 Parameter Optimization 4 1.4 The Problem Considered 5 1.5 Previous Work Done on I d e n t i f y i n g D i s t r i b u t e d Parameter Systems 5 1.6 Scope of the Thesis 6 2. DISCRETIZATION OF PARTIAL DIFFERENTIAL EQUATIONS 8 2.1 Introduction 8 2.2 S p a t i a l D i s c r e t i z a t i o n 8 2.3 Temporal D i s c r e t i z a t i o n 10 2.4 Space Time D i s c r e t i z a t i o n 11 2.5 Conclusion 12 3. SENSITIVITY ANALYSIS 14 3.1 Introduction 14 3.2 Meissinger's Method 14 3.3 A p p l i c a t i o n of Meissinger's Method to the One Dimensional 15 Di s t r i b u t e d Parameter System 3.3.1 A p p l i c a t i o n of Meissinger's Method to the Two Dimensional 16 Di s t r i b u t e d Parameter System 3.4 S t r u c t u r a l S e n s i t i v i t y 24 3.5 Conclusion i i i Page 4. IDENTIFICATION 25 4.1 Introduction 25 4.2 The Method of Bingulac and Kokotovic 25 4.3 I d e n t i f i c a t i o n of a One Dimensional D i s t r i b u t e d Parameter 26 System Using the Method of Bingulac and Kokotovic 4.3.1 Boundary and I n i t i a l Conditions 28 4.3.2 Example #1 29 4.3.3 Example #2 - S t r u c t u r a l S e n s i t i v i t y 33 4.3.4 Example //3 34 4.3.5 Example #4 34 4.4 I d e n t i f i c a t i o n of Two Dimensional D i s t r i b u t e d Parameter 40 Systems 4.5 Conclusion 42 48 5. CONCLUSION 5.1 Summary 4 ^ 49 5.2 Suggestions f o r Further Work APPENDIX A: CSMP 50 A . l Introduction 50 A. 2 The Method of Bingulac and Kokotovic Using CSMP . 50 APPENDIX B: ANALOGUE COMPUTER WORK 51 B. l Introduction 51 B.2 Track and Hold C i r c u i t 51 B.3.1 Meissinger's Method Implementation ill 51 B.3.2 Meissinger's Method Implementation #2 55 B.4 S t r u c t u r a l S e n s i t i v i t y 60 B.5 Use of the Hybrid Computer 60 LIST OF ILLUSTRATIONS Page Figure - Page 1.1 Process Optimization 2 1.2 Model B u i l d i n g 3 1.3 Parameter Optimization 5 2.1 Temporal D i s c r e t i z a t i o n 11 3.1 Three Loop Kokotovic System 19 3.2 S p a t i a l l y D i s c r e t i z e d D i s t r i b u t e d 20 3.3 Kokotovic D i s t r i b u t e d Parameter System 21 3.4 S t r u c t u r a l S e n s i t i v i t y 23 4.1 The Method of Bingulac and Kokotovic 27 4.2 Mode Control 27 31 32 35 36 37 38 39 43 44 45 46 B . l 52 Analogue Computer Implementation #1 of the Method of Bingulac B'.2 . . . . . 53 and Kokotovic B.3 54 4.3 Example #1 (a^ versus a ) 4.4 Example it 2 (a^ versus a^) 4.5 Example #3 (a^ versus a ) 4.6 Example //3 (a^ versus a^) 4.7 Example #4 (a^, a 2 versus no. 4.8 Example #4 ( a ^ a ^ versus no. 4.9 Example #4 (a,, versus no. of 4.10 Example #5 (a^ versus a^) 4.11 Example it 5 (a^ versus a^) 4.12 Example it 5 (& 1 versus B 2 ) 4.13 Example it 5 (B 2 versus B3> V Figure Page B.4 Track and Hold C i r c u i t 56 B.5 External Mode C o n t r o l l e r 57 B.6 Modes of Operation 56 ^"^ Analogue Computer Implementation #2 of the Method of Bingulac and*'* B.8 Kokotovic ... 59 v i .ACKNOWLEDGEMENT This work was supported by a National Research Council Scholarship i n 1968-69 and by a student a s s i s t a n t s h i p i n 1969-71. I am indebted to Dr. A.C. Soudack f o r h i s help and encouragement. I would l i k e to thank Dr. M.S. Davies for reading the manuscript and i n d i c a t i n g improvements. Special thanks are due to Mr. Koit Teng f o r hi s help with CSMP. F i n a l l y , I would l i k e to thank Miss L. Morris f or typing t h i s t h e s i s . v i i 1. INTRODUCTION 1.1 I d e n t i f i c a t i o n In order to con t r o l a process, i t i s necessary to obtain a mathematical model of that process. The procedure by which t h i s model i s obtained i s c a l l e d i d e n t i f i c a t i o n . When the input and the response of the system are used i n the i d e n t i f i c a t i o n p r o c e d u r e ^ , rather than some test s i g n a l , the procedure i s known as "on- l i n e " i d e n t i f i c a t i o n . Many procedures e x i s t f o r the i d e n t i f y i n g , tracking, and estimating of system states and (or) parameters^ 1 ^ Applications of parameter i d e n t i f i c a t i o n l i e i n the txvo areas of process optimization and model b u i l d i n g . In process optimization, the para-meters are chosen to optimize the response of a system.with respect to a ' p a r t i c u l a r c r i t e r i a . This can be seen i n Figure 1.1. For model b u i l d i n g , one would l i k e to determine a model of the system so that the model behaviour and that of the system approximate one another. The model parameters are chosen to accomplish the optimum match. This can be seen i n Figure 1.2. The model b u i l d i n g approach i s considered i n t h i s t h e s i s . Let the unknown system parameter vector be A. A model i s constructed with parameter vector rep l a c i n g the unknown system parameter vector. The tracking parameter vector a i s determined by minimizing an error function f ( e ) . M A * ... (12-17,27) Many procedures e x i s t to minimize an erro r function . The (28) gradient or steepest descent method i s used to i d e n t i f y the unknown parameter vector i n t h i s work. If y i s an element of a, the steepest descent equation i s dF = - K " a e ~ 17 f o r K > 0 (1.1) We w i l l define e as the output err o r i . e . the error between the model / / Pa ramete r Ad j u s tmen t Sys tem Pe r f o rmance c r i t e r i a c a l c u l a t i o n C o m p u t a t i o n o f t h e optimum pa ramete r v a l u e s F i g u r e 1 *1 / SYSTEM PERFORMANCE INDEX' MODEL PARAMETER UPDATE FIGURE 1.2 4 response and that of the system Y. Since the system response i s not a function of a, 8e/9a = 3Y /3a. m The next section w i l l introduce us to the methods f o r obtaining SY^/Sa, the gradient of the model response with respect to the parameters under study. 1.2 S e n s i t i v i t y Analysis S e n s i t i v i t y a nalysis arose from the need to consider the deviations from some set value of the parameters i n a co n t r o l system. It i s possible to e s t a b l i s h an i n d i r e c t r e l a t i o n s h i p between the parameter v a r i a t i o n s and the r e s u l t i n g deviation of the object function. This r e l a t i o n i s obtained by means of parameter s e n s i t i v i t y c o e f f i c i e n t s ^ " * " ^ , which are defined as the gradient of the model response with respect to the parameters under study. Meissinger's method for obtaining the parameter s e n s i t i v i t y c o e f f i c i e n t s i s the basis f or the work which uses the model structure to generate the s e n s i -. . ... . (8-11) t i v i t y c o e f f i c i e n t s With the knowledge of the s e n s i t i v i t y c o e f f i c i e n t s , the tracking parameters y i n equation (1.1) can be obtained, as w i l l be seen i n the next se c t i o n . 1.3 Parameter Optimization The parameters are obtained by using the i t e r a t i v e scheme developed ( 26 ) by Bingulac and Kokotovic . This method consists of successive i t e r a t i o n s to determine the parameters which minimize the performance function. This i s summarized i n Figure 1.3. 5 System M o d e l 1 Y 3^ Pe r f o rmance f u n c t i o n Ym S t e e p e s t d e s c e n t d e t e r m i n a t i p n o f optimum p a r a m e t e r s P a r ame te r a d j u s t m e n t between i t e r a t i o n s Figure 1.3 1.4 The Problem Considered This t h e s i s deals with the "on-line" i d e n t i f i c a t i o n of d i s t r i b u t e d parameter systems. The model b u i l d i n g scheme i s used, and so the form of the system i s assumed known. The d i s t r i b u t e d parameter systems under study are the one and the two dimensional parabolic d i f f e r e n t i a l equations i n which the unknown parameters are s p a t i a l l y dependent. The one dimensional parabolic d i f f e r e n t i a l equation with s p a t i a l l y dependent parameters i s the heat equation f o r a rod which i s constructed cf an inhomogeneous medium, i . e . the s p e c i f i c heat of the material i s s p a t i a l l y dependent. The two dimensional parabolic d i f f e r e n t i a l equation under consideration i s the heat equation for a slab of inhomogeneous mat e r i a l . The object of the work i s to i d e n t i f y the s p a t i a l l y dependent s p e c i f i c heat parameters. 1.5 Previous Work Done on I d e n t i f y i n g D i s t r i b u t e d Parameter Systems I d e n t i f i c a t i o n of d i s t r i b u t e d parameter systems has been attempted (18—2 3) (20) by various authors . Per d r e a u v i l l e and Goodson have extended / Shinbrot's method to make i t applicable to the i d e n t i f i c a t i o n of p a r t i a l d i f f e r e n t i a l equations. This scheme consists of obtaining an i n t e g r a l transform i n the s p a t i a l domain, which when applied to the p a r t i a l d i f f e r e n t i a l / equation and integrated by parts y i e l d s an algebraic equation. This method can not be applicable to a p a r t i a l d i f f e r e n t i a l equation with s p a t i a l l y dependent parameters, because no i n t e g r a l transform can be found i n the s p a t i a l domain. (22 23) Others ' have i d e n t i f i e d parameters i n d i s t r i b u t e d parameter systems by s t a t i s t i c a l means. An extensive l i t e r a t u r e search has been done and i t has been found that no one has i d e n t i f i e d s p a t i a l l y dependent parameters i n the p a r a b o l i c d i f f e r e n t i a l equation. 1.6 Scope of the Thesis Certain classes of ordinary d i f f e r e n t i a l equations have been i d e n t i f i e d using the method of Bingulac and Kokotovic. This method requires the s e n s i t i v i t y of the model response with respect to the parameters under study. Two methods are employed to obtain the s e n s i t i v i t y c o e f f i c i e n t s , namely, Meissinger's method, and that of s t r u c t u r a l s e n s i t i v i t y . This thesis extends the method of Bingulac and Kokotovic, Meissinger's method; and the method of s t r u c t u r a l s e n s i t i v i t y to the i d e n t i f i c a t i o n of s p a t i a l l y dependent parameters i n parabolic p a r t i a l d i f f e r e n t i a l equations. The d i s t r i b u t e d parameter systems were d i s c r e t i z e d s p a t i a l l y and the above method were employed. Chapter 2 deals with the schemes a v a i l a b l e to d i s c r e t i z e p a r t i a l d i f f e r e n t i a l equations. Chapter 3 develops Meissinger's method and the method of s t r u c t u r a l s e n s i t i v i t y f o r d i s c r e t i z e d d i s t r i b u t e d parameter systems. It i s found that 7 s t r u c t u r a l s e n s i t i v i t y i s of l i t t l e use when applied to the systems under study. Chapter 4 deals with the i d e n t i f i c a t i o n of d i s t r i b u t e d parameter systems using the method of Bingulac and Kokotovic and the s e n s i t i v i t y methods developed i n Chapter 3. The S/360 continuous modeling programme used on the IBM 360/67 computer i s applied to s p e c i f i c examples. Chapter 5 consists of a r e c a p i t u l a t i o n of the high points and a suggestion f o r future work. 8 2. DISCRETIZING PARTIAL DIFFERENTIAL EQUATIONS 2.1 Introduction The purpose of t h i s chapter i s to introduce the various schemes a v a i l a b l e for d i s c r e t i z i n g p a r t i a l d i f f e r e n t i a l equations, and to discuss the pros and cons of the methods studied. Let us consider the one dimensional heat equation with the s p a t i a l l y dependent parameter a ( x ) . ifl£Ltl-' t t(x) aVy) ( 2 . i ) o t „ 2 3x where xe(0,x^) and te(0,t^) with boundary conditions q(x,t)[ = L(t) and q ( x , t ) l = M(t) (2.2) and i n i t i a l conditions q ( x , t ) | t = ( ) = K(x) (2.3) Equation (2.1) has two independent v a r i a b l e ; the s p a t i a l v a r i a b l e x, and the temporal v a r i a b l e t. I f we were to d i s c r e t i z e equations (2.1) -(2.3) we could d i s c r e t i z e e i t h e r or both of these v a r i a b l e s . The basic d i s c r e t i z a t i o n schemes are: (i ) S p a t i a l d i s c r e t i z a t i o n . . The s p a t i a l v a r i a b l e i s d i s c r e t i z e d . ( i i ) Temporal d i s c r e t i z a t i o n . The temporal v a r i a b l e i s d i s c r e t i z e d . ( i i ) Space time d i s c r e t i z a t i o n . Both the s p a t i a l and the temporal v a r i a b l e s are d i s c r e t i z e d . Let us now consider the pros and cons of the various methods. 2.2 S p a t i a l D i s c r e t i z a t i o n To obtain the s p a t i a l l y d i s c r e t i z e d model we use the c e n t r a l d i f f e r e n c i n g (25) formula 3q 2(x,t)=q (t) - 2q (t) + q (t) Itl J J _ i ( 2 < 4 ) 3x (Ax) Z 9 where q ^ + ^ ( t ) = q(x + Ax,t) q (t) = q(x,t) q j _ 1 ( t ) = q(x - Ax,t) The d i s c r e t i z e d version of equations (2.1) - (2.3) becomes q,(t) = a . ( q , + 1 ( t ) - 2q (t) + q (t))/(Ax)' (2.5) where j = l , . . , n with boundary conditions q Q ( t ) = L ( t ) and q n + 1 ( t ) = M(t) and i n i t i a l conditions (2.6) (2.7) We now have n f i r s t order equations which represent the p a r t i a l d i f f e r e n t i a l equation. This can be represented by the following vector d i f f e r e n t i a l equation Q = AAQ + ABQB where Q i s an (nxl) vector, A i s an (nxn) matrix, X i s an (nxn) matrix, B i s an (nx2 matrix, and Q i s a (2x1) vector with the following form: ) Q = -2 1 0 1 0 1 -2 1 0 0 ; x = 1 -2 1 ; B = . , 0 ' ' 1 0 1 QB " I q0 V l To implement t h i s d i s c r e t i z a t i o n scheme on the analogue computer one requires n i n t e g r a t o r s . The more accuracy one desires, the greater the number of integrators needed, since/^c gets smaller and therefore n gets larger. 2.3 Temporal D i s c r e t i z a t i o n The temporal d i s c r e t i z a t i o n scheme i s a r r i v e d at by the backward d i f f e r e n c i n g d i s c r e t i z a t i o n Sq^q (x) - q^ C x ) 3.t • A t The d i s c r e t i z e d version of equations (2.1) - (2.3) becomes 2 a(x) 9 q (x) q.(x) - q (x) — ji.i. . —- J (2.10) 9x where q_. = q(x,.,t) and j=l, n with boundary conditions At (2.11) q.(0) = L. and q . (x.) = M. 3 3 3 f 3 and i n i t i a l conditions q Q(x) = K(x) (2.12) (2.13) Temporal d i s c r e t i z a t i o n requires that q^._^(x) be stored from the previous c a l c u l a t i o n so that q.j'(x) ^ s the only unknown i n equation (2.11). This i s the reason why backward d i f f e r e n c i n g i s used i n equation (2.10). The above equations can be represented by the following vector d i f f e r e n t i a l equation (2.14) 2 Ata(x)-^-| = AQ + ABQB dx where Q = 1 0 -1 1 -1 1 0 -1 B = VI qo V l -1 0 0 1 Implementing the temporal d i s c r e t i z a t i o n scheme on the analogue computer presents a few d i f f i c u l t i e s . Continuous memory i s needed, since the 11 q^_^ function i s required i n the s o l u t i o n of equation (2.11). Since t h i s i s an analogue computer scheme, i n i t i a l conditions on a l l the integrators must be obtained. Both the i n i t i a l conditions on 8q/Sx and q are not usually given and so they must be determined by using an i t e r a t i v e scheme. This method employs a feedback loop with four operational a m p l i f i e r s , as can be seen i n Figure 2.1. I f the loop gain exceeds unity, the c i r c u i t becomes unstable due to the p o s i t i v e feedback introduced into the c i r c u i t by (25) the four operational a m p l i f i e r s . The s t a b i l i t y requirement i s l/a(x)At<l The d i f f i c u l t i e s j u s t o u t l i n e d make.temporal d i s c r e t i z a t i o n unattractive for use with t h i s problem. c o n t i n u o u s memory T . A t ; l Figure 2.1 2.4 Space Time D i s c r e t i z a t i o n I f a l l the independent v a r i a b l e s i n a p a r t i a l d i f f e r e n t i a l equation are approximated by f i n i t e d i f f e r e n c e expressions, the p a r t i a l d i f f e r e n t i a l equation i s replaced by a system of simultaneous algebraic equations which can be solved on a d i g i t a l computer. The time d e r i v a t i v e can e i t h e r be forward or backward differenced. 12 Forward d i f f e r e n c i n g y i e l d s q. " q. • = o.(q... . -2q. . + q. , ,)At/Ax 2 (2.15) with boundary conditions and i n i t i a l conditions q o j = L j a n d V i j = M j ( 2 ' 1 6 ) q i 0 = K. .(2.17) (24) Backward d i f f e r e n c i n g y i e l d s q. . ,-j - q. . = a.(q... ... - 2q. .+1 + q. . . j A t / A x 2 (2.18) The forward d i f f e r e n c i n g method, i . e . equations (2.15) - (2.17), i s (29) u s u a l l y solved e x p l i c i t l y on the d i g i t a l computer. Forsythe and Wasow have proven that the forward d i f f e r e n c i n g method i s stable and the d i s c r e t i z a t i o n 2 e r r o r i s of the order of (Ax) provided that a^>0 and At$Ax 2/2a ± (2.19) I f these r e l a t i o n s are not s a t i s f i e d , i n s t a b i l i t y r e s u l t s due to the propagation of roundoff and truncation e r r o r s . The i m p l i c i t d i f f e r e n c i n g method using equations (2.16) - (2.18) has no s t a b i l i t y requirement, but requires the s o l u t i o n of nxm simultaneous equations (2.18). This can be accomplished by using Gauss E l i m i n a t i o n or P i v o t a l techniques. For many nodes t h i s i s very time consuming. 2.5 Conclusion It was o r i g i n a l l y hoped to use the Pace 231-R analogue computer f o r the work done i n t h i s t h e s i s . The analogue computer could not be used due to an i n s u f f i c i e n c y i n the number of m u l t i p l i e r s and inte g r a t o r s a v a i l a b l e i n the machine. This d i f f i c u l t y w i l l be explained f u l l y i n Appendix B. 13 To circumvent the' d i f f i c u l t i e s with the analogue computer, the s/360 continuous modeling programme (CSMP) was used on the IBM 360/67 computer. CSMP i s a problem oriented programme designed to f a c i l i t a t e the d i g i t a l simulation of continuous processes. The continuous process i s d i s c r e t i z e d by the programme. A complete explanation, of CSMP appears i n Appendix A. Since the d i g i t a l computer was used a l l the var i a b l e s had to be d i s c r e t i z e d , and so, space time d i s c r e t i z a t i o n had to be used. The problem was looked at as though the s p a t i a l v a r i a b l e was d i s c r e t i z e d and the temporal v a r i a b l e was maintained i n the continuous form. The CSM program was then used to d i s c r e t i z e the temporal v a r i a b l e . Since both the s p a t i a l and the temporal v a r i a b l e s were ul t i m a t e l y d i s c r e t i z e d , equation (2.19) had to be s a t i s f i e d . Though we w i l l c a l l the d i s c r e t i z a t i o n scheme used s p a t i a l d i s -c r e t i z a t i o n and w i l l use the equations developed i n Section (2.2), the space time d i s c r e t i z a t i o n scheme i s r e a l l y what i s obtained using CSMP, as can be seen from the preceeding argument. In the next chapter the parameter s e n s i t i v i t y c o e f f i c i e n t s are obtained f o r the s p a t i a l l y d i s c r e t i z e d p a r t i a l d i f f e r e n t i a l equations obtained i n t h i s chapter. 14 3. SENSITIVITY ANALYSIS 3.1 Introduction This chapter considers two methods to obtain the parameter s e n s i t i v i t y c o e f f i c i e n t s , namely; Meissinger's method, and the method of s t r u c t u r a l s e n s i t i v i t y . These methods w i l l be extended to obtain the parameter s e n s i -t i v i t y c o e f f i c i e n t s of the d i s c r e t i z e d parabolic d i f f e r e n t i a l equations obtained i n Section (2.2) . Though the method of s t r u c t u r a l s e n s i t i v i t y i s u s e f u l f o r ordinary d i f f e r e n t i a l equations, i t w i l l become evident that there i s l i t t l e advantage to i t over the other method when applied to d i s t r i b u t e d parameter systems. 3.2 Meissinger's Method Meissinger's method f o r the generation of the parameter s e n s i t i v i t y c o e f f i c i e n t s can be described by the following example. Consider the second order d i f f e r e n t i a l equation Y(t) + a^ C t ) + a Q Y ( t ) = X(t) (3.1) where X(t) i s the f o r c i n g function. Meissinger's parameter influence c o e f f i c i e n t s are defined as V1 = %Y/da1 and U 0 = 9 Y / 9 a Q (3.2) D i f f e r e n t i a t i n g equation (3.1) p a r t i a l l y with respect to a^ and aQ, we obtain u'l + a ^ + a 0"d 1 = -Y (3.3) "o + a A ) + ao uo - " Y ( 3 - 4 ) To obtain the s e n s i t i v i t y c o e f f i c i e n t s , equations (3.3) and (3.4) must be solved. These equations are s t r u c t u r a l l y the same as equation (3.1) but have d i f f e r e n t inputs. Let us now apply Meissinger's method to the d i s c r e t i z e d d i f f e r e n t i a l equation obtained i n the previous chapter. 15 3.3 A p p l i c a t i o n of Meissinger's Method to the One Dimensional Parabolic D i f f e r e n t i a l Equation In Section (2.2) i t was found that a one dimensional parabolic p a r t i a l d i f f e r e n t i a l equation can be put i n the following vector form Q = AAQ + ABQ_ B where A i s a matrix of the unknown parameters. Let a be an element of A. D i f f e r e n t i a t i n g equation (2.8) with respect to awe obtain |2 . | j [AXQ + ABQ B] - f W + A | i Q + AX |2 + | A BQ + i f n + AB ^ ,. 3 a 3 a 3 a 3 a 3 a B 3cj (3.5) Let U = 3 Q / 3 a be c a l l e d the s e n s i t i v i t y c o e f f i c i e n t vector. Since X, B, and Q_ are not functions of a. We obtain 3A U = ~ • AQ + AA U (3.6) O Oi Applying Meissinger's method to the d i s c r e t i z e d version of equation (2.8), namely equations (2.5) - (2.7) we obtain where U. . = 3 q . / 3 a . (3.8) with boundary conditions U .(t) | = 0. and U. , ( t ) | = 0. (3.9) 1 , J i=0. 1 , 3 i=n+l" and i n i t i a l conditions U .(0) = 0. (3.10) 1 > J 16 where i = l , . .., n and j = l , . . . , n and / 6 i j = l - f o r i=j / / 61j = 0 f o r i ^ j / I f equation (2.5) consists of n f i r s t order equations, Meissinger's method requires the generation of nxn f i r s t order s e n s i t i v i t y equations i n order to a r r i v e at a l l the parametric s e n s i t i v i t y c o e f f i c i e n t s . Let us now apply Meissinger's method to the two dimensional parabolic d i f f e r e n t i a l equation. 3.3.1 A p p l i c a t i o n of Meissinger's Method to the Two Dimensional Parabolic D i f f e r e n t i a l Equation Let us consider the two dimensional parabolic d i f f e r e n t i a l equation 2 2 | 1 (x,y,t) = a(x) (x,y,t) + g(y) (x,y,t) (3.11) 9x 9y where xe(0,x f), y e ( 0 , y f ) , t e ( 0 , t f ) with boundary conditions q ( x , y , t ) | x = 0 = L 1 ( y , t ) and q(x,y,t)| Q = L 2 ( x , t ) (3.12) 'q(x-,y,t)L = M (y,t) and q(x,y,t)| = M ( x , t ) ^ f y - ^ ^ z-and i n i t i a l conditions q(x,y,t) | t = ( ) = K(x,y) (3.13) We can d i s c r e t i z e equation (3.9) with respect to the two s p a t i a l v a r i a b l e s x and y. Using the d i s c r e t i z a t i o n developed i n Section (2.2) we obtain - ~2 - 2 « i + ' i - i ' j - 5 [vi - 2 q i + q j - i ] i ( 3- l 4 > 17 where i = l , . .., n and j = l , m. Let us pick m=n j u s t f o r s i m p l i c i t y . Applying Meissinger's method to equation (3.10) we obtain u. . , = —o[,u-+v- 2u. + u. J . + -^ L- tu... - 2u. + u. k i , j , k w 2 l + l l l - l 3,k Arr2 j+1 j j-1 -i> K . . (3.15) Ax Ax / Ay / and a • ft • U. . = t U. - 2U. + U. . 1 , + [ U. ,,-211. + U. . ] . i.J.m A x 2 i + l i i - l j'.,m ^ y2 j+1 j j-1 i,m (3.16) where U . . = 3q../3a, and U. , = 3q../33 are defined as the parameter i>j ,k ^ i j k i,k,m 4 i j m s e n s i t i v i t y c o e f f i c i e n t s , and where 6. ... =1. when i=i=k 13k J = 0. else and <S. . =1. when i=j=m ijm = 0. else I f equation (3.14) consists of (n xn) f i r s t order equations, xoe can 2 see that Meissinger's method requires us to generate (2><nxn ) f i r s t order s e n s i t i v i t y equations to a r r i v e at a l l the parametric s e n s i t i v i t y c o e f f i c i e n t s . The next section deals with the a p p l i c a t i o n of s t r u c t u r a l s e n s i t i v i t y to the preceding problem. 3.4 S t r u c t u r a l S e n s i t i v i t y The purpose of t h i s section i s to extend s t r u c t u r a l s e n s i t i v i t y so that i t w i l l be useful i n determining the s e n s i t i v i t y functions of a d i s c r e t i z e d p a r t i a l d i f f e r e n t i a l equation. 18 Kokotovic has shown that f o r a class of systems with one d i r e c t path a l l the s e n s i t i v i t y functions may be obtained simultaneously with the (8) use of one a d d i t i o n a l system model . This i s superior to Meissinger's method because Meissinger's method requires the generation of a system model for each parameter, as seen i n Section (3.2). The methods developed by K o k o t o v i c ^ , Vuscovic and C i r i c ^ " ^ t t l deal only with systems where the m branch transmittance i n a s i g n a l flow graph i s W (s, u^) and where ct^, the parameter of i n t e r e s t only appears i n that branch. S t r u c t u r a l s e n s i t i v i t y w i l l be extended to those systems where the parameters of i n t e r e s t appear i n more than one branch transmittance i . e . a. might appear i n W,(s, a.) and W (s, a., a.). i K. i . 1 i j] The class of systems considered are those i n which the outputs of the feedback paths are brought together and subtracted from the input and there i s one feedforward path. This class of systems i s c a l l e d the Kokotovic c l a s s . A three loop Kokotovic system i s shown i n Figure 3.1. A s p a t i a l l y d i s c r e t i z e d parabolic d i f f e r e n t i a l equation i s not of the Kokotovic c l a s s , as can be seen i n Figure 3.2. By applying s i g n a l flow graph theory to the d i s c r e t i z e d system, i t can become of the Kokotovic class provided that one of the boundary conditions i s set to zero (one feedforward path). The parameters now appear i n a more complicated form i n the transmittance i n the feedback paths, as can be seen i n Figure 3.3. The s e n s i t i v i t y functions considered w i l l be that of equation (3.8) i . e . A U. . = H /8a. With the p a r t i a l d i f f e r e n t i a l equation i n the Kokotovic class i t i s evident that the system equation becomes 19 X. (S> W 1(S,a 1) W 2(S,a 2) W 3(S,a 3) W 4(S,a 4) W 5(S,a 5) y(s) W 6(S,a 6) Figure 3.1 FIGURE. 3.2 FIGURE 3.3 to 22 q-j(s) = W.(s,a ,...,a )q (s) (3.17) J j n o where W.(s,a., ...,a ) i s the o v e r a l l t r a n s f e r function needed to obtain j 1 n the output q^(s), and where q Q ( s ) i s the input (one of the boundary condi t i o n s ) . Subsituting equation (3.17) i n t o equation (3.8) we obtain : 9W.(s,a * a ) U..(s)=rr-7 r- ' -1 - 1 — d.(s) (3.18) i j Wj(s a 1...a n) 9a. j I f the parameter cx\ appears i n more than one loop t r a n s f e r function, we must apply the chain rule to equation (3.18) to obtain a. n 9W.(s,a_...a -\ 9W, ( s , a ) x,j W (s.o^...a n) k = 1 9W k(s,a i) 9a ± n n = q i E B i k ' C k i = Z v i . J . k J k = 1 J , * fc.i k = 1 where and A wk 3 w i Bj } k = w7 w- . <3-20) J i k r A _1 k k , i W. 9a u ' (3.21) k l The functions E. .depend on the system's structure, and the functions .1>K C, ..depend on the structure of the l i n k s W. . k , i k I f W. (s) = ~ _ then C, k v ' s+2a. ct.(s+2a,) X - X i Equation (3.20) defines the t r a n s f e r functions from the system's input q. to the point S, ... These points are c a l l e d the s e n s i t i v i t y points. To determine the s e n s i t i v i t y functions U. ., i t s u f f i c e s to send the s i g n a l s from the s e n s i t i v i t y points S ' to the blocks C ..and sum the k,j k,x r e s u l t s , as i n equation (3.19). A l l the s e n s i t i v i t y functions are determined y - -C f^t- V A-l C>32 FIGURE 3.4 / 24 by varying the input q_. of the. s e n s i t i v i t y model. This can be seen i n Figure 3.4. Since the system under study does not have a s i n g l e output (as i s the case with ordinary d i f f e r e n t i a l equations), but has n outputs, the number of equations required to generate the s e n s i t i v i t y c o e f f i c i e n t s i s not reduced by using s t r u c t u r a l s e n s i t i v i t y . 3.5 Conclusion S t r u c t u r a l s e n s i t i v i t y can not be applied to p a r t i a l d i f f e r e n t i a l equations i n general. The systems considered using t h i s method must be of the Kokotovic c l a s s . This implies putting the systems i n t o the Kokotovic form, and hence l i m i t i n g the systems studied to those i n which only one boundary condition i s considered. Since the system equation has n outputs, s t r u c t u r a l s e n s i t i v i t y requires us to solve (n*n) equations i . e . the same number of equations as i n Meissinger's method and hence no saving r e s u l t s . The only advantage to the a p p l i c a t i o n of s t r u c t u r a l s e n s i t i v i t y i s that once the s e n s i t i v i t y model i s obtained we need only vary, the input to the model to obtain a l l the s e n s i t i v i t y c o e f f i c i e n t s . Meissinger's method has none of the above const r a i n t s . The next chapter w i l l deal with the i d e n t i f i c a t i o n of parameters using the method of Bingulac and Kokotovic and the s e n s i t i v i t y methods developed i n t h i s chapter. 25 4. PARAMETER IDENTIFICATION IN DISTRIBUTED PARAMETER SYSTEMS 4.1 Introduction In t h i s chapter the s e n s i t i v i t y methods developed i n the previous chapter are used i n conjunction with the method of Bingulac and Kokotovic i n order to i d e n t i f y parameters i n the f i r s t and second order parabolic d i f f e r e n t i a l equations. The method of Bingulac and Kokotovic i s presented i n the next se c t i o n . 4.2 The Method of Bingulac and Kokotovic The method of Bingulac and Kokotovic i s applicable to parameter i d e n t i f i c a t i o n problems using a mode c o n t r o l l e d analogue computer (or i n our case CSMP). T The parameter vector a = [a^ . . a^ ] i s adjusted by means of steepest descent i n order to minimize a performance index J(a) of the form ; T J(a) = ff F(e)dt dfi (4.1) so The s e n s i t i v i t y methods developed i n the previous chapter are used i n obtaining T | ^ = // | £ U dt dQ (4.2) 8a Se so It i s assumed that the parameter vector a i s held constant i n deri v i n g equation (4.2). For th i s reason the computer i s operated i n two modes. In the compute mode, the parameter vector i s held constant while the constitutents of equation (4.2) are ca l c u l a t e d . In the reset mode the parameter vector i s adjusted according to the steepest descent law a j ~ " K 9a. (4.3) 26 This method i s summarized i n Figures (4.1) and (4.2). The increments i n ou are required to be kept small because the d i r e c t i o n of steepest descent i s steepest only i n the v i c i n i t y of a ^ . The increments i n can be kept small by an appropriate choice of K i n the above equation. A compromise must be struck between the s t a b i l i t y of the method and the speed of convergence. The l a r g e r the value of K i s , the f a s t e r i s the rate of convergence, and the greater i s the r i s k of i n s t a b i l i t y . Let us consider now, the i d e n t i f i c a t i o n of the one dimensional p a r a b o l i c d i f f e r e n t i a l equation. 4.3 I d e n t i f i c a t i o n of the One Dimensional Parabolic D i f f e r e n t i a l Equation We would l i k e to i d e n t i f y the unknown parameters of the following d i s t r i b u t e d parameter system, | | (x,t) = 3(x) (x,t) (4.4) 3x with boundary and i n i t i a l conditions previously used i n equations (2.5) - (2.6) and where 3(x) i s the unknown s p a t i a l l y dependent parameter. Let the model be 2 |~j-(x,t) = <x(x) (x,t) (4.5) dx with the same i n i t i a l and boundary conditions as above and where a(x) i s a chosen parameter. Let the performance function be X f fcf 2 J = f f (q(x,t) - s ( x , t ) ) dxdt (4.6) o o which i s quadratic and covers the e n t i r e space time domain. D i s c r e t i z i n g equations (4.4) - (4.6) s p a t i a l l y as i n Section (2.2) T^ e obtain 27 System P e r f o r m a n c e F u n c t i o n Ym S e n s i t i v i t y Mode l E Pa r amete r a d j u s t m e n t between i t e r a t i o n s S t e e p e s t d e s c e n t d e t e r m i n a t i o n o f optimum p a r a m e t e r s F i g u r e 4 r l MODE BLOCK COMPUTE RESET A , B , C , D ; Compute a c c o r d i n g ^ 1 t o (4.2) R e t u r n t o i n i t i a l c o n d i t i o n s f o r t h e nex t i t e r a t i o n E H o l d <Xi c o n s t a n t A d j u s t ° i^ f o r t h e nex t i t e r a t i o n a c c o r d i n g t o (4.3) F i g u r e 4,2 28 d s i 2 ~dt = B i ( s i + l ( t ) " 2 s i ( t ) + S i _ l ( t ) ) / A X ( 4- 7> ~d7 = a . ( q i + 1 ( t ) - 2q.(t) + q . _ 1 ( t ) ) / A x 2 (4.8) J = A X / t f E ( q , - s . ) 2 d t (4.9) ° 1=1 D i f f e r e n t i a t i n g equation (4.9) with respect to a and s u b s t i t u t i n g i n t o equation (4.3) y i e l d s a j = _ K ^ f E (q. - s ) U dt . (4.10) i = l J where j = l , ..., n The s e n s i t i v i t y c o e f f i c i e n t s U „ are obtainable by using e i t h e r of the schemes developed i n Chapter 3. Before dealing with a few examples, a word i s i n order about the choice of boundary and i n i t i a l conditions. 4.3.1 Boundary and I n i t i a l Conditions The boundary and the i n i t i a l conditions are of the form q ( x , t ) | x L q = L ( t ) ; q ( x , t ) | = M(t) (4.11) q ( x , t ) | t = Q = K(x) (4.12) The conditions are assumed not to be functions of the parameters under study. I t i s also assumed that they are given f o r both the system and the model. D i s c r e t i z i n g equations (4.11) - (4.12) s p a t i a l l y we obtain q Q ( t ) = L(t) ; q n + 1 ( t ) = M(t) (4.13) q^(0) = I<L where i = i , n (4.14) 29 The form of the boundary conditions presents no problem provided that each value of ^ ( t ) and ^ ^ ^ ( t ) ^ s introduced into equations (4.7) and (4.8) at the correct point i n time i . e . as the time sweep progresses. This i s e a s i l y obtainable using CSMP. Since i t i s assumed that the parameters under study are not functions of the boundary and i n i t i a l conditions, l e t us choose the easiest conditions to deal with, namely q ( x ' t ) l x = 0 . = ° ' 5 q ( x ' t ) 1 x - x f = 0 (4.15) q ( x , t ) | t = Q = K(x) * (4.16) We are now i n a p o s i t i o n to deal with some examples. 4.3.2 Example #1 In t h i s example Meissinger's method i s used to obtain the s e n s i t i v i t y c o e f f i c i e n t s . Here n=3 i s chosen. The system equation i s dS, dt " = — 2 [ S 1 + 1 " 2S. + S 1_ 1] i = l n (4.17) Ax with boundary conditions S = S. = 0 and i n i t i a l conditions S.(0) =1. J o 4 i The model equation i s Ax with boundary conditions q = q'^ = 0 and i n i t i a l conditions q^(0) = 1. Ther performance function i s t n " J = / £ (q - S.) dt (4.19) 0 i = l 1 The Meissinger s e n s i t i v i t y model i s U i , j = f 2 f U i + 1 - 2 U i + U i - 1 > j + 74 tq-4-l - q , + q±_J (4.20) , J Ax Ax where • - 1 / . - - i , x — x) • • j ri j j JL j * • • ^ n / Where the boundary and i n i t i a l conditions are 1' U. . (0) = 0; U . = U .. = 0 • . (4.21) i j 03 ..4j The update algorithm i s / T n Act. = -K Z (q-_, - S.)U\. dt (4.22) J o ± = 1 i 1 i i The gain constant i n equation (4.23) i s set at K=.02. Since n=3. i s chosen and the increment i n x i s set at .2, the length of rod considered i s .8 u n i t s . Figure 4.3 shows the steepest descent of the c o n t r o l l e r parameters i n the (a^ a^) plane, p l o t t e d against contours of equal J. To obtain the contours of equal J , a.^ i s set at a constant. Figure 4.4 shows the steepest descent of the c o n t r o l l e r parameter parameters i n the (c^, ct.^ ) plane, p l o t t e d against contours of equal J , where the contours of equal J are obtained by setting'ct. to a constant. The i t e r a t i v e procedure was i n i t i a t e d from three locat i o n s to check for l o c a l minima. As can be seen i n Figure (4.3) and (4.4), the i n i t i a l guesses of (.2,.2,.2), (.05,.05,.05) and (.2,.05,.2) converged to the optimum values of (.1,.15,.18). An unfortunate c h a r a c t e r i s t i c of the steepest descent method i s i t s (28") slow convergence i n the region of the optimum I f J has a narrow ridge, i . e . i f the J contours can be approximated by e l l i p s e s with high e c c e n t r i c i t y , then the convergence rate i s even slower^"^ . This can be overcome somewhat by increasing the gain constant K as the optimum i s approached. I 33 4.3.3 Example #2 - S t r u c t u r a l S e n s i t i v i t y In t h i s example the method of s t r u c t u r a l s e n s i t i v i t y i s applied to the preceeding problem with the same choice of K, and i n i t i a l guesses of a. For n=3, we require s i x VI. l i n k s , as seen i n Figure 3.4, where W (S.o^) = a-j/s W 2(S) = 2 W 3(S,a 2) = a 2/(S+2a 2) V s ) = 1 W5(S, a 3) = 03/8+ 2 a 3 w6(s,0l) = s + 2 a i \ Using equation (3.21) we obtain a... = l/'a, C = C = . . C = 0 1 1 1 12 13 15 c 2 1 = ... = c 2 6 = 0 C 3 2=S/a 2 (s+2a 2) = ,.. = = 0 C41 = '* * C46 = ° C 5 3 = S/a 3(S+2a 3) S l = C52 = C54 = ' ' ' C56 = 0 C61 = - S / a l ( S + 2 a l ) C62 = - C66 " ° From equation (3.19) we obtain j=l,-»3 The nine s e n s i t i v i t y equations become U - . = q . [B. C 4-B. C. A = S . ' C. + S. , • C, 1»J J 1,1 11 1,6 J,6 A J , l U J,6 j,6 2, j H J j,3 3,2 j,3 3,2 U„ . = q..[B. • C_ ,] = S. • C, _ 3, J . j 1,5 5,3 3 , 5 5,3 34 The I d e n t i f i c a t i o n scheme using s t r u c t u r a l s e n s i t i v i t y produces p r e c i s e l y the same r e s u l t s and requires the same computing time as i n Example #1 i . e . Figure (4.3) and (4.4). 4.3.4 Example #3 In t h i s example, the same equations were used as i n example #1.. 2 Here 3(x) = x + ax + b where a = -.25, and b = .11. The d i s c r e t i z e d version of 3 becomes (.1, .17, .32). Figure 4.5 and Figure 4.6 show the steepest descent of the con-t r o l l e r parameters i n the (c^., a,.) and (a, , } plane r e s p e c t i v e l y . The i t e r a t i v e procedure was i n i t i a t e d from two locations to check for l o c a l minima. As can be seen i n the Figures, the i n i t i a l guesses of (.05, .05, .05), (..35, .35,.. 35) converged to the optimum value of (.1, .17, .32). In Figure 4.6 a great deal of o s c i l l a t i o n occurs due to the ridge at the minimum. 4.-3.5 Example #4 In t h i s example the same equations were used as i n Example #1. Here n = 5. was chosen. With n = 5. and the increment i n x set at .2, the length of the rod considered i s 1.2 u n i t s . 3(x) was set as a l i n e a r function of x i . e . 3(x) = ax + b where a = -.1 and b = .22. The d i s c r e t i z e d parameters 3 become (.2, .18, .16, .14, .12). The a values converged to these values from the i n i t i a l guess of (.1, .1, .1, .1, .1) w i t h i n twenty i t e r a t i o n s . This can be seen i n Figures (4.7) - (4.9). The greater the number of d i s c r e t i z a t i o n s ; the more equations need be solved i . e . f o r n = 5; 25 s e n s i t i v i t y equations are required to be solved, the beauty of using CSMP i s that there i s no constraint on the number 35 (X2 No. o f I t e r a t i o n s .1 — j 1 1 1 3 r m 20 F i g u r e U .7 ONE DIMENSIONAL PARABOLIC CASE n=5 ; (X(x)=ax+b ~ No. o f I t e r a t i o n s i 1 1 1 r — r — i r- i r 10 20 No. o f I t e r a t i o n s F i g u r e V-S ONE DIMENSIONAL PARABOLIC CASE n=5; cx(x)=ax+b : '. i r_ 1 r 1"0 No. o f I t e r a t i o n s .2 _ ONE DIMENSIONAL PARABOLIC CASE n = 5 ; c*(x)=ax+b <X5 .1 — i — i — < — r 10 No. o f I t e r a t i o n s ~i 1 i r 20 F i g u r e 4 -9 LO -40 of equations to be solved, other than the computer time required. 4.4 I d e n t i f i c a t i o n of the Two-Dimensional Parabolic D i f f e r e n t i a l Equation / We would l i k e to i d e n t i f y the unknown parameters of the following d i s t r i b u t e d parameter system 2 2 | | (x,y,t) = r(x) (x,y,t) + v(y) (x,y,t) (4.23) 9x 9x where xe(0,x^), ye(0,y^), tc(0,t^) with boundary conditions (4.24) s ( x , y , t ) | x = o > = M^y.t) ; s(x,y,t) 1 ^ = ^ ^ s ( x , y , t ) | y = o > = N^x.t) ; s(x,y,t) | y = y ^ = and i n i t i a l conditions s ( x , y , t ) | t = Q _ = P(x,y) Let the model be 2 2 | | (x,y,t) = a(x) (x,y,t) + B(y) -M*- (x,y,t) (4.25) 9x . 9x with the above i n i t i a l and boundary conditions. Let the performance function be t f x y J = f f f 2 o o o (q(x,y,t) - s(x,y,t)) dx dy dt (4.26) D i s c r e t i z i n g equations (4.23) - (4.26) with respect to the s p a t i a l v a r i a b l e s x and y, we obtain r. v, where i = 1, ..n and j = 1, p 41 with boundary conditions 0,J 1,3 n+l,j 2,j S i , 0 = M l , i 5 Si,n+1 = N 2 , i (4.29) and i n i t i a l conditions and performance function s. .(0) = P. . 1,3 i»3 t n m . J = AxAy / .Z. .Z (q, . - s. .) dt (4.30) J o x=l j = l 1,J D i f f e r e n t i a t i n g J with respect to ^ and we obtain f)T t n p = 2AxAy / 1 Z Z (q - s )U dt (4.31) 9ak ° i = l j = l X ' J 1 ' J ' | C Si t n p = 2AxAy / Z Z (q. . - s. .)U. dt (4.32) 3 6m ° i = l j = i x*3 1>3 ^-'J'" 1 S u b s t i t u t i n g equations (4.28) - (4.29) in t o the two dimensional update algorithms Aa. = -K — - • and A3 = -L ~ -k 9a r m 83 k m y i e l d s t n p Aa. = -K / • Z Z (q. . - s. .)U. . . dt (4.33) A3 m = -L• / f z E (q. . - s. .)U. . dt (4.34) ° i = l j = l ' J 1 , 3 1 ' 3 ' m The s e n s i t i v i t y c o e f f i c i e n t s U. . , and U. . are obtained by using equations (3.15) and (3.16). 42 4.4.1 Example #5 In t h i s example Meissinger's method was employed to obtain the s e n s i t i v i t y c o e f f i c i e n t s . Here n = m = 3. was chosen The equations considered were (4.27) and (4.28) with boundary conditions q 0 > J f o r i = l , . n and j =l^m.and i n i t i a l conditions q. . (0) 1. 1,0 These were chosen = o simply for convenience. r(x) and v(y) were set as a l i n e a r function of x and y i . e . r(x) = ax + b and v(y) = ay + b. With a = .25 and b = .05, the d i s c r e t i z e d parameters r^ and become (.1, .15, .2). Figures (4.10) - (4.13) show that the parameter guesses of (.2, .2, .2) and (.05, .05, .05) converged to the optimum. For the two dimensional case Meissinger's method requires the 2 s o l u t i o n of (2xnxn ) equations, i n our case t h i s leads to. 54 s e n s i t i v i t y equations. The problem associated with i d e n t i f i c a t i o n using the method of Bingulac and Kokotovic i s the number of equations that are required to be solved. For the two dimensional case as outlined above; nine system equations are required, nine model equations are required, as w e l l as the 54 s e n s i t i v i t y equations. For t h i s example CSMP requires 62 seconds to compile the program versus the 12 seconds of compilation time needed for example ill. This, of course, i s due to the number of equations that must be solved. 4.5 Conclusion In t h i s chapter f i v e examples of i d e n t i f y i n g parameters i n d i s t r i b u t e d parameter systems using the method of Bingulac and Kokotovic were attempted. Systems for n = 3 were e a s i l y i d e n t i f i e d using Meissinger's method to obtain the s e n s i t i v i t y c o e f f i c i e n t s . It was found that s t r u c t u r a l s e n s i t i v i t y y i e l d e d no economy i n the number of equations that had to be solved, and y i e l d e d exactly the TWO DIMENSIONAL PARABOLIC CASE .20H a c n -, r . 1 • | .10 .21 F i g u r e 4'10 .20 _ J .10 _ TV/0 DIMENSIONAL PARABOLIC CASE n=3 o((x)=ax+b; ft(v)=ay+fa TWO DIMENSIONAL PARABOLIC CASE cx(x)=ax+b; jg(y)-qy+P OPTIMUM -i 1 1 r 0 F i g u r e 4 . 13 0 l 1 1 1 same r e s u l t s as Meissinger's method as can be seen i n Examples //I and //2. For the nonlinear parameter choice i n Example #3, i t was found that a great deal of o s c i l l a t i o n occurred while approaching the optimum. This was due to the ridge i n the performance function. For Example #4 i n which n = 5, the f i v e parameters were i d e n t i f i e d using Meissinger's method, however, the compilation time of CSMP was double that of Example / / l . This resulted because 25 s e n s i t i v i t y equations had to be solved versus 9 s e n s i t i v i t y equations i n Example #1. For the two dimensional case, Example #5, the parameters were e a s i l y i d e n t i f i e d using Meissinger's method,: however compilation time of CSMP was s i x times that of Example //I. This resulted because 54 s e n s i t i v i t y equations had to be solved i n t h i s example. Though the method of Bingulac and Kokotovic can s u c c e s s f u l l y be employed to i d e n t i f y s p a t i a l l y dependent parameters i n d i s t r i b u t e d parameter systems, i t s ultimate f a i l i n g i s i n the number of s e n s i t i v i t y equations that must be solved. 48 • 5. CONCLUSION 5.1 Summary This work has been concerned with the i d e n t i f i c a t i o n of d i s t r i b u t e d / parameter systems of the parabolic type with s p a t i a l l y dependent parameters. The method of Bingulac and Kokotovic, Meissinger's method, and s t r u c t u r a l , s e n s i t i v i t y were extended to i d e n t i f y the above mentioned systems. Chapter 2 deals with the d i s c r e t i z a t i o n schemes a v a i l a b l e to d i s -c r e t i z e p a r t i a l d i f f e r e n t i a l equations. I t was seen that s p a t i a l d i s c r e t i z a t i o n could be considered with CSMP' provided that equation (2.19) i s s a t i s f i e d . Chapter 3 develops Meissinger's method and the method of s t r u c t u r a l s e n s i t i v i t y f o r d i s c r e t i z e d d i s t r i b u t e d parameter systems. I t was found that s t r u c t u r a l s e n s i t i v i t y can only be applied to those systems of the Kokotovic c l a s s ; and t h i s implies that one boundary condition must be set to zero. S t r u c t u r a l s e n s i t i v i t y does not y i e l d an economy i n the number of s e n s i t i v i t y equations that need be solved as i s the case i n ordinary d i f f e r e n t i a l equations. In f a c t , s t r u c t u r a l s e n s i t i v i t y requires the s o l u t i o n of the same number of s e n s i t i v i t y equations as that of Meissinger's method, and so s t r u c t u r a l s e n s i t i v i t y i s of l i t t l e use. Chapter 4 deals with the extension of the method of Bingulac and Kokotovic to the i d e n t i f i c a t i o n of d i s t r i b u t e d parameter- systems using the s e n s i t i v i t y methods developed i n the previous chapter. Both s e n s i t i v i t y methods were employed and i t was found that they y i e l d e d the same r e s u l t s . Though the method of Bingulac and Kokotovic can be s u c c e s s f u l l y employed; i t s ultimate f a i l i n g l i e s with the number of s e n s i t i v i t y equations that need be solved. 49 5.2 Suggestions for Further Research The problem under study i s suited to the hybrid computer and i t would be very i n t e r e s t i n g to approach i t from that point of view. This would avoid the need for analogue m u l t i p l i e r s , track and hold devices, and external mode con t r o l devices, so that the number of equations to be considered i n the analogue portion of•the computer could be increased. 50 APPENDIX A A . l The Continuous Modeling Programme The S/360 CSMP i s a problem oriented programme designed to f a c i l i t a t e the d i g i t a l simulation of continuous processes on large scale d i g i t a l computers. The general CSMP formulation of a model i s divided into three segments: I n i t i a l ; Dynamic; and Terminal; that describe the computations to be performed before, during, and a f t e r each simulation run. The i n i t i a l segment i s intended e x c l u s i v e l y f o r computation of i n i t i a l conditions. ' The dynamic segment consists of the system dynamics. The terminal segment i s used for those computations desired a f t e r completion of each run. A.2 The Method of Bingulac and Kokotovic Using CSMP The method of Bingulac and Kokotovic can e a s i l y be programmed using continuous system modeling. The two modes of operation discussed i n Section (4.2), i . e . the compute and the operate modes are obtained i n the following manner. The compute mode i . e . the mode i n which the parameters are held constant i s simply obtained by s e t t i n g the parameters to a set value i n the i n i t i a l segment of the programme. The reset mode i . e . the adjustment of the parameter values by means of the steepest descent algorithm (equation 4.3), i s obtained by programming the equation i n the terminal segment of the programme. The updated parameter values are then used as the i n i t i a l conditions for the next sx^eep. 51 APPENDIX B: ANALOGUE COMPUTER WORK B . l Introduction This appendix deals with the poss i b l e analogue computer imple-mentations f o r the i d e n t i f i c a t i o n of d i s t r i b u t e d parameter systems, and t h e i r f a i l i n g s . This r e s u l t s from a need f o r a large number of components i n the analogue computer, namely; m u l t i p l i e r s and in t e g r a t o r s . The large number of integ r a t o r s i s needed i n order to b u i l d the track and hold c i r c u i t s required. B.2 The Track and Hold C i r c u i t The track and hold c i r c u i t i s obtained by using a p a i r of mode co n t r o l l e d i n t e g r a t o r s . To obtain a mode co n t r o l l e d i n t e g r a t o r the Memory Logic Group (MLG) unit must be engaged. With the proper configuration on the MLG patch panel, an S = 1 mode pulse w i l l put a designated i n t e g r a t o r i n t o the set mode. An S = 1 pulse w i l l put the int e g r a t o r into the reset mode. A track and hold c i r c u i t can be obtained by using a p a i r of mode co n t r o l l e d i n t e g r a t o r s connected as i n Figure B.4. In the S = 1 mode, int e g r a t o r 1 follows the input and the relay i n the fi g u r e sets the output to ground. In the S = 1(S=0) mode, integrat o r 2 tracks i n t e g r a t o r 1. Since i n t e g r a t o r 1 i s held at i t s f i n a l value, i n t e g r a t o r 2 holds that value with the relay open. B.3.1 Meissinger's Method; Implementation #1 The analogue computer implementation of t h i s method can be accomp-l i s h e d by using the r e p e t i t i v e operation mode on the analogue computer. For a system with n=3. , three system, three model, nine s e n s i t i v i t y equations and three update routines are required. UPDATE SEGMENT Figure (B.3 tn •p-55 Putting the update routines into the continuous mode and a l l the other integrators into the r e p e t i t i v e operation mode, we can implement the scheme i n Figures (B.l) - (B.3) provided we have a s u f f i c i e n t supply of components i n the computer. Since three m u l t i p l i e r s are required i n the model equation, nine m u l t i p l i e r s are required i n the update routine, and nine m u l t i p l i e r s are required i n the s e n s i t i v i t y model, we can not implement the method on the analogue computer as the analogue computer has only twelve m u l t i p l i e r s . B.3.2 Meissinger's Method; Implementation #2 Since the structure of the s e n s i t i v i t y model equations i s very s i m i l a r , as can be seen i n equation (3.7), we can cycle through three s e n s i t i v i t y equations with three d i f f e r e n t inputs to obtain the nine s e n s i t i v i t y equations. This scheme involves f i v e modes of operation, namely; one mode (M and M) for the system and the model, three modes (SI, S2, S3) for the s e n s i t i v i t y model,, and one mode (S4) f o r the update routine. Since the analogue computer has at most two modes of operation, i t i s necessary to design an external mode c o n t r o l l e r . This device consists of three f l i p f l o p s , four AND gates and four i n v e r t e r s , as seen i n Figure. B.5. The mode control signals previously described appear i n Figure B.6. The analogue computer implementation can be seen i n Figures B.7 and B.8. A l l integrators with no mode con t r o l i n d i c a t e d i n both the system and the model are i n the M mode. A l l the integrators with no mode co n t r o l ind i c a t e d i n the s e n s i t i v i t y model are i n the M mode. Relay #1 and the track and hold devices ensure that a,b,c, enter the c i r c u i t once during the cycle and at the correct time. FIGURE B.5 FIGURE B.6 s4 FIGURE B.7 co 60 Since the implementation of t h i s scheme requires 30 in t e g r a t o r s , the analogue computer could not be used, as only 26 integ r a t o r s are a v a i l a b l e on the. computer. B.A S t r u c t u r a l S e n s i t i v i t y Implementation Since the s e n s i t i v i t y t r a n s f e r functions are i d e n t i c a l for each input, input switching was considered. The same mode co n t r o l as the previous example can be used. Using t h i s scheme, the number of integrators needed exceed the number a v a i l a b l e , and so the computer could not be used. B.5 Use of the Hybrid Computer' The problem with a l l the methods c i t e d i s e i t h e r t h e i r lack of m u l t i p l i e r s or the lack of storage elements (track and hold devices i . e . i n t e g r a t o r s ) . Using the hybrid computer f o r m u l t i p l i c a t i o n and storage as w e l l as i n t e g r a t o r mode cont r o l would eliminate these d i f f i c u l t i e s . 61 / REFERENCES 1. Hsia, T.S., and V. Vimolvanich, "An On-line Technique for System / ' I d e n t i f i c a t i o n " , I n s t i t u t e of E l e c t r i c a l and E l e c t r o n i c Engineers Transactions on Automatic Control, February, 1969. pp. 92-96. 2. Gavrilov, M., and R. P e t r o v i c , "Parameter I d e n t i f i c a t i o n v i a System L i n e a r i z a t i o n and Optimal Improvement of Parameter Values", 5th I n t e r n a t i o n a l Analogue Computation Meetings, pp. 545-549, 1967. 3. Kohr, R., "On the I d e n t i f i c a t i o n of Linear and Nonlinear Systems", Simulation, March, 1967. . pp. 165-174. 4. Nagumo,J.., and Noda, A., "A Learning Method for System I d e n t i f i c a t i o n " , I n s t i t u t e of E l e c t r i c a l and E l e c t r o n i c Engineers Transactions on Automatic Controls, v o l . AC-12, No. 3, June, 1967, 282-287. 5. Vuskovic, M.I., J.V. Medanic, and W. Schaufelberger, "Investigation of a Clas.s of Fast Adaptive Control Systems", 5th International Analogue Computation Meetings, pp. 567-574, 1967. 6. B u t l e r , R.E., and E.V. Bohn, "An Automatic I d e n t i f i c a t i o n Technique for a Class of Nonlinear Systems,"' IEEE Transactions on Automatic Controls, vol.'AC-11 no. 2, A p r i l , 1966. pp. 292-296. 7. • Mishkin, E. and Braun, L. "Adaptive Control Systems", Toronto, McGraw-H i l l . 1961. 8. Kokotovic, P., "Method of S e n s i t i v i t y Points i n the Investigation and Optimization of Linear Control Systems", Automatika i Telemekhaneka, v o l . 25 no. 12, pp. 1670-1676, December, 1964. 9. Kokotovic, P. and R.S. Rutman, " S e n s i t i v i t y of Automatic Control Systems (Survey)", Automatika i Telemekhanika, v o l . 26, no. 4, pp. 730-750, A p r i l , 1965. 62 10. Tomovic, R. , " S e n s i t i v i t y Analysis of Dynamic Systems", McGraw-Hill, N.Y., 1963. 11. Meissinger, H.F., "The Use of Parameter Influence C o e f f i c i e n t s i n Computer Analysis of Dynamic Systems", Proceedings of the Western J o i n t Computer Conference, pp. 181-192, May, 1960. 12. Rosenbrock, H.H., "An Automatic Method for Finding the Greatest or Least Value of a Function", Computer Journal, v o l . #3, pp. 175-184, 1960. 13. Beveridge,G. and R.S. Schechter, "Optimization Theory and P r a c t i c e " , McGraw-Hill, N.Y., 1970 pp. 363-367. 14. Powell, M.J.D., "An I t e r a t i v e Method f o r Finding Stationary Values of a Function of Several Variables", Computer Journal, v o l . #5, pp. 147-151, 1955. 15. Curry, H.B., "The Method of Steepest Descent f o r Nonlinear Minimization Problems", Quarterly of Applied Mathematics v o l . ill, pp. 258-261, October, 1944. 16. F l e t c h e r , R. , and CM. Reeves, "Function Minimization by Conjugate Gradients", B r i t i s h Computer Journal, pp. 163-171, July, 1964. 17. Bingulac, S.P., "Parameter Optimization on an Analog Computer", 5th Int e r n a t i o n a l Analogue Computation Meetings, pp. 562-566, 1967. 18. C o l l i n s , P.L., and H.C. Khatre, " I d e n t i f i c a t i o n of D i s t r i b u t e d Para-meter Systems Using F i n i t e Differences", Transactions of the ASME Journal of Basic Engineering. June, 1969, pp. 239-245. 19-. Tzafestas, S.G. , " I d e n t i f i c a t i o n of Stochastic D i s t r i b u t e d Parameter Systems, International Journal of Control, v o l . 11, no. 4, pp. 619-624, 1970. 63 20. P e r d r e a u v i l l e , F.J. and R.E. Goodson, A.S.M.E. Journal of Basic Engineering, pp. 463- , 1966. 21. Shinbrot, M. , NACA Tech. Note IN3288, 1954. 22. Fairman, F.W., and D.W.C. Shen, "Parameter I d e n t i f i c a t i o n f o r a Class of D i s t r i b u t e d Systems", Int. J . Control, v o l . #11, no. 6, pp. 929-940, 1970. 23. Bensoussan, A., " I d e n t i f i c a t i o n de Systemes Gouverne's per des Equation aux Derive'es! P a r t i c l e s " , Computing Methods i n Optimization Problems, no. 2. ed. Zadeh, L.A., Neustadt, L.W., and Balakmshian, A.V. pp. 25-34, A.P. 1969. 24. Crank, J . , and P. Nicholson, "A P r a c t i c a l Method for Numerical Evaluation of Solutions of P a r t i a l D i f f e r e n t i a l Equations of the Heat Conducting Type", Proc. Cambridge Phi l o s S o c , v o l . #43, pp. 50-67, 1947. 25. Bekey, G.A., and W.Ji Karplus, Hybrid Computation,Wiley, ".N.Y., 1968. pp. 214-221. . • ' 26. Bingulac, S. and P. Kokotovic, "Automatic Optimization of Linear Feed-back Control Systems on an Analogue Computer", Proc Internat. Assoc. f o r Analogue Computation, v o l . #7, pp. 12-17, January, 1965. 27. Korn, G.A., "Random Process Simulation and Measurements", McGraw-Hill, N.Y. , 1966. 28. Tomkins, C.B., "Modern Matematics for the Engineer", E.F. Beckenbach, ed. McGraw-Hill, Toronto, 1956, pp. 448-479. 29. Forsyt.he, G.E. and Wasow, W.Z. " F i n i t e Difference Methods for P a r t i a l D i f f e r e n t i a l Equations", Wiley, N.Y., 1960, pp. 88-139. 30. G i l b e r t , E.G., "A Selected Bibiography, on Parameter Optimization Methods Suitable for Hybrid Computation", Simulation, v o l . #8, pp. 350-52, June 1967. 64 31. Vuscovic, M. and C i r i c , V., " S t r u c t u r a l Rules for the Determination of S e n s i t i v i t y Fuctions of Nonlinear Nonstationary Systems", S e n s i t i v i t y Methods i n Control Theory ed. Radanovic L., Pergamon Press, 1964, pp. 153-165.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Identification of parameters in distributed parameter...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Identification of parameters in distributed parameter systems Bordan, Norman Stephen 1971
pdf
Page Metadata
Item Metadata
Title | Identification of parameters in distributed parameter systems |
Creator |
Bordan, Norman Stephen |
Publisher | University of British Columbia |
Date Issued | 1971 |
Description | This thesis deals with the identification of parameters in distributed parameter systems. Two sensitivity methods, namely; Meissinger's method and the method of structural sensitivity are extended to obtain the sensitivity coefficients of discretized distributed parameter systems. The method of Bingulac and Kokotovic is extended to identify parameters in the one and the two dimensional parabolic differential equations. CSMP (continuous system modeling programme) is used throughout to simulate the systems. Results for both sensitivity schemes are obtained, and it is found that although structural sensitivity is advantageous for parameter identification in ordinary differential equations, this is not the case for partial differential equations. |
Subject |
System analysis |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-05-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0101895 |
URI | http://hdl.handle.net/2429/34226 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1971_A7 B67_3.pdf [ 2.78MB ]
- Metadata
- JSON: 831-1.0101895.json
- JSON-LD: 831-1.0101895-ld.json
- RDF/XML (Pretty): 831-1.0101895-rdf.xml
- RDF/JSON: 831-1.0101895-rdf.json
- Turtle: 831-1.0101895-turtle.txt
- N-Triples: 831-1.0101895-rdf-ntriples.txt
- Original Record: 831-1.0101895-source.json
- Full Text
- 831-1.0101895-fulltext.txt
- Citation
- 831-1.0101895.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
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-0101895/manifest