Optimal Design of Distributed Parameter Measurement Systems Citation Kumar, Sudarshan (1978) Optimal Design of Distributed Parameter Measurement Systems. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/g5zj-8n68. https://resolver.caltech.edu/CaltechTHESIS:01262026-234414750 Abstract The measurement system design for a process governed by partial differential equations (so called distributed parameter system) consists of the choice of state variables to be measured, the frequency of measurements, and the optimal location of sensors. In this work, we consider the problem of optimal location of sensors in distributed parameter systems. The optimal sensor location problem is approached from the point of view of estimation theory. The objective chosen is the estimation of the system state variables with maximum possible accuracy. An upper bound for the error covariance of the state estimate is derived. Optimal sensor locations are determined by minimizing this upper bound of the error covariance. This approach offers significant savings in computation time over the conventional approaches based on direct minimization of the error covariance. A one-dimensional heat conduction system is considered as an example. The effect of process dynamics on optimal sensor locations is studied in detail for that example. The optimal locations are found to vary for spatially varying measurement noise. To apply the theory to a technologically important system we consider the optimal location of measurements for tubular reactors. The state of this system is described by two variables - concentration and temperature-and is governed by two nonlinear partial differential equations. For the first time, the optimal sensor location problem for a system with two state variables has been studied. For this problem, an approach based on discretization of the system equations is developed. The feasibility of estimating both concentration and temperature profiles from temperature measurements only is demonstrated. The effect of process noise and measurement noise on optimal sensor locations was studied and it was found that optimal sensor locations move towards region of lower measurement error variance. For the particular set of parameters chosen, it was found that two optimally located temperature sensors can give a sufficiently accurate estimate of the concentration and temperature profiles. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Chemical Engineering) Degree Grantor: California Institute of Technology Division: Chemistry and Chemical Engineering Major Option: Chemical Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Seinfeld, John H. Thesis Committee: Unknown, Unknown Defense Date: 1 September 1977 Record Number: CaltechTHESIS:01262026-234414750 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:01262026-234414750 DOI: 10.7907/g5zj-8n68 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 17839 Collection: CaltechTHESIS Deposited By: Benjamin Perez Deposited On: 28 Jan 2026 17:42 Last Modified: 28 Jan 2026 17:44 Thesis Files PDF - Final Version See Usage Policy. 29MB Repository Staff Only: item control page

OPTIMAL DESIGN OF DISTRIBUTED PARAMETER MEASUREMENT SYSTEMS Thesis by Sudarshan Kumar In Partial Fulfillment of the Requirements for the degree of Doctor of Philosophy California Institute of Technology Pasadena, California 1978 (Submitted September 1, 1977) ii This thesis is dedicated to my parents. iii ACKNOWLEDGEMENTS I wish to thank my advisor, Professor J. H. Seinfeld, for his guidance and encouragement throughout the course of this work. I sincerely appreciate his help and advice on various technical matters related to the work presented here. I would also like to thank Ken Bencala and Tom Peterson, with whom I shared an office, for their friendship and various discussions regarding this work. The discussions with Deepak Dhar and Greg McRae have also been most invaluable. I would also like to thank Lenore Kerner, who typed most of this thesis, for her superb typing and cheerful help. I would like to express my deep gratitude to my brother Ramesh who has gladly shouldered so many responsibilities that were properly mine. Finally, the support and patience of my wife Anjali has been particularlY,- helpful. I appreciate her ·friendship very much. iv ABSTRACT - The measurement system design for a process governed by partial differential equations (so called distributed parameter system) consists of the choice of state variables to be measured, the frequency of measurements, and the optimal location of sensors. In this work, we consider the problem of optimal location of sensors in distributed parameter systems. The optimal sensor location problem is approached from the point of view of estimation theory. The objective chosen is the estimation of the system state variables with maximum possible accuracy. An upper bound for the error covariance of the state estimate is derived. Opti- mal sensor locations are determined by minimizing this upper bound of the error covariance. This approach offers significant savings in com- putation time over the conventional approaches based on direct minimizaft tion of the error covariance. A one-dimensional heat conduction system is considered as an example. The effect of process dynamics on optimal sensor locations is studied in detail for that example. The optimal lo- cations are found to vary for spatially varying measurement noise. To apply the theory to a technologically important system we consider the optimal location of measurements for tubular reactors . The state of this system is described by two variables - concentration and temperature-and is governed by two nonlinear partial differential equations. For the first time, the optimal sensor location problem for a system with two state variables has been studied. For this problem, an approach based on discretization of the system equations is developed. V The feasibility of estimating both concentration and temperature profiles from temperature measurements only is demonstrated. The effect of process noise and measurement noise on optimal sensor locations was studied and it was found that optimal sensor locations move towards region of lower measurement error variance. For the particular set of parameters chosen, it was found that two optimally located temperature sensors can give a sufficiently accurate estimate of the concentration and temperature profiles. vi TABLE OF CONTENTS ACKNOWLEDGEMENTS iii ABSTRACT iv LIST OF FIGURES viii LIST OF TABLES X CHAPTER l: INTRODUCTION l CHAPTER 2: OPTIMAL LOCATION OF MEASUREMENTS FOR DISTRIBUTED PARAMETER ESTIMATION ll 2. l Introduction 12 2.2 System Description 15 2.3 Optimal Filters for Linear Distributed Systems 18 2.4 Integral Form for the Linear Dynamical System 26 2.5 Optimality Index 29 2.6 Upper Bound for the Covariance Matrix 2.7 Approximation of Pin Finite-Dimensional Ir Subspace 2.8 30 33 Optimal Location of Sensors in a One-Dimensional Diffusion System 38 2.8. l Effect of Boundary Conditions 40 2.8.2 Effect of R(x) 43 2.8.3 Effect of P0 (x,y) 48 2.8.4 Comparison of the Results to Solution of The Riccati Equation 2.9 Conclusions 48 51 vii TABLE OF CONTENTS (Continued) CHAPTER 3: OPTIMAL LOCATION OF MEASUREMENTS IN TUBULAR REACTORS 53 3. l Introduction 54 3.2 Reactor Model 56 3.3 Discretization by Orthogonal Collocation 62 3. 4 State Space Formulation and Optimal Estimate 69 3.5 Optimality Index and the Minimum Principle 75 3.6 Numerical Example and Results 78 3.6. l Steady State Solution 79 3.6.2 Determination of Optimal Locations 83 3.6.3 Effect of Varying Measurement Noise 86 3.6.4 Effect of System Noise 91 3.6{5 Number of Sensors 95 3.6.6 Accuracy of Concentration and Temperature Estimates 3.7 Surrrnary CHAPTER 4: SUMMARY AND RECOMMENDATIONS l 01 107 108 REFERENCES 112 APPENDIX A 118 APPENDIX B 122 APPENDIX C 125 viii LIST OF FIGURES Figure l. l Page Block diagram description of the general problem of estimation and control of dynamic systems. 1.2 3 Approaches to distributed parameter estimation problems. 8 2. l Contour plot of J for R(x) = 0.8 exp(-2x) 45 2.2 Contour plot of J for R(x) = 0.8 [-2(1-x)] 46 2.3 Contour plot of J for R(x) = 0.8 exp[+4lx-0.5I] 47 3. l Steady state concentration as a function of s. 81 3.2 Steady state temperature as a function of s. 82 3.3 Trace Pas a function of time for two sensors located at (A) positions s 4 and s 6 (B) optimal positions s 2 and s 3; R(s) = 0.0004 88 ~ 3.4 Trace Pas a function of time for two sensors located at (A) positions s 4 and s 6 (B) optimal positions s 2 and s 3; R(s) = 0.004 3.5 89 Trace Pas a function of time for two sensors located at (A) positions s 3 and s 6 (B) optimal positions s 3 and s 4 ; R(s) = 0.0001 exp(-lOls-0.51) 3.6 90 Trace Pas a function of time for two sensors located at (A) positions s 4 and s 6 (B) optimal positions s 2 and s 3 ; Q = 0. l I 3.7 93 Trace Pas a function of time for two sensors located at (A) positions s 4 and s 6 (B) optimal positions s 2 and s ; Q+(s,s') = 0.01 exp(-lOls-s'I) 3 94 ix LIST OF FIGURES (Continued) Figure 3.8 Trace Pas a function of time for one sensor located at (A) position s8 (B) optimal position 97 s2 3.9 Trace Pas a function of time for three sensors located at (A) positions s 4, s 6 , s 7 (B) optimal positions s 2, s 3 , s4 3.10 98 Trace Pas a function of time for five sensors located at (A) positions s 4, s 5, s , s 7 , s8 6 (B) optimal positions s 2, s 3 , s 4 , s 5, s 6 3. 11 Trace Pas a function of time for (A) one (C) three (D) five, and (E) 99 (B) two seven optimally 100 located sensors -~ 3. 12 (A) Trace P (B) trace Pc, and (C) trace PT as functions of time for one optimally located sensor 3. 13 (A) Trace P (B) trace Pc, and 102 (C) trace PT as functions of time for two optimally located sensors 103 1 3.14 (A) Trace P (B) trace Pc, and (C) trace PT as functions of time for three optimally located sensors 3. 15 (A) Trace P (B) trace Pc, and (C) trace PT as functions of time for five optimally located sensors 3. 16 (A) Trace P (B) trace Pc, and 104 105 (C) trace PT as functions of time for seven sensors. 106 X LIST OF TABLES Table Page 2. 1 Optimal filters for various measurement models. 2.2 Eigenfunctions and Green's functions for different boundary conditions. 2.3 49 Values of J(t) obtained by solving the full Riccati equation for various combinations of sensor locations. 3. l 44 Effect of initial covariance P0 (x,y) on optimal sensor locations. 2.6 42 Effect of measurement noise covariance R(x) on optimal sensor locations. 2.5 41 Optimal sensor locations for different boundary conditions. 2.4 23 50 Orthogonal polynomials and the interior collocation poirhs. 80 3.2 Effect of measurement noise. 87 3.3 Effect of process noise. 92 3.4 Optimal locations for one, two, three, and five sensors. 96 CHAPTER 1 INTRODUCTION 2 The general problem of estimation and control of dynamic systems can be described with the help of Figure 1. 1. In practice, it is often desired that the system state follow a particular path, and control action is necessary to guide the system along this specified path. Before one can decide on the required control action, it is necessary that adequate knowledge of the system state be available. Knowledge of the system state can be acquired through comparison of the predictions of a system model and measurements on the system. Since most physical systems are subject to random disturbances and uncertainties in system parameters, deterministic models are often inadequate. A common approach to account for the stochastic nature of a system is to include a random dynamical disturbance in the system equations. In addition measurements of the system state contain errors that can be described as a random process. The process of obtaining an estimate of the system state from noisy - known as filtering and the algorithm used for obtaining measurements~ 1s the estimate is called a filter. If a performance measure is introduced into the estimation problem, we can define an optimal estimate obtained from an optimal filter. Mathematical models of physical systems can be classified as lumped parameter or distributed parameter. Systems described by ordinary differential equations are termed lumped parameter systems. (These systems are also known as finite-dimensional.) On the other hand, systems whose state varies spatially as well as temporally are described by partial differential equations or integral equations and are called distributed parameter systems. dimensional.) (These systems are known as infinite- Desired State State Estimate Control Inputs Estimator or Filter Physical System ' ' Svstem Variables , Statistica Information on Errors and Disturbances Measurement System Measurement Errors FIGURE 1.1: BLOCK DIAGRAM DESCRIPTION OF THE GENERAL PROBLEM OF ESTIMATION AND CONTROL OF DYNAMIC SYSTEMS Controller . Random Input Disturbances ., w 4 Kalman and Bucy (1961) developed the optimal filter for linear lumped parameter systems. The state estimate depends on the system dynamics and the measurements. The estimate error is measured by the error covariance which is governed by a Riccati differential equation. Balakrishnan and Lions (1967), Falb (1967) and Kushner (1970} extended the Kalman-Bucy filter for linear distributed systems. Tzafestas and Nightingale (1968a, 1968b, 1969), Thau (1969), Meditch (1970), and Sakawa (1972) obtained filters for linear distributed systems using various approaches. Estima- tion of states in nonlinear distributed systems has been considered by Seinfeld et al. (1971} and Hwang et al. (1972). Yu et al. (1974} also derived filters for nonlinear time-delay systems. The accuracy of the state estimate obtained from optimal filters depends on the measurements available on the system. The measurement sys- tem design problem in lumped system is concerned basically with the question of whif h states are to be measured. On the other hand for distribu- ted systems, the measurement system design problem is more complicated. Since it is physically impossible to measure the state over the entire spatial domain, the question of just where to locate these sensors is an important one. In addition, to specify the measurement system design com- pletely, the questions of choice of states to be measured, frequency of measurements, and number of measurement sensors should also be dealt with. In most physical processes, however, there will normally be little choice in the states to be measured because of the difficulty in continuously measuring certain state variables. The number of sensors is usually governed by economic considerations. Assuming that the state variables 5 being measured can be measured continuously, the major question in measurement system design is the selection of sensor locations. Let us discuss briefly two engineering systems of interest to be considered in this work. The first system is a bar undergoing heat con- duction in which heat is lost to the surroundings along the length of the bar. The unsteady state temperature u is described by (Section 2.8 defines the nomenclature) We assume that this idealized model does not represent the physical process exactly. The temperature profile in the bar is to be estimated by measuring the temperature at certain points. Recognizing that the tempera- ture measurements are not necessarily perfect, we wish to select the location of temperature sensors so that a "best" estimate of the system state J! (i.e. the temperature) is obtained. • Some interesting questions are: Should the sensors be located near the ends of the bar or near the middle? t Do the locations change if the boundary conditions on the system change? • How do the measurement locations change if the measurement errors and model errors are functions of space? The second system to be considered is a non-adiabatic tubular reactor in which an exothermic chemical reaction is taking place. In this system, there are two state variables - concentration and temperature. The equations defining the system are nonlinear and coupled (and will be given 6 in Chapter 3). Assuming the model and the measurements to be imperfect, we wish to estimate the concentration and temperature profiles in the reactor. Since we can make measurements at only a finite number of points along the reactor length, the choice of sensor locations is very important. In general, it is difficult to make continuous concentration meas- urements and therefore we would like to explore the possibility of estimating both concentration and temperature profiles based on temperature measurements only. Measurement system design encompasses all these and other related questions. The optimal design of lumped parameter measurement systems has received considerable attention although comparatively little work has been directed towards the specific problem of measurement locations in distributed parameter systems. For distributed parameter systems Bensoussan (1972) derived necessary conditions for optimality of sensor locations. He posed the problem ~ s an optimal control problem with the state described by the infinite dimensional Riccati equation and the sensor locations considered as control variables. Based on a modal (eigenfunction expansion) approach, Yu and Seinfeld (1973) developed an algorithm for determining suboptimal sensor locations. Brewer and Moore (1974) considered the use of estima- tion theory for minimizing the cost of measurements while meeting certain estimate accuracy criteria. Pimentel (1975) considered the infrequent sampling problem and employed a modal representation of the distributed system to determine optimal sensor locations. Chen and Seinfeld (1975) derived necessary conditions for optimal sensor locations for a general linear distributed system and developed an algorithm for determining the set of optimal locations. Aidarous (1976) considered discrete-time 7 measurements and determined sensor locations such that a scalar measure of the error covariance matrix over one time increment is minimized. The problem of optimal sensor locations, approached as an estimation problem, involves a finite dimensional approximation at some point. are basically two approaches. There The distributed parameter system can be approximated right at the beginning, and the sensor location problem is attempted by employing finite-dimensional filters. called approximation at the beginning. This approach can be On the other hand, one can retain the distributed nature of the problem by employing infinite-dimensional filters and use finite-dimensional approximation for numerical implementation. This approach can be termed approximation at the end. Of course, the techniques for finite-dimensional approximation in both the cases are the same. Figure 1.2 illustrates the concepts of approximation at the beginning and approximation at the end. Various techniques that can be used for fin1te-dimensional approximation are finite-difference methods, orthogonal collocation, Galerkin method and eigenfunction expansion. Sein- feld and Koda (1977) have discussed these techniques for numerical implementation of distributed parameter filters. The basic problem in determining optimal sensor locations is the computational requirement. In all the approaches available, it has been necessary to solve the infinite-dimensional Riccati equation. In this work, our objective is to develop computationally attractive approaches for the optimal sensor location problem. In Chapter 2, we consider a general stoch- astic linear distributed parameter system and develop an upper bound for the estimate error covariance matrix. The optimal sensor locations are 8 DISTRIBUTED PARAMETER SYSTEM (PDE) ~, ~, DISTRIBUTED PARAMETER FILTER FINITE-DIMENSIONAL APPROXIMATION Finite Differences Eigenfunction Expansion -Methods of Weighted Residuals ~· LUMPED PARAMETER FILTER ~ ... APPROXIMATION AT THE BEGINNING Integro-Partial Differential Equations ~ FINITE-DIMENSIONAL APPROXIMATION Finite Differences, etc. STATE (AND PARAMETER) ESTIMATES ~ , APPROXIMATION AT THE END FIGURE 1.2: APPROACHES TO DISTRIBUTED PARAMETER ESTIMATION PROBLEMS 9 determined by minimizing a scalar measure of this upper bound. The prob- lem requires the solution of a linear partial differential equation rather than the nonlinear Riccati partial differential equation. The linear equation is solved by employing the system Green's function. approach falls into the category approximation at the end. This The optimal sensor locations are obtained by minimizing a function of the variables defining the sensor locations. An example involving one dimensional diffusion system is considered. The effect of process dynamics and vary- ing measurement noise covariance on optimal sensor locations is studied. The results are compared to those obtained from the conventional approach of minimizing a measure of the filter covariance matrix. In Chapter 3, we consider the optimal sensor location problem for non-adiabatic tubular reactors. The system state is described by two variables, namely concentration and temperature. This is apparently the first tiine that an optimal sensor location problem for a system with more than one state variable has been attempted. An alternate approach involving approximation at the beginning is developed. The system equa- tions are discretized by orthogonal collocation and the resulting nonlinear equations are linearized around the steady state. The finite- dimensional optimal filter is used, and the optimal sensor locations are determined by employing the matrix minimum principle (Athans, 1968) and an algorithm similar to the one used by Chen and Seinfeld (1975). The effect of process dynamics on the optimal sensor locations is studied. Variations in the optimal sensor locations with different process and measurement noise statistics is also studied. The feasibility of estimat- ing both concentration and temperature profiles based on only temperature 10 measurements is clearly demonstrated. The question of number of sensors to estimate both concentration and temperature profiles is also investigated. A summary of the work done and recominendations for further work are given in Chapter 4. CHAPTER 2 OPTIMAL LOCATION OF MEASUREMENTS FOR DISTRIBUTED PARAMETER ESTIMATION 12 The general optimal measurement design problem for linear systems can be approached with the objective of obtaining the best possible state estimates, in which case the problem becomes one of minimizing a functional depending on the covariance of the optimal filter. The specific problem of optimal measurement locations for linear distributed parameter systems is considered. Bensoussan (1972) and Chen and Seinfeld (1975) have shown that the optimal sensor location problem for distributed systems can be posed as an optimal control problem for a system described by the infinite-dimensional matrix Riccati equation for the filter covariance. In this chapter, a more efficient approach based on an upper bound of the filter covariance is developed. The optimal sensor locations are determined by minimizing a nonlinear function of the variables defining the sensor locations. The relationship between the~present approach and that of minimizing a measure of the filter covariance is studied. A detailed example is considered, and the results of the two approaches are compared. 2.1 Introduction Measurement system design for a physical process generally in- volves the choice of the states to be measured and the frequency and location of measurements. In this chapter, we consider the problem of the detennination of optimal sensor locations for linear distributed parameter systems with continuous time measurements. Measurements at 13 certain points in the spatial domain of the system may yield more information about the system than the measurements at other points and, therefore, the accuracy of the state estimate depends on the number and locations of sensors. Since the number of sensors is generally governed by economic considerations, it is desirable to locate the given number of measurement sensors at points that lead to a best estimate of the system state. These optimal sensor locations are a function of the process dynamics and boundary conditions under consideration, and of the spatial dependence of measurement and system model uncertainties. The problem of dynamic measurement system design has been of particular interest. Johnson (1969), Muller and Weber (1972), and Mehra (1976) considered the problem for lumped parameter systems in which various structural parameters in the observability matrix are chosen so that a giv en measure of the observability matrix is optimized. Athans (1972) considered the problem of selecting one measurement provided by one out of many sensors, at each instant of time, during the observation. Herring and Melsa (1974) generalized Athans' results to selecting the best combination of measurements at each instant of time. Bensoussan (1972) considered the problem of optimal sensor loca- tion in distributed parameter systems. He posed an optimal control problem for the system described by the infinite-dimensional matrix Riccati equation governing the filter covariance with the sensor locations considered as control variables. Yu and Seinfeld (1973) con- sidered a modal approach to a general linear distributed parameter 14 system and developed an algorithm for determining suboptimal sensor locations. Following Bensoussan (1972), Chen and Seinfeld (1975) derived necessary conditions for optimality based on the Riccati equation governing the estimate error covariance, and developed an algorithm for determining the optimal sensor locations. Aidarous (1976) considered a discrete-time measurement distributed parameter filter and determined sensor locations such that a scalar measure of the error variance over one time increment is minimized. The dual problem of optimal pointwise control for distributed systems has also received attention (Amouroux and Babary, 1975). In previous work (Aidarous, 1976; Chen and Seinfeld, 1975), it has been necessary to solve nonlinear matrix Riccati differential equations to determine optimal sensor locations. In this work, we develop an approach which requires the solution of a linear matrix differential equation, thus making it computationally more attractive. We derive an upper bound for the error covariance matrix of a general stochastic linear distributed parameter system. The upper bound is a sum of two terms, only one of which is a function of the sensor locations. The optimal locations are determined by minimizing a scalar measure of this term with respect to sensor locations. For lumped parameter systems, Mehra (1976) has used a similar approach to select measurement schedules and sensor designs. In Section 2.2, we define the general linear distributed parameter system under consideration. Section 2.3 describes the general forms of the optimal distributed parameter filter. In Section 2.4, 15 we develop the integral fonn for the linear dynamical system described in Section 2.2. Section 2.5 is concerned with the optimality index for defining the optimal sensor locations. In Section 2.6, we develop an upper bound for the error covariance matrix P, involving the error covariance matrix P for the optimal filter for the system in the absence of dynamic noise. In order to work with the solution of P, in Section 2.7 we develop an approximation to the matrix ff in a finitedimensional subspace of the original infinite-dimensional space. When expressed in the final fonn, it is possible to obtain the optimal sensor locations by optimizing a nonlinear function of the variables defining these sensor loc~tions. In Section 2.8, we consider an ex- ample involving one-dimensional heat transfer. The effect of boundary conditions, varying measurement noise, and varying initial covariance on the optimal sensor locations is studied. The full nonlinear Riccati equ ation for the covariance matrix is also solved, and it is shown that the sensor locations obtained by the method developed here are approximately the same as those based on the full nonlinear Riccati equation. . 2.2 System Description In this work, we consider linear systems which can be described by an equation of the fonn au(x,t) = A u(x,t) + B(x,t) w(x,t) at X • (2. 1) on a connected open domain D of ad-dimensional Euclidean space Ed, with 16 boundary ao. by x. The d-dimensional spatial coordinate vector is denoted The state u(x,t) is an n-dimensional vector, and AX is an n x n matrix linear, spatial differential, or integro-differential operator. It is assumed that the operator Ax is well posed. B(x,t) is an n x s dimensional matrix and the process noise w(x,t) is ans-dimensional zero-mean white Gaussian process with covariance E[w(x,t)wT(y,T)] = Q(x,y,t)o(t-T) (2.2) where Q(x,y,t) is a synmetric positive semi-definite matrix. The initial and boundary conditions for the system are E[u(x,t 0 )J= u0 (x) (2.3) E[{u(x,t 0 )-u 0 (x)} {u(y,t0 ) - u (y)}T] = P (x,y) 0 and 0 Jt ax u(x,t) = 0 xE ao (2.4) The initial state u(x,t 0 ) is an n-dimensional white Gaussian process with mean u0 (x) and covariance matrix P0 (x,y). ax is an nxn matrix, linear, spatial differential operator of suitable order over ao. The boundary conditions have been assumed to be homogeneous, although the case of inhomogeneous boundary conditions poses no difficulties. It is well known (Courant and Hilbert, 1966) that a homogeneous differential equation with nonhomogeneous boundary conditions is essentially equivalent to a nonhomogeneous differential equation with homogeneous boundary 17 conditions. Brogan (1968) considers several examples and establishes in detail the above equivalence through a general approach. For con- venience we may express the system equation (2. 1) as the stochastic evolution equation on a real Hilbert space v d~~t) = A u(t) + B(t) w(t) (2.5) where A is a linear closed operator on V , and is the infinitesimal generator of a strongly continuous semigroup. B(t) denotes a bounded linear transfonnation, and w(t) denotes a white noise process. The system output, in general, is the result of a continuous transformation of the state functions. It is important to note that in a distributed parameter system, the measured variable can either be spatially dependent or spatially independent. The most important case is that in which measurements are made at a finite number of points in the spatial domain. ~ It is possible to represent fonnally both of these cases by representing the measurements as z(t) = H(t) u(t) + v(t) (2.6) where z(t) is an r-dimensional vector and H(t) is a bounded linear operator that may transfonn the state into spatially dependent or independent outputs. The measurement noise v(t) is assumed to be zero-mean white Gaussian process with covariance R(t). We note that the measure- ment noise covariance R(t) is spatially dependent or independent based on the type of ·measurements. We wish to obtain an estimate of the state u(t) based on all the past measurements z(t') 0 .2.. t' .2.. t. An 18 estimate u(t) which minimizes (2. 7) is called the optimate estimate. The estimate error u(x,t) = u(x,t) - u(x,t) (2.8) has the covariance P(x,y,t) defined by ~ ~T (y,t)J P(x,y,t) = E[u(x,t)u (2.9) The covariance matrix P(x,y,t) is symmetric and positive semi-definite. 2.3 Optimal Filters for Linear Distributed Systems Let us define the operators P(t) and Q(t) to be integral operators with kernels P(x,y,t) and Q(x,y,t) respectively. The equation ,ft . defining the optimal estimate is (Balakrishnan, 1976) d~it) = A u(t) + P(t)H*(t)R- 1(t) [z(t) - H(t) u(t)J (2.10) with the i ni ti al con di ti on (2. 11) The kernel matrix of the operator P(t)H*(t)R- 1(t) is called the gain matrix for the filter. The covariance matrix P(x,y,t) obeys the kernel form of the operator equation P(t) = AP(t) +P(t)A*+B(t)Q(t)B*(t)-P{t)H*(t)R-l(t)H(t)P(t) (2.12) 19 where* denotes the adjoint operator. The operator P(t) in (2.12) is a bounded integral operator with kernel P(x,y,t), and is defined as P(t)f = l P(x,y,t) f(y) dY (2.13) for fC L2(D;Rn), where L2 (D;Rn) denotes the vector space of all realvalued functions, square-integrable over the domain D. The covariance operator is trace class and nonnegative definite (Balakrishnan, 1976). Symbolically, we express the nonnegative definiteness as P(t) .2:.. 0, which implies that for any f C L 2 (D;Rn), ff fT(x)P(x,y,t)f(y) dydx ~ o (2. 14) DD The initial conditions and the boundary conditions for the covariance P(x,y,t) ar§ (2.15) ax P(x,y,t) = O xeao yE DxaD (2. 16) P(x,y,t)a; = 0 Xe DxaD where the transpose of a matrix differential operator ay[•] is defined T and is an operator to the left. by the relation (ayP) T = PT ay, example, if a[·] y then For 20 where a 0 , a 1 , and a 2 are assumed to be square matrices. For spatially dependent measurements, the measurement noise covariance is given by E[v(x,t)yT(y,T)] = R(x,y,t) o(t-T) (2. 17) while for spatially independent measurements, we have (2.18) The operator R(t) occurring in (2.10) and (2. 12) is, therefore, an integral operator with kernel R(x,y,t) in the first case, and simply the matrix R(t) in the second case. For spatially dependent measurements, therefore, we have R( ( )( • ) = J R( x,y, t )( • ) dy (2.19) • D The inverse operator R- 1(t) is an integral operator, with the property R- 1(t)[R{t)f] = f for all functions f in the domain of R( t). (2.20) In expanded form, this gives J Rt(x,y,t) J R(y,z,t) f(z) dzdy = f(x) D (2.21) D where Rt(x,y,t) is the kernel of the integral operator R- 1(t). From equation (2.21), it follows then, that Rt (x,y,t), the generalized 21 inverse of R(x,y,t), must obey the integral equation ! Rt(x,y,t) R(y,z,t) dy = I O(x-z) (2.22) R- 1(t) is thus an i~tegral operator with kernel Rt(x,y,t) satisfying the integral equation (2.22). For spatially independent measurements, the inverse operator R- 1 (t) is the ordinary inverse of matrix R(t). We now proceed to obtain the explicit equations for the optimal estimate u(x,t) and the estimate error covariance matrix P(x,y,t) for various kinds of measurement systems. The measurement systems are specifi~d for four cases and corresponding optimal filters are presented in Table2.l. ror all the cases, the initial condition and boundary conditions on the optimal estimate u(x,t) are (2.23) x € ao (2.24) The covariance matrix P(x,y,t) has the initial condition (2 . 25) and boundary conditions ax P( x ,y, t) = 0 x € ao , y € oxao (2.26) T p ( X ,y , t) ay = 0 y € ao , x € Dxao 22 Case 1: Continuous Measurements over the Entire Spatial Domain The measurement process is defined by z(x,t) = H(x,t) u(x,t) + v(x,t) (2.27} where v(x,t) is a white Gaussian process with properties E[v(x,t)J = a E[v(x,t)vT(y,t)J = R(x,y,t) o(t-T) (2.28) Therefore, the operator R( t) in the optimal filter is given by R(t)(•) = J R( x ,y, t }( • ) dy D R- 1(t)(•) = ( 2. 29) f Rt (x,y, t }( •) dy (2.30) D where Rt{x,y: t) is defined by the relation (2.22). Case 2: Continuous Integral Measurements over the Entire Spatial Domain The measurement process is defined by Z ( X, t) = J H( ,y, t) U(y, t) dy + X V ( X, t) (2.31) D where the measurement noise v(x,t) has the same properties as in Case 1. Consequently R(t) and R- 1(t) are also defined as in equations (2.29) and (2. 30). DOD III P(x,x 1 ,t)HT(x',x 11 ,t)Rt(x 11 ,y 11 ,t) [z(y 11 ,t) - DODD 11 11 H(y 11 ,y 1 ,t)u(y 1 ,t)dy'] dy 11 dx 11 dx 1 D {I 11 11 1 -[l 1 P(x,x' ,t)H (x' ,t)dxJ R-l (t{z(t) - 11 D I H(y' ,t)G(y' ,t) dy] aP(x,y,t) -- Ax P( x,y,t ) + P( x,y,t ) AYT + B( x,t )Q( x,y,t ) BT( y,t ) at 1 1 P(x,x',t)H (x',t) dx]R- ( t { [ H(y',t)P(y',y,t) dy] aui~·t) = Ax U(x,t) Case 3 D I •P(y' ,y,t) dy'dx' - fJf f P(x,x' ,t)HT(x' ,x" ,t)Rt(x ,y ,t)H(y ,y ,t)P(y' ,y,t) dy'dy dx dx' aP(~ty,t) = Ax P(x,y,t) + P(x,y,t)A~ + B(x,t)Q(x,y,t)BT(y,t) au~~,t) = Ax u(x,t) + Case 2 DD aP(~ty,t) = Ax P(x,y,t) + P(x,y,t)A~ + B(x,t)Q(x,y,t)BT{y,t) - ff P(x,x' ,t)HT(x' ,t)Rt(x' ,y' ,t)H(y 1 ,t) DD au~~,t) = Ax ~(x,t) _ + ff P(x,x 1 ,t)HT(x 1 ,t)R~Jx 1 ,y' ,t) [z(y' ,t) - H(y' ,t)~(y' ,t)] dy'dx' Case 1 Table 2. 1 Optimal Filters for Various Measurement Models N w m P(x,n , t)H1 (t) [ R- l (tJij [z/ t) - Hj ( t)C(xj, t ] 1 ~~ (continued) l,J=l . _I P ( x, n; , t) H~ ( t) [ R- l ( t !_J,J ~ . . FfJ. ( t) P ( nJ. ,y, t) aP(x,y,t) = A P(x y,t) + P(x,y,t)AT + B(x,t)Q(x,y,t)BT(y,t) at X ' Y au(x,t) = A ~(x,t) + . 2 at X , ,j=l Case 4 Table 2.1 ~ N 25 Case 3: Spatially-Independent Integral Measurements The measurement equation takes the fonn z(t) = f H(x,t) u(x,t) dx + v(t) (2.32) D where v(t) is a white Gaussian process with properties E[v(t)] = 0 (2.33) The operator R{t) in this case is the matrix R(t) and R- 1(t) is the ordinary inverse of matrix R(t). Case 4: Continuous Measurements at m Discrete Points in Space The measurements are represented as z.(t) = H.(t) u(x.,t) + v.(t) l l l l i=l,2,···,m (2.34) where vi(t) is a white Gaussian process with properties i = 1,2,··•,m E[v.(t)J =0 l . (2.35) T E[v.(t) v.(T)J = c.. (t) o(t-T) l J lJ i .s_ i ,j .s_ m The measurements given by (2.34) can be represented in the form (2.32) by considering the mr-dimensional vectors z(t) and v(t) defined by T(t)QT z(t) = zT (t) zT (t) • • • • • zm 2 [ 1 and (2. 36) 26 (2. 37) The kernel matrix H(x,t) in equation (30) takes the form (2.38) The operator R(t) is the mrx mr matrix R(t) with Ci/t) as the (i ,j)th In the filter equations [R- 1(t)]ij is the (i ,j)th submatrix submatrix. of dimensions rx r in the matrix R- 1 (t). It may be pointed out that Cij(t), the covariance mat r ix of t he measurement errors at locations xi and xj may be obtained by evaluating a given function C(x,x ,t) at 1 x=x ,. and x =x J.. The matrix function C.( x,x ,t) specifies the measure - 1 1 ment error covariance as a continuous function of any two points x and x 1 in the spatial domain D. 2.4 Integr~l Form for the Linear Dynamical System The general solution of equation (2.1) with homogeneous boundary conditions can be expressed (Courant and Hilbert, 1966) in the form u(x,t) = G(x,t;x ,t0 ) u(x ,t 0 ) dx 1 1 f 1 D t + G(x,t;x ,t 1 f f 1 B(x , t 1 ) 1 w(x ,t 1 ) 1 dx dt 1 ) 1 (2.39) t0 D where G(x,t; x ,t 1 1 ) is the Green's function matrix, and satisfies the equation aG(x,t ·,x 1 ,t 1 -) = A G(x t ·x t at X ' ' 1 1 ) (2. 40) 27 with the terminal condition G(x,t';x',t') = I o(x-x') (2.41) and boundary condition x e ao ax G(x,t;x' ,t') = 0 (2.42) If the operator Ax is time and space-invariant, it is sufficient to consider G(x,t;O,O) and then use the relation G(x,t;x',t') = G(x-x',t-t';0,0) (2.43) For convenience, we can express (24) as t u(t) = ¢(t,t0 )u(t0 ) + f ¢(t,t')B(t' )w(t') dt' (2.44) to where ¢{t,t ) and ¢(t,t') define linear transfonnations, the domain of 9- which is the set of bounded functions of x and the range of which is a subset of the domain. The transformation ¢{t,t') is defined by the relation [~(t,t')f] = l G(x,t;x',t')f{x') dx' for any f ~ L2 (D;R"), and obeys the relations (2.45) lim ¢(t,t') = I t-+t I (2.46) The equivalent concept of a transition matrix in lumped systems has the important property of invertibility. However, in distributed parameter 28 systems, the operator ¢(t,t 1 ) may or may not be invertible, depending on the system under consideration. The transformations that have the property of invertibility are termed group operators, whereas the transformations lacking this property are ca 11 ed semi -group operators. Parabolic systems, for example the diffusion equation, give rise to semi-group operators. (As a consequence of this prope'.ty, solving the heat equation backwards in time is not a well-posed problem). However, if the operator ¢(t,t 1 ) has the properties of a group, the complete past history of the system can be determined. For example, in wave processes, a complete knowledge of the state at a given tiJTE t, and of the inputs applied to the system is sufficient to reconstruct all future and past states. ~1e now list some important properties of semigroup operators (Hille and Phillips, 1957; Yosida, 1966). (1) ¢{t,t 1 ) is a bounded operator J! (2) lim ¢{t,t t+ t (3) 1 ) 1 =I ¢(t,t )¢{t ,t 11 (4) ¢(t,t 1 11 ) ) = ¢(t,t (2. 47) 1 ) > tll > t I ~ 0 t is strongly continuous, i.e., lim 1i(¢(t+o,t 0~ 1 ¢(t,t ))u(t)II 1 1 ) - • = 0 for all u in the domain An abstract approach for the development of Green's functions has been developed (Hille and Phillips, 1957; Yosida, 1966) in relation to semigroup theory and parallels the development of transition matrices in lumped systems. However, we shall not dwell here on this abstract approach but will concentrate on the problem of constructing Green's 29 function for the well-posed physical problem, assuming existence and uniqueness properties. Green's functions for distributed parameter dynamical systems can be constructed in a variety of ways. Often it is possible to use analytical methods, such as integral transfonns, to obtain the Green's function for linear systems. Expansion in spatial eigenfunctions and constructing the Green's function from eigenvalues and eigenfunctions is a very powerful method for linear systems (Morse and Feshbach, 1953; Stakgold, 1968). For these systems, one can find the eigenvalues and eigenfunctions if the spatial domain is a well defined geometrical region like a rectangle, a circle, or an ellipse, and the boundary conditions are specified in the corresponding principal coordinate directions. Twirand three-dimensional problems can also be treated and the eigenfunctions are generally obtained as separable products of functions of the coordinate variables. For systems in which it is not possible to find analytical expressions for the Green's functions, one can construct Green's functions by numerical methods. 2.5 Optimality Index The covariance matrix P(x,y,t) is _a function of the measurement locations xi' i = l,2,···,m. In order to define optimal locations, we must choose an optimality index that is a suitable measure of the covariance matrix P(x,y,t). For the case of a scalar state, a robust criterion is to minimize the spatial integral of the estimate error variance, i.e. J(t) = l P(x,x,t) dx (2.48) 30 For finite-dimensional systems, the analogue is to minimize the trace of the covariance matrix. For n-dimensional state vectors in distri- buted parameter systems, a natural extension is to l J(t) = Tr P(x,x,t) dx (2.49) The indices (2.48) and (2.49) arise naturally because using Mercer's theorem (Stakgold, 1968; Balakrishnan, 1976) they can be proved to be the sum of all the eigenvalues of the corresponding covariance operators. 2.6 Upper Bound for the Covariance Matrix Because of the nonlinearity and dimensions of the partial differential Riccati equation governing the covariance matrix P(x,y,t), it is difficult to obtain results for the optimal sensor location problem without goi~g into extensive numerical computations (Chen and Seinfeld, 1975). In order to develop a computationally tractabl~ approach, we prove a theorem relating to an upper bound on the covariance matrix of the optimal linear distributed filter. We now state and prove the basic theorem, based in part on the work of Sorenson (1968) in which a similar upper bound is derived for a discrete-time finite dimensional system. We prove the theorem for the case of spatially independent integral measurements from which the case of measurements at discrete points in space is readily tlerived, as shown in Section 2.3. It may be pointed out that the result derived here holds for all types of measurement models which can be described by the equation (2.6). Before stating the theorem, 31 let us define the inequality notation to be used. If r 1(x,y) and r 2(x,y) are two nxn nonnegative definite matrices, then (2.50) implies that for any n-dimensional vector f(x), ff fT{x)r 1{x,y)f{y) dydx ..::_ff DD fT(x)r 2{x,y)f{y) dydx (2.51) DD From the above definition, it follows that J Tr r 1{x,x} dx . ::_ f Tr r2{x,x) dx D (2.52) D Now we state the theorem for the upper bound for the covariance matrix P(x,y,t). A proof of this theorem is given in Appendix A. !~:~~:~= The covariance matrix for the optimal filter has an upper bound described by P{x,y,t):. P{x,y,t) + W(x,y,t) (2.53) where P(x,y,t) is the optimal filter covariance matrix for the system (2. 1) in the absence of dynamic noise, and W{x,y,t) is the matrix de- fined by t W(x,y,t} = J (J G(x,t;x' ,-r}B{x' ,-r}Q{x' ,y' ,-r}BT(y' ,-r)GT(y,t;y' ,-r} dy'dx'd-r to Proof: bo See Appendix A. {2.54) 32 The solution of the optimal filter covariance equation P(x,y,t) for the process in the absence of dynamic noise can be written in a closed form in terms of the system Green's function matrix. The equation for P(x,y,t) for the case of point measurements can be obtained by setting Q(x,y,t) = 0 for the covariance equation in case 4. Thus, we have aP(x,y,t) = A P(x y t) + P(x y t) AT at X ' ' ' ' Y m T -1 i-(x,n.,t)H.(t)(R (t)) .. H.(t)P(nJ.,y,t) I 1 l lJ J i ,j=l (2.55) F\ I.C. B.C. ax P{x,y,t) = 0 xE ao ye Dxao (2. 56) -P(x,y,t)ayT = 0 ye ao x€ DxaD (2.57) fr The solution '· of equations (39)-(41) can be written in the form P{x,y,t) = J J G(x,t;x ,t0 ) [P~(x ,y•) + M(x ,y• ,t)]tGT(y,t;y' ,t0 ) dy'dx' 1 1 1 DD (2.58) where G(x,t;x' ,t 1 ) is the system Green's function matrix, and M(x,y,t) The matrix M(x,y,t) defined above is closely related to the concept of observability. The system (2.1)-(2.6) is observable in the sense of 33 Wang (1964), if and only if M(x,y,t) has a bounded generalized inverse for some finite t 2:.. t 0 • However, for (2.58) to be valid we need not assume the existence of a generalized inverse for M(x,y,t). Indeed, for M(x,y,t) = 0, (2.58) is still a valid solution of (2.55)-(2.57) because in this case [P:(x',y') + M(x',y',t)Jt = P0 (x',y'). The validity of (2.58) as a solution of (2.55)-(2.57) is verified in Appendix B. 2.7 Approximation of Pin Finite-Dimensional Subspace The expression (2.58) for P(x,y,t) involves generalized inverses. Since the generalized inverse is defined by the integral equation (2.22), it is convenient to approximate it in a finite-dimensional subspace of the original infinite-dimensional space. Let us consider a real scalar valued symmetric function k(x,y) whose generalized inverse is denoted by kt(x,y). We can represent k(x,y) as a doubled series (Courant and '1' Hi 1bert, 1966; Stakgol d, 1968) 00 k(x,y) = I k . . cp. (x) cp .(y) i,j=l lJ l J where the basis functions {<Pi} fonn a complete orthononnal set. The expansion coefficients kij are given by kij = JJ k(x,y)cpi(x)cpj(y) dydx (2.60) DD The function k(x,y) can be approximated in the nonn by a finite double series (2.61), 34 We can also write equation (2.61) as k(x,y) = ¢T(x) K¢(y), where ¢(x) is the N-dimens ional column vector of elements <t>i (x) 1 ~ 1 ~ N, and K is the Nx N matrix with elements k ... The generalized inverse kt(x,y) lJ can then be approximated by (2.62) To verify this we may consider k(x,y) to be the kernel of the integral operator K, and kt (x,y) to be the kernel of the integral operator K-1 . Let us also consider an arbitrary function f(x) which can be approximated as (2.63) ' where a is the N-dimensional column vector of elements al.. Then by the definition of inverse operator K- 1(Kf) = f kt(x,y) f k(y.z)f(z) dzdy = f(x) D (2.64) D Substituting the expressions for kt(x,y), k(x,y), and f(x), we find that the definition (2.62) is, indeed, satisfiei. We may also note that if k{x,y) is time-dependent, the elements of the matrix K will be timedependent. For the case of an n x n dimensional matrix-valued function k(x,y), the approximation can again be written as k(x,y) = ¢T(x) K¢(y); however, in this case ¢(x) is an nN xn matrix given by tThe basis functions {<f>i} are not necessarily the eigenfunctions of the operator K 35 0 • • • • • 0 0 • • . . . • 0 <I>(x) = T 0 0 . . 0 • • • • • • 0 and K is nNxnN composite matrix made of NxN submatrices Kij l < i,j < n, where Kij is the NxN matrix of expansion coefficients for the (i,j)th element in k(x,y). We may note that the orthononnality of {¢i} gives f <I>(x) <I>T(x) dx = I (2.65) D for both the scalar and matrix cases. We now consider the case of scalar P(x,y,t) but the results can be interpreted for matrix case with the above definitions ip mind. In order to obtain an approximate solution to P(x,y,t) for a scalar case, we consider the complete orthononnal set ¢.(x) i=l,2,• .. ,oo. 1 In the N-dimensional subspace, the basis functions considered are ¢i(x), i= 1,2,···,N. Denoting the N-dimensional column vector of the basis functions ¢i(x~ by <I>(x), we can approximate the scalar Green's function G(x,t;x',t') as G(x,t;x' ,t') = <I>T(x)A(t,t')<I>(x') (2.66) where A(t,t') is the NxN symmetric matrix of time dependent coefficients. We shall also approximate the initial covariance P (x,y) as 0 36 Po(x,y) = ¢T(x) ~o ¢(y) so that Pt (x,y) = ¢T(x) ~o-1 ¢(y). 0 (2.67) The expression ( 2.59 ) for M( x,y,t ) can then be written as M(x,y,t) = t f t Jm T(x)A(T,t · )qi(ni)(R- l (T)) 1J.. ¢T(n.)A(T,t )qi(y) dT J . qi 1,J - 1 0 O o (2.68) or, M(x,y,t) = qiT(x) S(n.;t) qi{y) 1 (2.69) where S(ni;t) depends on the sensor locations ni' i = 1,2, ... , m and is given by Therefore, and, With the help of this expression, we can express P(x,y,t) given in (2. 58) as 1 P(x,y,t) = Jf qiT(x)A(t,t )qi(x')qiT(x') [~~l + S(ni;t)]- ¢(y')qiT(y') 0 DD •A(t,t 0 )¢(y) dy'dx' 37 Using the orthononna l i ty property, we get (2.73) Thus, we have obtained an expression for P(x,y,t) in terms of the N basis functions and the expansion coefficients of Green's function. The dependence on sensor locations ni' i= 1,2,·· •,mis through the dependence of matrix Son these sensor locations. 2.6, P(x,y,t) ~ P(x,y,t) + W{x,y,t). As proved in Section It then follows that l P(x,x,t) dx 2 l ~(x,x,t) dx + l W(x,x,t) dx or, J(t) 2 f P(x,x,t) dx + r(t) (2.74) D r(t) = where f W(x,x,t) dx D Of the two tenns in the upper bound, r(t) is independent of the sensor locations. Therefore to minimize this upper bound of J(t) we need only to minimize J P(x,x,t) dx with respect to sensor locations. From equa- tion (2.73) D J P(x,x,t) dx = J ¢T(x)A(t,t0 ) D D r~~l + S(ni;t~-1(t,t0 )¢{x) dx (2.75) The right hand side is a nonlinear function of the sensor locations ni, i = 1 ,2, • • • ,m. The optima 1 sensor 1ocati ons are the values ni that minimize this nonlinear function, and can be determined by using a 38 suitable gradient method. A slightly modified version of the Fletcher- Powell method (Fletcher and Powell, 1963) was used to obtain the optimal sensor locations for the example considered in this work. The optimality index (2.49) neglects the cross-correlation of estimate errors at different points in space. at different If the cross-correlation spatial points is important, a suitable optimality index is J 1(t) =ff Tr P(x,y,t) dydx (2 . 76) DD Then utilizing the upper bound ll Tr P(x,y,t) dydx_'.'. ll Tr P(x,y,t) dydx ll + Tr W(x,y,t) dydx ( 2. 77) and (2.73) fqr P(x,y,t), optimal sensor locations for minimum J 1(t) can be readily detennined. 2.8 Optimal Location of Sensors in a One-Dimensional Diffusion System We cons i der heat conduction through a long bar of uniform crosssection, along which heat is lost to the surroundings. The unsteady state temperature u(x,t) is described by (2.78) Assume that the initial temperature is not known exactly, and we have E[u(x,O)] = O 39 (2. 79) (Equation (2.78) also describes the process of diffusion of a species A in a liquid medium B where A reacts homogeneously with B. In some cases, neutron transport can also be described by the same equation.) Let us suppose that equation (2.78) does not model the process exactly, and we introduce a random dynamic error w(x,t), a Gaussian process with properties E[w(x,t)] = 0 (2.80) E[w(x,t)wT(y,T)] = Q(x,y,t) o(t-T) where Q(x,y,t) is a symmetric non-negative scalar-valued function. Assume that the temperature u(x,t) is measured continuously at a finite number of points, xi i=l,2,· 00 ,m, and the measurements zi(t) are described by "~ i = 1,2,···,m (2.81) where the measurement noise vi(t) is a Gaussian process with properties E[vi(t)] = 0 i 1 j E[v.(t)v.(T)] 1 J . ={o r i ( t) 0 ( t-T) (2. 82) i = j The variance of the measurement noise can be evaluated from a continuous function R(x) of the measurement location x, i.e., ri = R(x)lx=xi· It is desired to select optimal sensor locations xi, i = 1,2,· 00 ,m, 40 As outlined before, we minimize an upper bound of J(t) by minimizing J(t) given by J(t) = '{ P{x,x, t) dx (2.83) D We have studied the case of two sensor locations for the parameter values K = 0. 1, and 8 = 0.25. For studying the effect of different boundary conditions, measurement noise covariance functions, and initial covariances, sensor locations for minimizing J(t) at t = 1 were considered. 2.8. 1 Effect of Boundary Conditions Three different boundary conditions are considered. The boundary conditions and corresponding eigenfunctions and Green's functions are given in Table 2.2. tions {¢i(x)} . length. The eigenfunctions were employed as the basis func- Two sensor locations were determined for a bar of unit For ease of computation, the initial covariance P0 (x,y) was chosen to be 2o(x-y). Although this choice of P0 is not trace class, it can be viewed as the limiting case of the covariance matrix of a random process that is correlated over a small interval. Moreover, for compu- tations we must use a finite dimensional representation of P0 (x,y), which is always trace class. Therefore, the difficulties associated with P0 not being trace class do not arise. The measurement noise co- variance was taken to be unity independent of the sensor location, i.e., R(x) = 1. The results are sunmarized in Table 2.3. (C) (B) au~I = o ax x=l au(x,t)I =O ax x=o u(l,t) = 0 u(O,t) = 0 u(O,t) = 0 (A) au~ =0 ax x=l Boundary Conditions 2 n = 1,2, ... l ; ./l cos n1rx n = 1,2, .... ./l sin n1rx n = l ,2, ... n=l I 00 n=l n=l 2 2 1 ) I exp[-(n 1r K+B)(t-t')]cos(n1rx)cos(n1rx') 1 •sin [ ( n- l /2) 1rx' ] exp[-(n~1r 2 K+B)(t-t )]sin(n1rx)sin(n1rx 2 2 I exp[-((n-l/2) 1r K+B)(t-t')]sin[(n-l/2)1rx] 00 Green's Function G(x,t;x' ,t') e-B(t-t') + 2 2 ~ ./l sin(n-l/2)1rx Eigenfunctions Table 2.2 Eigenfunctions and Green's Functions for Different Boundary Conditions .... ~ 42 Table 2.3 Optimal Sensor Locations for Different Boundary Conditions Initial Covariance Measurement Noise Covariance Boundary Conditions (A} u(O,t) = 0, (B) u(O,t) = 0, Sensor Locations u(x,t) X R(x) = l = 0 x=l u(l,t) = 0 = au(x,t) (C) au(x,t) • ax x=l = 0 ax x=O 0.57 , 0.87 0.43 , 0.57 0.21 , 0.79 43 It is observed that the optimal sensor locations are, in general, removed from the boundaries when the boundaries are kept at fixed temperature, but are close to the boundary when the boundary condition is of zero flux. When the boundaries are held at specified temperature, we do not expect to gain much information about the system by making measurements near the boundary. On the other hand, if the boundary condition is specified to be zero flux, the temperature near the boundary is constantly adjusting because of the heat flow through the bar. Consequently, measurements near the boundary convey more information about the dynamics of the heat transfer process. 2.8.2 Effect of R(x) For the boundary conditions of case A, the effect of varying R(x) on the optimal sensor locations was studied. .'! In these experiments, the initial covariance P0 (x,y) was maintained at 2o(x-y). The results are given in Table 2.4. For constan~ R, the optimal locations are found to be 0.57 and 0.87. For cases (b) and (d), R(x) decreases as we move towards the end x = l; therefore the second sensor location moves to x = 1, where R(x) is minimum. In case (d), R(x) is relatively high at the point x =0.57, and therefore the optimal location for the first sensor changes to x= 0.37 where R(x) has a lower value. In cases (c) and (e), the end x=l has very high value of R(x) and therefore the second sensor moves to positions of lower R(x). The first sensor moves to the position x = 0.50 where R(x) achieves its minimum for case (e) and is quite low 44 Table 2.4 Effect of Measurement Noise Covariance R(x) on Optimal Sensor Locations Boundary Conditions: f:J(O,t) = 0 au(x,t)I =0 ax x=l Initial Covariance: R(x) Sensor Locations (A) 1.0 0.57 , 0.87 (B) 0.8 exp[-2x] o. ~o , 1.0 (C) 0. 8 exp[ -2 ( 1-x) J 0. 50 , 0.85 (D) 0.8 exp[-~lx-0.51] 0.37 , 1.0 (E) 0.8 exp[+4lx-0.5!] ,'! 0. 50 , 0. 71 45 z 0.80 0.07 0 0.06 ~ (.) 0.08 0.06 0 _J 0.60 a::: 0.20 0 (/) z w (/) 0.40 0 ,ff z 0.40 0 (.) w (/) 0.20 0.60 0.00 Q20 Q40 OEO Q80 FIRST SENSOR LOCATION FIGURE 2.1: CONTOUR PLOT OF J FOR R(x) = 0.8 EXP (-2x) 1.00 46 0.12 0.16 z 0.80 0 0.12 ~ 0 0.12 0 _J 0.60 er 0.165 0 (/) z 0.22 w (/) 0.40 0.29 0 z 0 0 w (/) 0.20 0.00 0.20 0.40 0.60 0.80 FIRST SENSOR LOCATION FIGURE 2.2: CONTOUR PLOT OF J FOR R(x) ,,;, 0.8 EXP[-2<1-x)] 1.00 47 z 0.80 0 0.30 ~ (_) 0 _J 0.60 0:: 0 Cf) z w Cf) 0.40 0.45 ,., 0 z 0 (_) 0.60 w Cf) 0.20 0.00 0.20 0.40 0.60 0.80 FIRST SENSOR LOCATION 1.00 FIGURE 2.3: CONTOUR PLOT OF J FOR R<x) = 0.8 EXP[+41x-0.5J] 48 for case (c). Although it is di f ficult to separate out the effects of the dynamics of the process, and the measurement noise covariance, the sensor locations, in general, move towards points of lower measurement error variance. • Contours of constant J(t) are shown in Figures 2.1, 2.2, and 2.3 for cases (b), (c), and (e). The plots are symmetric about the 45° line as the two sensors are equivalent. The optimal locations are in - dicated by crosses. 2.8.3 Effect of~ (x,y) To study the effect of in i tial covariance P0 (x,y), the boundary conditions of case A were considered, and R(x) was specified as R(x) = 0.8 exp[-2(1-x)]. Different values of P0 (x,y) were assumed and sensor locations were obtained. The results are given in Table 2.5. 2.8.4 Comparison of the Results to Solution of the Riccati Equation Now we consider tDe full equation for the covariance matrix P(x,y,t) and study the effect on the sensor locations of considering the upper bound of ~(x,y,t) rather than P(x,y,t) itself. The full nonlinear Riccati equation governing the covariance matrix P(x,y,t) was solved using an alternating-direction implicit method. The values for the optimality index J(t) = J P(x,x,t) dx (2.84) D were evaluated at t=0.10 for varying Q(x,y,t). Initial covariance 49 Table 2.5 Effect of Initial Covariance P0 (x,y) on Opti mal Sensor Locations Boundary Conditions: u(O,t) = O au(x,t)I =0 ax x=l Measurement Noise Covariance: R(x) = 0.8 exp[-2(1-x)J P0 (x,y) Sensor Locati ans (A) 0.2 o(x-y) 0.60 , 0.75 (B) 2o (x-y) 0. 50 , 0.85 (C) 200 (x-y) 0.48, 0.84 ,'! 50 Table 2.6 Values of J(t) Obtained by Solving the Full Riccati Equation for Various Combinations ·of Sensor Locations Q{x,y,t) = 0 x1/x 2 0. 1 o. 1 0.3 0.5 0.7 0.9 2.817 2.698 2.758 2.714 2.460 2.364 2.310* 2.544 2.370 0.3 2.817 0.5 2.698 2.460 0.7 2.758 2.364 2.544 0.9 2.714 2.310* 2.370 2.558 2.558 Q{x,y,t) = 0.6 exp [-lO(x-y)] x1tx 2 0.1 0. 1 0.3 0.5 0.7 0.9 2.839 2.719 2.780 2.736 2.481 2.384 2.330* 2.566 2.391 0.3 2.819 0. 5 2.719 2.481 0.7 2.780 2.384 2.566 0.9 2.736 2.330* 2. 391 2.580 0.3 0.5 0.7 0.9 2.986 2.864 2.926 2.881 2.622 2.524 2.469* 2.709 2.530 2.580 Q(x,y,t) = 0.5 o(x-y) x1/x 2 0. l 0. 1 0.3 2.986 0.5 2.864 2.622 0.7 2.926 2.524 2.709 0.9 2.881 2.469* 2.530 2. 723 2. 723 51 P0 (x,y) was maintained at 2o(x-y), and R(x) = 0.8 exp[-2(1-x)] was used. The sensor locations x1 and x2 were restricted to five points 0. l, 0.3, 0.5, 0.7, and 0.9. The values of J(t) for different combinations of sensor locations, and different Q(x,y,t) are given in Table 2.6. For the three cases considered, J(t) for the optimal sensor locations is indicated by* . In all three cases, optimal sensor locations , are seen to be 0.3 and 0.9 based on the prespecified locations. The optimal sensor locations determined by using the method developed in this chapter are 0.296 and 0.866. Thus, we see that the method works quite well for reasonable covariance functions Q(x,y,t) and provides considerable savings in computation time. 2. 9 Conclusions In this chapter we propose a computationally attractive approach for the optima~ sensor location problem in distributed parameter systems. An upper bound to the filter covariance matrix for linear distributed parameter sys terns is deve 1oped. When the optimal sensor location prob- lem is posed as a minimization problem for a suitable measure of this upper bound, the computational effort involved is significantly smaller compared with the effort necessary for minimizing a measure of the filter covariance matrix (Aidarous, 1976; Chen and Seinfeld, 1975). An example involving one-dimensional heat conduction is considered and the results obtained by the two approaches are compared. The upper bound in (2.53) may indeed be a weak one for certain cases. However, the term P(x,y,t) in the upper bound takes into account the effect of both the 52 system dynamics and the measurement noise covariance on the sensor locations. The numerical results of the example show that for spati- ally independent Q(x,y~t), the optimal sensor locations obtained by considering the full matrix Riccati equation are not ~ifferent from those obtained by using the upper bound approach developed here. CHAPTER 3 OPTIMAL LOCATION OF MEASUREMENTS IN TUBULAR REACTORS 54 3.1 Introduction In this Chapter. we consider the problem of determining optimal sen- sor locations for a meaningful chemical engineering system. Most systems of interest are described by more than one state variable. In practice, it may be difficult to obtain continuous measurements of some of the state variables while measurements of other variables may be disturbed by noise. Often. for control purposes, it is desired that the state vari- ables be estimated accurately from a knowledge of the system model and some measurements on the system. In a situation where we wish to esti- mate the system state using only a limited number of measurements. an important question is the selection of sensor locations. It may also be desirable to determine the number of sensors required in order to estimate the state to a given accuracy. The particular system considered in this Chapter is a non-adiabatic tubular reactor in which a first order exothermic chemical reaction is taking place. The reaction term introduces nonlinearity in the mass and energy balance equations and makes the system very interesting. The non- adiabatic reactor model was chosen because it is more realistic than the isothermal or adiabatic models. The steady state profiles in non- adiabatic reactors are often characterized by hot spots and the random input disturbances may cause considerable deviations away from the steady state concentration and temperature profiles. It is. therefore. desirable to obtain good estimates of the concentration and temperature profiles. In this work, we deal with the problem of determining optimal sensor locations with the objective of estimating concentration and temperature profiles as accurately as possible. 55 Keeping in mind the random nature of input disturbances and inherent uncertainties in the system parameters and the measurements, the sensor location problem is approached from a stochastic point of view. We use the framework of estimation theory and consider the location of sensors so that a best estimate of the concentration and temperature profiles is obtained. Since the tubular reactor is a nonlinear system with two state variables, it is difficult to directly apply the theories available (Aidarous, 1976; Chen and Seinfeld, 1975). The approach developed in Chapter 2 will require the numerical computation of Green's function matrix for the linearized problem which can be very time consuming. An alternate approach based on discretization of the reactor equations is developed and the optimal sensor locations are selected from a specified set of locations. In Section 3.2, we introduce the equations describing the tubular reactor and nondimensionalize them by introducing appropriate dimensionless variables. The steady state equations are obtained and the reactor equa- tions are linearized around the steady state. In Section 3.3, we give a brief outline of the method of orthogonal collocation, and following Georgakis et al. (1975), the reactor equations are discretized using the orthogonal collocation method. The nonlinear discrete equations are linearized around the steady state. Equations for the optimal state esti- mate and the corresponding covariance matrix are obtained in Section 3.4 using the familiar Kalman filter (Jazwinski, 1970). The optimality index for sensor locations is defined in Section 3.5 and an iterative algorithm for determining optimal sensor locations is developed based on Pontryagin's minimum principle. In Section 3.6 we describe in detail the results of a numerical example. The parameters in the model are chosen so that a unique 56 stable steady state is guaranteed (Varma and Amundson, 1973b). The effect of process dynamics, process noise, and the measurement noise on the optimal sensor locations is studied. 3.2 Reactor Model The system under consideration is a non-adiabatic tubular reactor in which a single first order chemical reaction is taking place. We have chosen to deal with a non-adiabatic reactor model as against the isothermal or adiabatic models because it is more representative of the practical situation. Varma and Amundson (1972, 1973a, 1973b) have considered the ques- tion of uniqueness and stability of the steady state of non-adiabatic tubular reactors. Georgakis et al. (1975, 1976) have studied the problem of stabilization of unstable tubular chemical reactors by modal control and observer design. Fluid consisting of species A1 enters the tubular reactor of radius R and length Lat concentration c~1n and temperature T!1n , and undergoes a first order chemical reaction with reation rate r'(c' ,T') and heat of reaction 6H. The concentration of species A1 is denoted by c', and T' denotes the temperature of the reaction mixture at any point along the reactor length. The reactor is cooled by a heat exchanger in which the coolant is maintained at an essentially constant temperature T~. constant. The over-all heat transfer coefficient U is assumed The density, heat capacity, and thermal conductivity of the reacting mixture are denoted by p, Cp' and k respectively. In order to 57 write down the transient heat and mass balance equations, we assume that the fluid density, heat capacity, thermal conductivity and heat of reaction all remain constant. The axial diffusivity is denoted by D and radial gradientsareassumed to be negligible. Under these assumptions, the transient equations governing the mass balance for species A1 and the overall energy balance are ac• - D a2 c' arr "'T' 'l. 2 T' ac' a"s7"7 - V asr - r'(c',T') 'l.T' 1 1 pCp;t' = kt-rr- pVCp;s' + {-tiH)r {c ,T (3 ·1) 2U 1 ) +R(T~ -T') (3.2) where t' is the time, s' is the distance along the reactor length, and Vis the velocity of the reaction mixture. The boundary conditions for (3. 1) and (3.2) are at s' = O Ve'.in = (vc' - D ac') asr s'=O+ t' > 0 (3.3) pVC T~ =(pVCPT' - k aT:) p 1n as s'=O+ t' > 0 (3.4) at s' = L ac'= 0 as' t' > 0 (3.5) aT'= 0 as' t' > 0 (3.6) Equations (3. 1)-(3.6) can be non-dimensionalized by introducing T' T~,n c' T =- s' t'D t = .--r- C =c'. in s = L (3. 7) L 58 The dimensionless reactor equations are (3.8) aT = "7T a2 T - P -aT + aSr + y(T -T) Le -at as eh as a (3.9) where Pe and P are the heat and mass transfer Peclet numbers, Le is h em the Lewis number, and S, y and rare the dimensionless heat of reaction, heat transfer coefficient, and reaction rate respectively. The dimension- less parameters are defined as follows pVC L p k pC D p eh Le = _:__p_ = k Pe L2 r'( c ! D c ,n '. m (-tiH) D c '. ,n S = --=-k--=T. --!- - T '. ) 1n' 1n a = ------...--- (3.10) 2UL 2 y=~ ,n r = r'~c' ,T'~ r 1 C ,n ! , T.,n ) For a first order reaction, the dimensionless reaction rate r(c,T) takes the form r(c,T) = c exp[o(l-1/T)] (3.11) where o( = E'/RgTin)' the dimensionless activation energy, is defined in terms of the activation energy E', the gas constant Rg' and the inlet temperature Tin of the reactant stream. The boundary conditions for (3.8) 59 and (3.9) are at s = 0, (3.12) aT _ as - peh ( T - l) (3.13) ~ =0 as (3.14) clT = O (3.15) at s = 1, as The steady state concentration c and temperature Tare governed by (3. 16) d2T P p ddT + aBr + y(T - T) = 0 - eh s a (3.17) with boundary conditions at s = 0, (3.18) dT = p (T _ 1) ds eh at s = l, de= 0 ds dT = 0 ds (3.19) 60 The dimensionless reactor equations (3.8) and (3.9) are now modified to account for the random inputs to the system. We assume that these random disturbances can be modeled as additive white Gaussian processes wc(s,t) and wT(s,t) with properties E[w(s,t)J = O (3.20) T + E[w(s,t)w (s' ,T)] = Q (s,s' ,t)o(t-T) where E[•] is the expectation operator and w(s,t) is the two dimensional column vector [wc(s,t) wT(s,t)]T. The white noise assumption implies that the disturbances at two different times t and Tare uncorrelated. The 2x2 covariance matrix Q+(s,s' ,t) is syrTmetric ins ands', and is nonnegative definite. The modified equations describing the reactor dynamics are (3.21) (3.22) with boundary conditions (3. 12)-(3.15). Let ~(s,t) and e(s,t) define the random deviations in concentration and temperature away from the steady state profiles c(s) and T(s) ~ ~(s,t) = c(s,t) - c(s) (3.23) e(s,t) = T(s,t) - T(s) where c(s,t) and T(s,t) are the concentration and temperature values in equations (3.21) and (3.22). The equations for ~(s,t) and e(s,t) can be 61 linearized around the steady state for small deviations. These linear- ized equations are (3.24) _J Le~~=~:~ - Pe :~ + as[w ( ~~)- _ + e / :~)+ wT(s,t) h c,T \ c,T (3.25) with boundary conditions at s = 0 ~= p tjJ as em (3.26) ae p e as - eh at s = l ( 3. 27) ae _ as -o The initial values for t/J and e are assumed to be white Gaussian processes with properties E [t/1 (s ' 0 )] = 0 e(s,O) (3.28) 62 E ~ l/J(s,O)Jr, ~(s' ,0) e{s' ,o~ = 8(s,O) (3.29) The initial covariance matrix P+(s,s') is a 2x2 synvnetric positive 0 definite matrix. 3.3 Discretization by Orthogonal Collocation The non-adiabatic tubular reactor has been described by a set of nonlinear partial differential equations. discretized to enable their solution. These equations need to be In conventional finite difference methods, the partial differential equations need to be discretized on a relatively fine mesh. The high dimensionality of the resulting lumped parameter system makes the computations time consuming. Thus, a discreti- zation method which results in a lower order system is highly desirable. The method of orthogonal collocation is one such method and can be advantageously used to reduce the dimensionality of the problem. We shall now give a brief outline of the method of orthogonal collocation for solving the transient concentration and temperature profiles governed by the deterministic equations (3.8) and (3.9). The supporting theory and mathematical justification are treated in detail by Villadsen (1970) and Finlayson (1972). The orthogonal collocation method is a special case of the method of weighted residuals. The unknown exact solution is approximated by a trial solution consisting of a series of known functions of the spatial variable(s). These known functions are chosen to be a family of suitable 63 orthogonal polynomials. The trial solution is substituted into the par- tial differential equation and the unknown time-dependent coefficients are determined so that the partial differential equation is satisfied at a set of grid points called the collocation points. The number of terms in the trial solution can be increased arbitrarily, and the unknown timedependent coefficients are determined by requiring that the partial differential equation be satisfied at more and more collocation points. The orthogonal collocation method utilizes a continuous representation of the approximate solution over the whole region of space. Therefore, the spatial derivative at any point is expressed in terms of the values at all other collocation points, rather than in terms of just the neighboring points as is the case in finite difference methods. Thus, the ortho- gonal collocation method requires fewer equations, although the advantages of dealing with a tridiagonal matrix are lost. However, as the number of collocation points is increased, the rate of convergence increases very fast for the orthogonal collocation method. The rate of convergence for the usual second order finite difference methods is l/(N-1) 2 where N is the number of grid points. On the other hand, for orthogonal collocation method the rate of convergence is (1/N)N (Ciarlet et al., 1967), where N is the number of interior collocation points. The choice of trial func- tion depends on the particular problem under consideration. For problems in which the solution is symmetric about a central axis, the solution can be expanded in terms of powers of s 2. However, for problems whose solutions do not possess synmetry properties, we must include both odd and even powers of s. Since for the tubular reactor, we do not expect the pro- files to be symmetric with respect to s, we choose a trial function of 64 the form c(s,t) = d(t) + e(t)s + s(l-s) N I g.{t) Y~p,q)(s) i=l l 1-l (3.30) where Y~p,,q)(s) are the first N orthogonal Jacobi polynomials over the ,_ interval (0,1) with weight function sq(l-s)P. The trial function (3.30) includes both odd and even powers of s, and has N+2 constants. The re- quirement that the solution satisfy the partial differential equation (3.8) at N interior collocation points provides N conditions, and the other two are provided by the boundary conditions at s = 0 ands= 1. The !!!_th polynomial Y~p,q)(s) in the family of orthogonal Jacobi polynomials is a linear combination of powers of s with highest power m m I b. sj j=o J starting with y{p,q)(s) = 1, the coefficients in the successive polyo nomials are determined by requiring Y~p,q)(s) to be orthogonal to all the polynomials of order less than m. The Jacobi polynomials y(p,q)(s) m are determined by solving the equation (3.31) 0 n = 0, 1,2, ... , m-1 The polynomials y(p,q)(s) are determined to a multiplicative constant, m which is determined by requiring the first coefficient to be l. For 65 p=q=O, these polynomials are the shifted Legendre polynomials. In the trial functions, we have chosen p=q=O and the polynomials Y.(s) are l written without any superscript. nomial of degree N+l ins. The trial function (3.30) is a poly- We rewrite (3.30) as a polynomial in s with time dependent coefficients f.(t) J N+2 I c(s,t) = fJ.(t)sj-l (3.32) j=l At the N collocation points si; i = 2,3, ... , N+l and the two boundary points s 1 = 0 and sN+ 2 = 1, we can write the concentration c(s,t) and its first two spatial derivatives as N+2 c(s.,t) = l ac as I j=l N+2 = (j-l)f.(t)s~- 2 J l I j=2 Si f .(t)s~-l J l (3.33) N+2 a2 c 3 I (j-l)(j-2)fJ.(t)s~asT s. = j=3 l l In terms of the matrix elements Gij' Kij' and Lij' (3.33) can be written as N+2 c(si,t) = l j=l G.. f.(t) lJ J G.. = sj-l lJ s.l 66 N+2 ac as s. l = I j=l dsj- l ds K.. f .(t) K.. = ---,--, L. .f.(t) L.. lJ J lJ N+2 = I j=l lJ J lJ (3.34) s. l One can write a similar trial solution for T(s,t) and substitute the two trial solutions into (3.8) and (3.9). The coefficients d(t), e(t), and gi(t) are then determined by satisfying the partial differential equations (3.8} and (3.9) at the N interior collocation points, chosen to be the zeros of the Nth order polynomial YN(s). However, for computa- tional purposes it is more convenient to solve for the concentrations and temperatures at the collocation points. The first and second derivatives -in (3.34) can be expressed in terms of ci(t) = c(si,t) as follows, N+2 ac = l A.. c.(t), ass. j=l lJ J i = 1,2, ... , N+2 l (3.35) N+2 l j=l B .. CJ. ( t) i = 1,2, ... , N+2 lJ The (N+2) x (N+2) dimensional matrices A, B with elements Aij' Bij' are given by (3.36) 67 where K, L, and Gare the (N+2)x(N+2) matrices with elements K.. , L.. , lJ lJ and G. . in ( 3. 34) lJ The boundary conditions at s=O and s=l can also be approximated and their linearity allows the concentrations and temperatures at the two ends to be eliminated. tions. Georgakis et al. (1975) give the detailed equa- Because of different non-dimensionalization, our equations are slightly different and, for the sake of completeness, are given in Appendix C. The discrete equations, for the case Le= 1, are de. dt __J_ = N+l I k=2 ( 3. 37) D. kck - ar( c., T.) + P· J J J J (j = 2,3, ... , N+l) dT. dt __J_ = N+l I k=2 S.kTk + a6r(c.,T.) + y(T -T.) + q. J J J a J J (3.38) (j=2,3, ... , N+l) where Djk' Sjk' pj, and qj are defined in the Appendix C and are functions of the elements of matrices A and Band P and P . The solution of eh em (3.37)-(3.38) with proper initial conditions gives us the concentration and temperature values at the interior collocation points. values are determined by equations (C.4) and (C.5). The boundary Once the con- centration values at collocation points are known, one can evaluate the coefficients f.'s in the polynomials approximating the concentration proJ file. The concentration at any other point can be determining readily by using (3.32). a similar way. The temperature at any point can also be determined in 68 The steady state profiles can be detennined by solving (3.37)-(3.38) with arbitrary initial conditions over a sufficiently long period of time so that the left hand sides reach an arbitrarily small constant. The steady state values cj and fj can also be determined by solving the algebraic equations, N+l l D.kck - ar(c.,i.) + p. = 0 k=2 J N+l L k=2 J J J s.kTk + aSr(c.,T.) + y{T -i.) + q . = 0 J J J a J J (3.39) (3.40) The solution of (3.39) and (3.40) gives the steady state concentrations c.J and temperatures T.J at the N interior collocations points s.,J j = 2, ... , N+l. The steady state concentration at any other point can then be determined by evaluating the coefficients fj's in (3.33) and using (3.32). The steady state temperature profile can be determined in a similar way. Let us now consider the random deviations w(s,t) and e(s,t) of (3.23). At the collocation points sj, these deviations are denoted by xj(t) and yj(t) x.(t) = w(s.,t) J J j = 2, ... , N+l y.{t) = e(s.,t) J J The dynamical equations governing xj and yj, j = 2, ... , N+l can be written as (3.41) 69 Equations (3.42) and (3.43) can also be obtained by linearizing (3.37) and (3.38) around the steady state values cj and Tj' and introducing the random disturbances wc(sj,t) and wT(sj,t). The initial values for xj and y. are white Gaussian processes with properties. J (3.44) (3.45) where P:(sj,sk) is the covariance matrix P:(s,s') evaluated at s=sj and s'=sk. 3.4 State Space Formulation and Optimal Estimate In conventional state space notation, the state vector u is the 2Ndimensional column vector, [x 2 y2 x3 y3 ..... xN+l YN+l]T. The transient equations (3.42)-(3.43) can be written as du _ dt - F u(t) + F,:(t) (3.46) 70 (3. 47) A A where the 2N x 2N matrix Fis a matrix the (i,j)th element of which is A A the 2x2 matrix Fij' i,j = 1,2, ... , N given by :] D.. lJ lJ ( -ar) ac c., T. ar ) , , D ii ij a a8 ( ac" c.j'. 1 1 - a ( -ar) aT (3.48) c.1 'i.1 A - 5; i + a8 ( ~~l c ., T. 1 A A i ; j .. [ 0 F"'" = A A i =j 1 A where i = i+l and j = j+l. The process noise ~(t) is the 2N-dimensional vector {3.49) The properties of w(s,t) assumed in Section (3.2) imply that ~(t) is a zero-mean white Gaussian process with 2N x 2N covariance matrix Q(t). The (i,j)th element of Q(t), i, j = 1,2, ... , N is the 2x2 matrix + Q (si+l'sj+l't). Thus we have E[~(t)J =0 (3.50) where 71 Q(t) = (3. 51) + • Q (sN+l'sN+l't) The initial state u0 is also a zero-mean white Gaussian process with 2N x 2N covariance matrix P0 (t). The (i,j)th element of P0 (t), i,j = 1, 2, .•. , N is the 2x2 matrix P:(si+l'sj+l't). We assume that point measurements of concentration and temperature can be made in principle at any given point along the reactor length. Let us consider the case when both concentration and temperature measurements are made at the N interior collocation points. The measurement vector is described by z.(t) = H. u.(t) + v.(t) l l l l i = 1,2,3, ... , N (3.52} where zi(t) is the vector of concentration and temperature measurements at si+l' and ui(t) is the two dimensional column vector with elements xi+l and Yi+l· The measurement matrix Hi can be chosen to specify the type of measurements - concentration or temperature or both. For meas- uring both concentration and temperature at the collocation points, Hi is a 2x2 matrix and is assumed to be time invariant without loss of generality. The measurement noise vector vi(t) is a zero-mean white Gaussian process with covariance 72 ' i = j E[v.(t)v.(T)] = l J (3.53) ' i ' j The above assumption implies that the measurement noise at two different locations si and sj is not correlated. The 2N-dimensional total meas- urement vector z(t) can be written as z(t) = Hu(t) + v(t) (3.54) where His a 2N x 2N matrix and v(t) is a 2N dimensional vector. In de- tailed form, (3.54) is v1(t) v2(t) = (3.55) + . vN{t) where Ei is a 2x2N matrix composed of N 2x2 submatrices, of which the ith submatrix is the identity matrix and all the rest are zero. The measure- ment noise vector v(t) is a zero-mean white Gaussian process with positive definite covariance matrix R(t). The 2N x 2N matrix R(t) is a quasi- diagonal matrix with the .:!_th diagonal block Rii(t), i = 1,2, ... , N. We now wish to obtain the best estimate of the state u(t), given the measurements z(T), 0 ~ T ~ t. Let us denote an estimate of the system state by ue(t). mate error covariance matrix is defined by The esti- 73 This error covariance matrix is nonnegative definite. tem and measurement scheme, the For a given sys- estimate ue(t) for which the covariance matrix is minimum is called the optimal estimate. Let us denote this A optimal estimate by u(t) and corresponding covariance matrix by P(t). Thus T P(t) = E[{u(t) - u(t)}{u(t) - u(t)}] A A (3.56) A It has been shown by Kalman that the optimal estimate u(t) for the system (3.46)-(3.54) is governed by the differential equation (Jazwinski, 1970} A d~~t) = F~(t) + P(t)HTR- 1(t}[z(t) - H~(t}] (3.57) The covariance matrix P(t) is governed by the nonlinear Riccati equation (3.58) P(O) = P 0 A The equations governing u(t) and P(t) are also known as the Kalman filter or the optimal filter equations. Using the notation of (3.55), the equa- tions defining the optimal filter can be written as A T T 1 du(t) = Fu(t) + N I P(t)E .H.R:.(t)[z.(t) - H.E.u(t)] dt i=l l l 11 l l l A A u(O) = 0 A (3.59) 74 and, dP(t) = FP{t) + P(t)FT + Q(t) - NL Pt ( ) E.H.R TT -1( .. t)H.E.P(t) dt 1111 11 i=l (3.60) P(O) = P 0 It should be noted here that the equation for P(t) is independent of the measurements z(t) but depends on the system matrix F, the measurement matrix Hi' and the noise covariances Q(t) and R(t). Consequently, we can evaluate the error covariance P(t) off-line. Let us assume that the concentration and temperature measurements can be made only at the interior collocation points and the optimal sensor locations are to be chosen from this set of N points. This may seem to be an unnecessary restriction, but the number of interior collocation points can be increased arbitrarily. Let us introduce the param- eters a.(t), i = 1,2, ... , N such that l , if a measurement is made at si+l at time t (3.61) , otherwise If the number of sensors is specified to be m, only m of the a.(t) values l will be 1. Thus the equations for u(t) and P(t) for measurements at only m locations can be written as N A ~t = F~(t) + l\ a.(t)P(t)E'.H'.R:~(t)[z.(t) - H.E.~(t)] UL l l l 11 l l l i=l A u(O) = 0 (3.62) 75 and N dP(t) = FP(t) + P(t)FT + Q(t) - I a.(t)P(t)E'.H'.R:~(t)H.E.P(t) dt . _1 1 1 1 11 1 1 ,(3.63) P(O) = P 0 3.5 Optimality Index and The Minimum Principle The measurement locations are specified through the parameters ai(t), i = 1,2, ... , N. Different measurement policies will lead to dif- "' ferent estimates u(t) and covariance matrices P(t). Since the covariance matrix P(t) indicates the accuracy of the estimate, it is natural that an optimal measurement policy be defined as one that minimizes a specified scalar measure of P(t). optimality index. This scalar measure will be called the Since P(t) is a nonnegative definite matrix, its trace is a good measure of its magnitude, and we can choose the optimality index J to be tf J = a 1Tr P(tf) + s1 J Tr P(t) dt (3.64) 0 where tf is the specified final time and a 1 and s1 are two positive constants specifying the relative wei,ghts of the two terms. The optimal sensor location problem can now be viewed as a deterministic control problem in which the state matrix P(t) is governed by the matrix Riccati equation (3.62) and the optimality index J is to be minimized with respect , to the parameters a.(t); i = 1,2, ... , N regarded as control variables. matrix minimum principle corresponding to Pontryagin's minimum principle A 76 has been formulated by Athans (1968) for optimal control problems in which the state is governed by a matrix differential equation. We shall now develop the necessary conditions for minimizing J by employing this matrix minimum principle. The algorithm used for determining the opti- mal sensor locations is based on the work of Chen and Seinfeld (1975). Let ~(t) denote the 2N x 2N adjoint variable matrix associated with the state matrix P(t). Then the Hamiltonian H(t) is defined as H(t) = s1TrP(t) + Tr[P(t) ~T(t)J (3.65) or H(t) • B1TrP(t) + Tr[(FP(t) + P(t)FT + Q(t) - P( t) ( I a . (t) H R: \ t) H. i=l E"!" "!" 1 1 1 11 E: ) 1 1 ~ P( t )) T( t )] The adjoint matrix ~(t) is governed by the matrix differential equation ¢(t) = - 1fiiil (3.66) aPITT +[ J ( 3. 67) TT 1 (t)H.E. P(t)~(t) I a.(t)E.H.R:. 1 1 1 11 1 1 N i=l with the terminal condition (3.68) 77 , The optimal values for the a. are then determined by the application of minimum principle which states that H(P*(t),~*(t),a~(t),t) , -< H(P*(t),~*(t), a.(t),t) 1 (3.69) where ai(t) is the set of optimal parameters ai(t), P*(t) and ~*(t) are the corresponding state and adjoint variable matrices. By using the definition of the Hamiltonian, the inequality (3.69) simplifies to N I i=l (3.70) This inequality forms the basis of the numerical algorithm used here for determination of optimal set of parameters a~(t). 1 It should be noted here that in the inequalities (3.69) and (3.70) the parameters a.(t) 1 have been assumed to be time dependent and therefore will lead to the time dependent optimal sensor locations. In most situations however, one is interested in obtaining time-independent sensor locations. For fixed sensor locations the parameters ai are independent of time and the inequality corresponding to (3.70) is N > I - i=l •; f \{P*(t)E;H;Rj)(t)H;E;P*(t)W*(t)] dt 0 (3.71) 78 Let us define the switching functions a.; i = 1,2, ... , N as l (3.72) In order to minimize J, we must choose them a.'s corresponding tom l largest cri's as equal to l and rest zero. Therefore, a simple iterative algorithm for determining optimal sensor locations is (l) Choose m ai's corresponding to the initial guess of sensor locations as l and rest zero. (2) Evaluate the covariance matrix P(t), the adjoint variable matrix ~(t), and the switching functions cri; i = 1,2, ... , N (3) Choose them ai's corresponding tom largest cri's as l and rest zero (4) If there is no change in the values of a l• IS, stop; otherwise go to step 2. The set of ai's determined by this iterative procedure gives us the optimal sensor locations. 3.6 Numerical Example and Results In this Section, we present the results for a numerical example. In Section 3.4, we have formulated the problem for the case where both concentration and temperature measurements are made simultaneously. How- ever, in practice it is costly and difficult to make continuous concentration measurements. Therefore, we consider the case in which only point measurements of temperatures are made, and it is desired to estimate both concentration and temperature profiles from these measurements. 79 3.6.l Steady State Solution The determination of optimal sensor locations requires the knowl- edge of the steady state values of concentration and temperature at the collocation points. Varma and Amundson (1972, 1973a, 1973b) have ex- plored in detail the steady state multiplicity and stability aspects of nonadiabatic tubular reactors. In this example, the parameters for the reactor model were chosen to guarantee a unique stable steady state. The parameter values are p =p em eh a= 0.875 =5 S = 0.5 y = 12.5 o = 25 Ta= l The number of interior collocation points N was chosen to be 7. collocation points are the zeros of the polynomial v7(s). The Table 3. l gives the polynomials Yi(s), i = 0,1, ... , 7 and the interior collocation points si, i = 2, ... , 8. The steady state values of concentration and temperature at the collocation points were determined by solving the 14 nonlinear algebraic equations (3.39)-(3.40). These algebraic equa- tions were solved by using a combination of Powell's method (Powell, 1970) and the Newton-Raphson method. Powell's method worked well for initial guesses far away from the true values but did not converge to the true solution. well. At this point, the Newton-Raphson method worked very The concentration and temperature values at other points were determined by interpolation. The steady state concentration and tempera- ture profiles are plotted in Figures 3.1 and 3.2. The steady state 80 Table 3. 1. Y (s) = l v1 (s) = v2 (s) v3(s) = 0 2s - 1 6s 2 - 6s + l 20s 3 - 30s 2 + 12s - 1 70s 4 - 140s 3 + 90s 2 - 20s + 1 252s 5 - 630s 4 + 560s 3 - 210s 2 + 30s = = Y4 (s) Orthogonal Polynomials and the Interior Collocation Points v5 (s) = - l v6(s) = 924s 6 - 2772s 5 + 3150s 4 - 1680s 3 + 420s 2 - 42s + 1 v7(s) = 3432s 7 - 12012s 6 + 16632s 5 - 11550s 4 + 4200s 3 - 756s 2 + 56s - 1 The interior collocation points s 2, s 3 , ... , s 8 are the zeros of v7 (s) given below: s2 = 0.0254 S3 = 0.1292 S4 = 0. 2971 S5 = 0.5000 s6 = 0.7029 S7 = 0.8708 s8 = 0.9746 u C) :z:: u w :z: 0::: I- c:::c I- C) ....... :z:: • ~ C) C) C) ('\J C) . C) (T") C) :::r . C) C) . C) Lf) o.o I \ 0.40 s 0.60 0.80 1 .00 SOLITTION OF ALGEBRAIC EQUATIONS X SOLUTION OF DIFFERENTIAL EQUATIONS FIGURE 3,1: STEADY STATE CONCENTRATION AS A FUNCTION OF S 0.20 - OJ ...... w I- w :E: 0.... w 0:::: c::r: :=:> I- 0:::: . . b 0 • . - o LI) • -- a - LJ') ..... (\I - a •. (\I - LJ') (TJ - a 0.20 0.40 s 0.60 ~ 0.80 FIGURE 3,2: STEADY STATE TEMPERATURE AS A FUNCTION OF S 0.0 I 1.00 - - SOLUTION OF ALGEBRAIC EQUATIONS X SOLUTION OF DIFFERENTIAL EQUATIONS CX) N 83 concentration decreases monotonically along the reactor length while the temperature profile shows the exi stence of a hot spot near the entrance. The steady state temperature is greater than the inlet temperature at all points along the reactor length. As mentioned in Section 3.3, the steady state profiles can also be determined by solving the 2N differential equations (3.37)-(3.38) over a sufficiently long period of time. These equations were solved with proper initial conditions, and the resulting concentration and temperature values at a finite number of points are shown in Figures 3. land 3.2. The results obtained by the two methods are in good agreement. For the initial conditions cj(O) = 0.20, Tj(O) = 1.10 j = 2,3, ... , N+l, the steady state was reached approximately at t = 1.4. 3.6.2 Determination of Optimal Locations The initial covariance P~(s,s') and the process noise covariance Q+(s,s' ,t) are assumed to be diagonal matrix functions of the form (3.73) (3.74) where k and k are constants and ii' i = 1, ... , 4 are the characteristic 2 1 lengths for correlation between differential spatial points. The matrix Q+(s,s' ,t) is assumed to be time invariant . Let us further assume, for simplicity, that the largest characteristic length i is much smaller 84 compared to the distance between any two collocation points j, k = 2,3, ... , N+l j 'f k In this case, the 2N x 2N covariance matrices P0 and Qare approximated by P o = kl I (3.75} Q= k I 2 This choice of P0 implies that any correlation between the uncertainties in the initial state at different collocation points is negligible. Similarly the choice of Q assumes that the correlation between system model uncertainties at various collocation points is negligible. The time invariant meas- urement noise variance at point s 1. is given by the scalar R11 ... In practice, one can evaluate Rii from a continuous function R(s) of the variables, i.e. Rii = R(s)I s=s.· For all cases 1 P0 was chosen to be 0.04 thus giving trace P0 = 0.56. the constant k1 in The constants a 1 and s1 in the optimality index were both chosen to be unity. The adjoint variable matrix x(t) is governed by a differential equation with final value specified. In order to change it to an initial value problem, we define and The matrix x(T) is governed by (3.76) 85 TT -1 - x(T)P(tf-T) [ N I a . (tf-T)E.H.R .. H.E. ] i =l l l l 11 l l TT - [ i=lIN a.(tf-T)E .H.R -1 .. H.E.] P(tf-T)X(T) l l l 11 l l (3.77) with the initial condition x( 0) = a. 1 (3.78) I The switching functions ai of (3.72) are (3.79) The advantage of expressing ai in the form (3.79) is that one need not store the matrix x(T). The covariance matrix P(t) is calculated and stored . The adjoint variable matrix x(T) is calculated but only the scalar integrand in the expression for ai need be stored for use in integration. For seven interior collocation points, the matrices P(t) and x(t) are governed by 14xl4 matrix differential equations. Taking the advantage of the sylTllletry property of these matrices, one needs to solve 105 nonlineardifferential equations in each case. A variety of computer programs are available for solving large systems of differential equations. The implementation of the Gear technique by Hindmarsh and Byrne (1975) was chosen because of the ability of the program to handle stiff differential 86 equations. This program uses a variable step, variable order method and employs linear multistep backward difference formulae. The optimal sen- sor locations are determined by using the algorithm of Section 3.5. 3.6.3 Effect of Varying Measurement Noise The effect of varying measurement noise was studied by determining optimal locations for a fixed Q and varying R(s). variance Q was taken to be 0.01 I. The process noise co- In the first two cases, R(s) was kept constant at 0.0004 and 0.004 respectively. In the third case R(s) was taken to be a function of s, namely R(s) = 10- 4 exp(+l0!s-0.51) The results for the optimal sensor location problem for these three cases are surrmarized in Table 3.2 which gives the initial and the optimal locations and corresponding values of the optimality index J. For the case of constant R(s), the optimal locations are governed mainly by the dynamics of the system and remain same in the two cases although the magnitude of R(s) is ten times higher in the second case. The magnitude of J is higher for the second case because of the larger uncertainty in measurements. In both cases, the optimal value of J is smaller compared to the initial value by a factor of 2 to 4. It is observed that the optimal locations for the temperature sensors are near the hot spot in the steady state temperature profile. This would seem to be quite natura.l because the un- certainty in temperature resultjng from the process noise is expected to be high near this hot spot. 87 Table 3.2 Effect of Measurement Noise P0 = 0.04 I Q = 0.01 I l 2 3 Initial Optimal No. of IteraLoca- LocaR(s) tions tions . tions Initial J -2 0.0004 1 11. 78xl 0 s2,s3 S4,S6 0.004 s ,s6 l 27 .46x10 -2 s2,s3 4 14. l 9xl0- 2 .0001 exp(+lOls-.51) 2 S3,S6 S3,S4 Optimal J 3.42xlo- 2 10.55xlo- 2 11. OOxlD- 2 . - a • f-- c:c 0:::: . . a l~ \ R(S) = 0.0004 0 .10 0.20 o.3o o.~o o.so o.so TIME FIGURE 3,3: TRACE PAS A FUNCTION OF TIME FOR TWO SENSORS LOCATED AT (A) POSITIONS S4 AND S6 (B) OPTIMAL POSITIONS S2 AND S3 o.o a.___ _ _____._ _ _ __.__ _ _ ___.__ _ _ __..__ _ _ __.__ _ _ __. ("\j a a I =:fl a a a w u • CD 0... a . a co .a a a ("\j CX> CX> I- c::c 0:: u w o.__ • o.o I 0.10 \ 0.20 o.3o o.qo RCS)= 0,004 o.so o.so TIME FIGURE 3.4: TRACE PAS A FUNCTION OF TIME FOR TWO SENSORS LOCATED AT CA) POSITIONS S4 AND S5 CB) OPTIMAL POSITIONS S2 AND S3 • 0 - - -- -- - - - - - - - - - - - - - - - - - - - - - - - ----- C) 0 0 Ln • C) - C) • -, 0 Ln (\J 8• (\J • Ln C) en • g co I.O . CJ 0 I- 0::: <C w u I II \ \ R(S) = 0,0001 EXP(+lOIS-0.51) 0.10 0.20 0.-30 0.1!0 0.50 0.60 TIME FIGURE 3.5: TRACE PAS A FUNCTION OF TIME FOR TWO SENSORS LOCATED AT (A) POSITIONS s3 AND s6 (B) OPTIMAL POSITIONS S3 AND S4 o.o lo CJ'--"-_ _ _ _ _ _ _ _ _ __...__ _ _____.___ _ _-'--_ ______J CJ c-- N Lf) C lo LJ') CJ r-. LJ') 0... _.lo 8 .-.• N LJ') _. LJ') I.O 0 91 For the third case, the measurement noise variance assumes its lowest value near the center points= 0.5 and increases exponentially as we move towards either end. The optimal sensor locations change. One sensor is still located near the hot spot while the other moves to a position of low measurement uncertainty. The plots of trace P vs. time for the initial locations and the optimal measurement locations are shown in Figures 3.3, 3.4 and 3.5. 3.6.4 Effect of System Noise The measurement noise variance R(s) was kept fixed and the optimal sensor locations were determined for three covariance matrices Q. first two cases Q was taken as 0.01 I and 0. 1 I, respectively. In the As dis- cussed before, these choices of Q assume negligible correlation between the system model errors at different collocation points. In the third case these errors were assumed to be correlated and Q was determined from the matrix Q+(s,s') given by Q+(s,s') 0.01 exp(-lO !s-s'I) 0 (3.80) = 0 0.01 exp(-lO!s-s'I) This choice of Q+ assumes that there is no correlation between the disturbances wc(s,t) and wT(s,t). The results for the three cases are summarized in Table 3.3. The optimality index J is high in the second case because of high system noise covariance. The optimal locations remain unchanged in the three cases because the system noise Q is not a function of the spatial variables even though the magnitude of Q is varying. Figures 3.6 and 3.7 show the 92 Table 3.3. Effect of Process Noise P0 = 0.04 I R(s) = 0.0004 lniti al Optimal Loca- Locations tions Q No. of Iterations Initial J Optimal J 3.42xl0- 2 7. 13xl0-2 3.8lxl0 -2 l 0.01 I s4,s6 s2,s3 l 2 0. l s4,s6 s2,s3 l ll.78xl0 -2 27.82xl0 -2 3 Q evaluated from S4,S6 s2,s3 l 12.85xl0-2 I (3.80) I- a::: c:::c w u CL . C) • • C) C) . C) ('J . C) 0 :::t' C) C) • CD C) C) (X) . C) C - C) ...-4 ('J \ 0.30 0.20 0.25 a.a a.as 0.15 0.10 TIME FIGURE 3,6: TRACE PAS A FUNCTION OF TIME FOR TWO SENSORS LOCATED AT (A) POSITIONS S4 AND S5 (B) OPTIMAL POSITIONS S2 AND S3 ~ "' Q = 0.1 I \.0 w . - a <C I- 0:::: . • a 0 a . C\I a a a =1' \ Q+CSJS') = 0,01 EXPC-lOIS-S' I) 0. 05 0 .10 0. 15 0. 20 0. 25 TIME FIGURE 3.7: TRACE PAS A FUNCTION OF TIME FOR TWO SENSORS LOCATED AT CA) POSITIONS S4 AND S5 CB) OPTIMAL POSITIONS S2 AND S3 0. 0 0. 30 L-----------1------1---------+-----+------4-------' a~/ u w . CD o._ a . a co a .....• a a C\I ~ .j:::,, 95 plots of trace P vs. time for the second and third case, while the plot for the first case is shown in Figure 3.3 3.6.5 Number of Sensors Optimal locations for 1, 3, and 5 sensors were determined for the case of constant Q and R. The plots of trace P In each case only one iteration was required. vs. time for the initial and optimal sensor loca- tions are shown in Figures 3.8, 3.9, and 3.10. for the case of l, 2, 3, and 5 sensors. Table 3.4 shows the results From the plots it is observed that for optimally located sensors, trace Pis a monotonically decreasing function of time. This implies that the concentration and temperature estimates become more accurate as more temperature measurements become available for use in estimation. function of time and space. The errors in the state estimate are a The optimally located sensors are near the region where the errors are expected to be highest. However, when the sensors are not optimally located, the measurements cannot give accurate information about the state at points far away,and therefore the total uncertainty, measured by trace P , may increase at a particular time. That is why for nonoptimally located sensors trace Pis, in general, not a decreasing function of time. The question of number of sensors required can be answered from the data in Table 3.4. We plot trace P vs. time t for l, 2, 3, 5, and 7 optimally located sensors in Figure 3.11. From this Figure, we conclude that for the tubular reactor under consideration, very little improvement in the accuracy of the state estimate is possible by using more than two sensors. One optimally located sensor can also give good state estimates. 96 Table 3.4. Optimal Locations for One, Two, Three, and Five Sensors P 0 = 0.04 I Q = 0.01 I R(s) = 0.0004 Number of Initial Sensors Locations Optimal Locations No. of Iterations l s8 s2 l 2 S4,S5 s2,s3 l 3 S4,S5,S7 S2,S3,S4 l 5 S4,S5,S5 S2,S3,S4 l s7's 8 S5,S5 Initi al J Optimal J 36.57xl0 -2 11. 78xl0 -2 4.72xl0 -2 3.42xl0- 2 11. 28xl0 -2 11. l OxlO -2 3.20xl0 -2 3.llxlO -2 . . 0 0 0:::: I- . I I ,/ \ \ A 0. 10 0. 20 0. 30 0. l!O O. 50 0. 60 TIME FIGURE 3,8: TRACE PAS A FUNCTION OF TIME FOR ONE SENSOR LOCATED AT (A) POSITION s8 CB) OPTIMAL POSITION S2 0. 0 C) L__ _--=:t=========±==========±========:::::f========:±::::=========f C) - 0 ~ w u C) 0 c::t: C) o..._ ('I") C) :::r C) '.J ID a a f-- 0:::: <C a w u . • a a a . \ 0 .0 0 . 10 0 . 20 0 . 30 0 . l!O O. SO L - - - - - - L - - - - - - L - - - - - - - - L - - - - . . l - - - - - - " -- N J O. 60 -----' TIME FIGURE 3.9: TRACE PAS A FUNCTION OF TIME FOR THREE SENSORS LOCATED AT (A) POSITIONS S4J S5J S7 (B) OPTIMAL POSITIONS s2J s3 s4 . ('J a a =t' a CD CL a ro. a .....• a a ...... ('J I.O o:> . 0 .10 0. 20 0. 30 0. 40 0. 50 _ 0. 60 ____. TIME FIGURE 3.10: TRACE PAS A FUNCTION OF TIME FOR FIVE SENSORS LOCATED AT CA) POSITIONS s4J s5J s6J s7J s8 CB) OPTIMAL POSITIONS S2J S3J S4J S5J S5 0. 0 . L - - -- - - - - - l - - - - - ' - - - - - - - ' - - - - - _ . . _ __ _ __.___ C) a C) ("\I a :::r . a C) C) CD . C) C) co . C) I - C) C) - . C) ("\I I.O I.O • . . e:::: I- <C I I Q I Cl 0 0 __. 0 N Cl 0 w u . ("I") 0... 0 0 0 ::,, 0 LJ') 0 0 0 to 0 .10 0.20 o.30 o.40 a.so a.so TIME FIGURE 3.11: TRACE PAS A FUNCTION OF TIME FOR CA) ONE CB) TWO CC) THREE CD) FIVE~ AND CE) SEVEN OPTIMALLY LOCATED SENSORS 0.0 ___ _ _ _ _ _ _ _____.__ _ _ ___.___ _ _ _~ - - - - - - - - - - - - - ll \} 0 0 - 101 However, as is evident from the results in Table 3.4, even five sensors with arbitrary locations may not give an estimate as accurate as provided by just one optimally located sensor. 3.6.6 Accuracy of Concentration and Temperature Estimates The tubular reactor is a system with two state variables, concentration and temperature. The trace of the covariance matrix P(t) is a measure of the total uncertainty at time tin the estimates of concentration and temperature profiles. The uncertainties in the concentra- tion and temperature estimates are measured by trace Pc(t) and trace PT(t) respectively, where N = l p2i-l,2i-l(t) i=l and, (3.81) N = i=l I P2.l, 2.l ( t) Since the concentration and temperature equations are coupled, we can obtain good estimates of the concentration profiles even though only temperature is being measured. In Figures 3. 12-3.16, we plot trace P, trace Pc and trace PT vs. time for optimally located 1, 2, 3, 5 and 7 sensors. In each case all the three variables being plotted are decreas- ing functions of time. The concentration estimates also become better with time as more temperature measurements are made available. steady state, trace Pc is larger than trace PT. At the This is an expected result :E: I<C a:: >< ....... I- a::: <C w u . . I . . I a a a . ...... a a (\I a 0 en a a :::r a a U1 a a CD a 0.10 0.20 0.30 0.40 0.50 0.60 TIME FIGURE 3.12: (A) TRACE P (B) TRACE Pc AND (C) TRACE PT AS FUNCTIONS OF TIME FOR ONE OPTIMALLY LOCATED SENSOR 0.0 ~ ;\ \ I-' 0 N ~ I<C 0::::: ........ X I- 0::::: <C u w • . • 0 0 a C) - a (\/ 0 C) (T) C) C) :::r . C) a Lf) C) C) . 0 C.D 0. 10 0.20 0 , 30 0.4 0 0.50 0. 60 TIME FIGURE 3. 13: (A) TRACE P(B) TRACE PcJ AND (C) TRACE PT AS FUNCTIONS OF TIME FOR TWO OPTIMALLY LOCATED SENSORS 0 .0 l\ w 0 - . . a ~ c::::c I- • a a a . ....a a N a >< ........ c:::: a . ("') a a =1' • a l.n a a I- c:::: c::::c w u a CD O. 10 0. 20 0. 30 0. l!O O. 50 0. 60 TI M E FIGURE 3.14: (A) TRACE P (B) TRACE Pc (C) TRACE PT AS FUNCTIONS OF TIME FOR THREE OPTIMALLY LOCATED SENSORS 0. O L__ ~:::±==========f:::=========:::±::=:=========±=:========:±=========:=::! l\ ~ 0 ...... C) . a ~ 0::: Ic::I: I a 0 a 0 ~ \ 0 .10 0.20 o.3o o.4o a.so a.so TIME FIGURE 3.15: (A) TRACE p (B) TRACE Pc, AND (C) TRACE Pr AS FUNCTIONS OF TIME FOR FIVE OPTIMALLY LOCATED SENSORS 0.0 . L~---===t:==========±::=========±=========:±::===========±:========::::! -. a C'\I C) C) >< ...... . a en C) :::1' • a lJ') c::I: 0::: I- u w . C) CD a 0 1--1 u, . C) . CY) C) C) ::r :;:: c::c: 0::: C) CJ . C) C) -. C) • ("\J CJ .,._ CJ >< ....... 0::: .,._ u c::c: LLJ C) • C) l[) C) . C) (0 \., 0 .10 0.20 0.30 0.40 0.50 0.60 TIME FIGURE 3.16: (A) TRACE P (B) TRACE Pc1 AND (C) TRACE PT AS FUNCTIONS OF TIME FOR SEVEN SENSORS 0.0 I I "' ...... 0 107 since we are making only temperature measurements and estimating both temperature and concentration from these measurements. For all the cases in Fi gures 3. 12 - 3.16, the steady state trace Pc is about 4 to 5 times the steady state trace PT. 3.7 Summary In this Chapter, we have studied the problem of optimal sensor loca- tions for a non-adiabatic tubular reactor. The concentration and tem- perature equations are discretized by the orthogonal collocation method and the optimal sensor l ocat i ons are determined by employing the matrix minimum principle. A detailed example is considered and the effect of process noise and measurement noise on the optimal locations is studied. We have also looked at the question of number of sensors required for estimating the concentration and temperature profiles to a desired accuracy. CHAPTER 4 SUMMARY AND RECOMMENDATIONS 109 In this work, we have considered the prob l em of optimal measurement system design in distributed parameter systems . The measurement system design problem has been treated within the framework of estimation theory with the objective of estimating the system states as accurately as possible. An upper bound for the error covariance of state estimate has been developed. The optimal sensor locati on problem is then posed as a mini- mization problem for a scalar measure of this upper bound. The use of this upper bound enables us to simplify considerably the computational aspects of t he problem. This approach has been illustrated by an example involving one-dimensional heat conduction. The effect of process dynamics on optimal sensor locations has been clearly demonstrated, and it is shown that the state estimates can be considerably improved by proper choice of sensor locations. The optimal sensor location problem for a tubular chemical reactor has also been considered. Chemical reactors are an integral part of the process industry and the knowledge of concentration and temperature profiles is often required. The sensor location problem for this system is very interesting because it has two state variables, concentration and temperature, which are coupled through the mas s and energy balance equations. The optimal locations of sensors in this system have been deter- mined based on an in i tial discreti zation of the system equations. The accuracy of concentration and temperature estimates depends strongly on the sensor locations. For spatially invariant error statistics, the opti- mal sensor locations are governed mainly by process dynamics. The 110 feasibility of estimating both concentration and temperature profiles based on temperature measurements only has been investigated. It has been shown that the concentration profiles can indeed be estimated quite accurately on the basis of temperature measurements. The question of the num- ber of temperature sensors necessary for an accurate estimation of the profiles has also been investigated. Most systems of interest in the process industry are distributed in nature and are often described by more than one state variable. The opti- mal sensor locations for these processes can be determined by using the approach developed in Chapter 3. The equations describing the system can be discretized using a low-order discretization method and the sensor locations can be selected using the same approach. The approach developed here can be used for selecting optimal sensor locations in a general process system. The optimal control problem for a dynamic process is the dual of the optimal estimation problem. Therefore, the problem of optimal point- wise controls can be approached in a way similar to the optimal sensor location problem. It is hoped that the approach developed here will find potential applications in the optimal location of pointwise controls. In addition to the problem of optimal location of pointwise controls, some other interesting areas of research are the following: The problem of the number of sensors required for a given accuracy in state estimate should be further investigated. By properly assigning a value to improved estimate accuracy, and comparing it to the cost of adding more sensors, the cost-effective number of sensors can be determined. 111 Sensor location problem can be investigated based on other objectives. The location of sensors with the objective of 11 best" estimating certain spatially varying parameters in the systems is an important area of research. In this work, we have dealt with the case of measurements continuous in time. However, in certain situations continuous measure- ments may be very costly. For such systems, the determination of opt imal timing and freq uency of measurements should be investigated. 112 REFERENCES Aidarous S.E., (1976) Optimal Allocation Strate..[_ies in Stochastic Distributed Systems, Ph.D. Thesis, Universite Catholique de Louvain, Louvain, Belgium. Amouroux M. and Babary J. P. (1975) Optimal Pointwise Control for a Class of Distributed Parameter Systems. Proceedings of the IFAC 6th World Congress, Part I. Athans M. (1968) The Matrix Minimum Principle. Inf. Control, ll_, 592-606. Athans M. (1972) On the Determination of Optimal Costly Measurement Strategies for Linear Stochastic Systems. Automatica, ~' 397-412. Balakrishnan A. V. and Lions J. L. (1967) State Estimation for InfiniteDimensional Systems. J. of Computer and System Science, 1_, 391-403. Balakrishnan, A. V. (1976) Applied Functional Analysis, SpringerVerlag, New York. Bensoussan A. (1972) ' Optimization of Sensors' Location in a Distri- buted Filtering Problem,in Stability of Stochastic Dynamical Systems, Springer-Verlag, New York, 62-84. Brewer J. W. and Moore S. F. (1974) Monitoring: An Environmental State Estimation Problem. J. Dynamic Systems Measurement and Control, 96, series G, 363-365. ; Brogan W. L. (1968) Optimal Control Theory Applied to Systems Described by Partial Differential Equations,in Advan. in Control Systems,§_, 221-316. Chen W. H. and Seinfeld J. H. (1975) Optimal .Location of Process Measurements. Int. J. Control, 11, 1003-1014. 113 Ciarlet P. G., Schultz M. H. and Varga R. S. (1967) Numerical Methods of High-Order Accuracy for Nonlinear Boundary Value Problems. Numerische Mathematik, 1, 394-430. Courant R. and Hilbert D. (1966) Methods of Mathematical Physics, Vols. I and II, Interscience, New York. Falb P. L. (1967) Infinite-Dimensional Filtering: Filter in Hilbert Space. The Kalman-Bucy Inf. Control,..!.!_, 102-137. Finlayson B. A. (1972) The Method of Weighted Residuals and Variational Principles, Academic Press, New York. Fletcher R. and Powell M. J. D. (1963) A Rapidly Convergent Descent Method for Minimization. Comput. J., §_, 163-168. Georgakis C, Aris R. and Amundson N. R. (1975) Modal Control of Unstable Reactors. Presented at the 68th Annual Meeting of the AIChE, Los Angeles. Georgakis C., Aris R. and Amundson N. R. (1976) Stabilization of Unstable Nonlinear Distributed Chemical Reactors by Observers Design, Proceedings of the Joint Automatic Control Conference, 161-170. Herring K. D., and Melsa J. L. (1974) Optimum Measurements for Estimation. I.E. E. E. Trans. Auto. Control, AC-19, 264-266. Hille E. and Phillips R. S. (1957) Functional Analysis and Semi-Groups. Am. Math. Soc. Providence, Rhode Island. Hindmarsh A. C. and Byrne G. D. (1975) EPISODE: An Experimental Package for the Integration of Systems of Ordinary Differential Equations, LLL, Univ. of California, Report UCID-30112. 114 Hwang M., Seinfeld J. H. and Gavalas G. R. (1972) Optimal Least Square Filtering and Interpolation in Distributed Systems. J. Math Anal. Appl., 39, 49-74. Jazwinski A. H. (1970) Stochastic Processes and Filtering Theory, Academic Press, New York. Johnson C. D. (1969) Optimization of a Certain Quality of Complete Controllability and Observability for Linear Dynamic Systems. J. of Basic Eng., 91, series D, 228-238 . . Kalman R. E. and Bucy R. S. (1961) New Results in Linear Filtering and Prediction Theory. J. of Basic Eng., 83, series D, 95-108. Kushner H.J. (1970) Filtering for Linear Distributed Parameter Systems. SIAM J. Control, 8, 346-359. Meditch J. S. (1970) On State Estimation for Distributed Parameter Systems. J. of the Franklin Institute, 290, 49-59. Mehra R. K. (1976) Optimization of Measurement Schedules and Sensor Designs for Linear Dynamic Systems. I. E. E. E. Trans. Auto. Control, AC-21, 55-64. Morse and Feshbach (1953) Methods of Mathematical Physics, Parts I and II, McGraw-Hill, New York. Muller P. C. and Weber H. I. (1972) Analysis and Optimization of Certain Qualities of Controllability and Observability for Linear Dynamical Systems. Automatica, _§_, 237-246. 115 Pimentel K. D. (1975) Monitoring: Toward a Mathematical Theory of Environmental The Infrequent Sampling Problem. Ph.D. Thesis, LLL, Univ. of California, Livermore, CA. Powell M. J. D. (1970) A Hybrid Method for Nonlinear Equations,in Numerical Methods for Nonlinear Algebraic Equations, Ed. P. Rabinowitz, Gordon and Breach, New York, 87-114. Sakawa Y. (1972) Systems. Optimal Filtering in Linear Distributed-Parameter Int. J. Control, J&, 115-127. Seinfeld J. H., Gavalas G. R. and Hwang M. (1971) in Distributed Parameter Systems. Nonlinear Filtering J. Dynamic Systems, Measurement and Control, 93, series G, 157-163. Seinfeld J. H., and Koda M. (1977) Numerical Implementation of Distri- buted Parameter Filters with Application to Problems in Air Pollution. Lecture Notes in Control and Information Sciences: Working Conference on Distributed Parameter Systems Modeling and Identification, Ed. A. Ruberti, Springer-Verlag, Berlin, in press. Sorenson H. W. (1968) Controllability and Observability of Linear, Stochastic, Time-Discrete Control Systems in Advan. in Control Systems, §_, 95-158. Stakgold I. (1968) Boundary Value Problems of Mathematical Physics, Vols. I and II, Macmillan, New York. Thau F. E. (1969) On Optimum Filtering for a Class of Linear Distributed- Parameter Systems. J. of Basic Eng., 91, 173-178. 116 Tzafestas S. G. and Nightingale J. M. (1968a) Optimal Filtering, Smoothing and Prediction in Linear-Distributed Parameter Systems. Proc. IEE , 115, 1207-1212. Tzafestas S. G. and Nightingale J.M. (1968b) Concerning Optimal Filtering Theory for Distributed Parameter Systems. Proc. IEE , 115, 1737-42 Tzafestas S. G. and Nightingale, J.M. (1969) Maximum-Likelihood Approach to Optimal Filtering of Distributed-Parameter Systems. Proc. IEE, 116, 1085-1093. Varma A. and Amundson N. R. (1972) Adiabatic Tubular Reactor. Some Problems Concerning the Non- Can. J. Chem~ Eng., 50, 470-485. Varma A. and Amundson N. R. (1973a) Some Observations on Uniqueness and Multiplicity of Steady States in Non-Adiabatic Chemically Reacting Systems. Can. J. Chem. Eng., g, 206-226. Varma A. and Amundson N. R. (1973b) Stability Considerations. Villadsen J. (1970) ing Problems. The Non-Adiabatic Tubular Reactor: Can. J. Chem. Eng., g, 459-467. Selected Approximation Methods for Chemical Engineer- Inst. for Kemiteknik D.T.H., Copenhagen. Wang P.K.C. {1964) Control of Distributed Parameter Systems in Advan. in Control Systems, l, 75-172. Yosida K. (1966) Functional Analysis, Springer-Verlag, New York. Yu T. K. and Seinfeld J. H. (1973) Observability and Optimal Measurement Location in Linear Distributed Parameter Systems. Int. J. Control,~, 785-799. 117 Yu T. K., Seinfeld J. H. and Ray W. H. (1974) Time Delay Systems. Filtering in Nonlinear I.E.E.E. Trans. on Auto. Control, AC-19, 324-333. 118 APPENDIX A This appendix gives the proof of the inequality (2.53). The optimal filter for the process with spatially independent integral measurements is given by (2.32) and (2.33), where u{x,t) is the optimal estimate of the state u(x,t), and P(x,y,t) is the error covariance matrix. Let us n<M consider two processes described by the following equations: (A.1) (Process I) z1(t) = J H(x,t)u 1(x,t) dx (A.2) 0 I. C. u1 (x,t ) = 0 B.C. ax u1(x,t) = 0 (A.3) 0 au II (x,t) at x e ao . (A.4) = A uII(x t) X (A.5) ' (Process II) z11 (t) = J H(x,t)u 11 (x,t) dx + v(t) (A.6) D I.C. E[u 11 (x,t 0 )] = u0 (x) E[{u 11 (x,t )-u (x)}{u 11 (y,t 0 )-u 0 {y)}T] = P0 (x,y) 0 B.C. 0 ax uII(x ' t) = O xe ao {A.7) (A.8) 119 We note that u(x,t) = uI (x,t) + uII {x,t) (A.9) and process II is simply the original process with no dynamic noise. Let us now define the suboptimal estimates GI(x,t) and Qlt(x,t) for the two processes -, utilizing the gain matri _x for the optimal estimate u(x,t) of the original system. The suboptimal estimates for the two processes are governed by AI au (x,t) = A uI(x t) at X ' (A. 11) I.C. xe ao B.C. (A.12) {A. 13) I. c. u11 (x,to) = uo(x) B.C. Cl., uII ( x, t) = 0 (A.14) xe ao (A. 15) 120 Note that in equations (A.10) and (A.13), P(x,x',t) is the error covarian~e for the optimal estimate u(x,t) in (2.32). Th~ optimal estimate u(x,t) of the original process and the s~boptimal estimates uI(x,t) and uII(x,t) of the processes I and II obey the relationship u(x,t) = uI(x,t) + uII(x,t). ~ ) u(x,t) = ~I u (x,t) + ~II u (x,t. In terms of the estimate errors, we have The covariance matrix P{x,y,t) is, by defi ni ti on P(x,y,t) ~ E[u(x,t)uT(y,t)] I II ~I ~II (y,t)} T] = P (x,y,t) + P (x,y,t) + E[{u (x,t)}{u ( A. 16) where P1{x,y,t) and PI 1{x,y,t) are the error covariance matrices for the estimates QI(x,t) and oII(x,t) of the states uI(x,t) and uII(x,t), • ) T] and respectively. The cross-variance terms, E[ {u~I( x,t ) }{u~II( y,t} 11 E[{u (x,t)Hii(y,t)}T] are zero each because the two processes are forced by uncorrelated white-noise processes and have uncorrelated initial conditions. Therefore P(x,y,t) = PI (x,y,t) + PII (x,y,t) (A. 17) If we use a suboptimal gain to obtain an estimate for the state u(x,t) with corresponding error covariance matrix Ps(x,y,t), we must have P(x,y,t) ~ Ps(x,y,t). Using the same suboptimal gain for obtaining estimates of the states uI(x,t) and u1I(x,t) with error covariances 121 pis(x,y,t) a~d PIIs(x,y,t), we shall have the relation (A.18) We choose the suboptimal gain to be the gain that yields an optimal estimate for the process II, and use the same gain for state estimation in process I. Let us denote the error covariance P11 s(x,y,t) for the optimal estimate of process II by P(x,y,t). The error covariance matrix pis ( x ,y, t) for process I cannot be 1arger than the error covariance W(x,y,t) for an estimate obtained by using zero gain. The matrix W(x,y,t) is governed by the equation aw(~ty,t) = Ax W{x,y,t) + W{x,y,t) A~+ B(x,t)Q(x,y,t)BT(y,t) (A. 19) I.C. W(x,y,t 0 ) = 0 B.C. ax W{x,y,t) = O x e ao y E DxaD W(x,y,t) a;= 0 ye ao xe DxaD {A.20) {A.21) It can be easily verified that the solution to (A. 19)-(A.21) is given by (2.54). The inequality (2.53) is thus proved. 122 APPENDIX B In this appendix we verify the validity of (2.58) as a solution of the equations (2.55)-(2.57). From equation (2.58) P(x,y,t 0 ) = II G(x,t ;x' ,t )[P~(x 1 ,y 1 ) + M(x' ,y' ,t )JtGT(y,t ;y' ,t DD 0 0 0 0 ) dy 1 dx 1 0 (B.l) Using the fact that M(x,y,t ) = 0, and G(x,t ;x',t 0 get P(x,y,t0 ) = P0 (x,y) . 0 ) = I o(x-x'), we 0 The initial condition (2.56) is thus verified. Boundary conditions (2.57) can be easily verified by using the properties of the Green's function matrix G(x,t;x',t'). To verify that the l'"(x,y,t) given by (2.58) obeys the partial differential equation (2.55) we shall need, in addition to the properties of the Green ' s function matrix, the folla.t1ing identity aAt(xaty ' t) = - II At (x,x' ,t) aA(x' ,y',t)- At (y' ,y,t) dy'dx' - at DD (B. 2) Di fferentiating (2.58) with respect tot, and using the properties of the Green's function matrix, we have + II G(x,t;x' ,to) ;t [P~(x',y') + M(x',y',tl] t GT(y,t;y ' ,to) dy'dz' DD Using the identity (B.2), the last term on the r.h.s. of (B.3) is (8.3) 123 =---u G(x,t;x',to) u [P~(x',x") + M{x',x",t)]t C.Jl GT(n;,t;x",to) Hi(t) •(R-\t));j H/t)G(nj,t;y",to)] [P~(y 11 ,y 1 ) + M(y 11 ,y 1 ,t)]t GT(y,t;y',to) dy 11 dx 11 dy 1 dx 1 •{Hi{t){R-l{t));j Hj(t)}{JI G(nj't;y",t 0 ) [P~(y'',y') + M(y",y',t)r DD •GT(y,t;y',t ) dy"dy'} 0 =-i j 1 { P(x,ni' tl}{ Hi {t)(R-\t)) ij Hj (tl}{ P(nj ,y, t)} Thus we get aP(x,y,t) at ) -( T = A~( x x,y,t + P x,y,t) AY - m T 1 L ~(x,n;,t)H.(t)(R(t)) .. 1 i,j=l lJ (B. 4) 124 The expression (2.58) for P(x,y, t ) thus obeys the partial differential equation (2.55) with initial condition and boundary conditions given by (2.56) and (2.57). 125 APPENDIX C In this appendix, we derive the discrete equations (3.37) - (3.38) for the tubular reactor. As expressed in (3.35) the first and second deriva- tives of concentration at the collocation points si, i = 1,2, ... , N' where N' = N+2, are given by N' ac = I as s. j=l , A.. C. ( t) lJ J (C.1) N' 2 a c I as7 s. = j=l 8 .. cj ( t) lJ 1 The boundary conditions at s=O and s=l provide the equations N' I j=l N' I j=l Alj c.J = pem (c,- n (C.2) AN'j C.J = 0 (C.3) Equations (C.2) and (C.3) allow us to express the concentrations c1 and cN' at the boundaries in terms of the concentrations at the interior collocation points. (C.4) (C.5) 126 where Using {C. 1) and {3.8), we find that the concentrations cj at the interior collocation points are governed by the equations de. N' N' crt1= I BJ.k ck - P I k=l em k=l A.k ck - ar{'c.,T.) J J J (C.6) {j = 1,2, ... , N') The equations (C.4) and (C.5) can be used to eliminate c1 and cN'' Thus, we get de. N+l _J_dt = I D.k ck - ar(c.,T.) + p. k=2 J J J J ( j = 2 , 3 , ... , {C. 7) N+ 1) where D.k = (B.k - p A.k + n,w.k) J J em J J (C.8) wjk = (Bjl - pem Ajl)(AlN' AN'k - AN'N' ' A,k) + (BjN' - pem AjN')[AN'l Alk - (All - pem) AN'k] 127 The equations for temperature T. at the interior collocation points can J similarly be obtained. dT. N+l Le _J_dt = ' k~2 We get S.k Tk + aSr(c.,T.) + y(T -T.) + q. J J J a J J (C.9) (j = 2,3, ... , N+l) where s.k = B.k - p A.k + n2 W' "k J J eh J J ( C. 10) W'jk = (Bjl - peh Ajl)(AlN' AN'k - AN'N' Alk) + (BjN' - peh AjN') [AN'l Alk - (All - Peh) AN'k] and for the case of Le= 1, we get dT. N+l --1.dt = I k=2 s.k Tk + a.Sr(c.,T.) + y(T -T.) + q. J J J a J J In fact for this case Sjk = Djk and qj = pj. (C. 11)