Boundary identification and control of distributed parameter systems using singular functions

Chemical Engineering Science 54 (1999) 1181—1204

Boundary identification and control of distributed parameter systems using singular functions Shrikar Chakravarti, W. Harmon Ray* Department of Chemical Engineering,1415 Engineering Drive, University of Wisconsin-Madison, Madison, WI 53706, USA Received 19 September 1997; accepted 6 September 1998

Abstract Singular functions are fundamental properties that capture the spatial nature of distributed parameter systems (DPS). This paper presents a framework for the dynamic identification of singular functions of an unknown DPS. The input to a DPS can be spatially distributed or applied at the boundary. Previous work [Gay and Ray, Chemical Engineering Science 50 (10) (1995) 1519—1539] considered the former. The current effort is more concerned with the case of boundary inputs. The concept of boundary topos is also introduced. Both singular functions and boundary topos can be identified through simple boundary input response tests and provide alternate bases for a pseudo-modal control scheme. Thus, this work also provides a systematic means for identification and control of distributed parameter systems through the use of boundary inputs alone. Examples of self-adjoint and non-self-adjoint systems are considered to demonstrate the efficacy of the identification and control schemes.  1999 Elsevier Science Ltd. All rights reserved.

1. Introduction The work presented here is aimed at developing a systematic framework for the boundary identification and control of linear distributed parameter systems (DPS) using empirically determined singular functions. A DPS is characterized by spatially varying outputs, states and parameters. The input, i.e. the control variable, can be distributed along the spatial domain or applied at the boundary. Examples of distributed parameter systems in chemical engineering include tubular or packed bed reactors, metallurgical heating processes, fluidized bed reactors, paper machines and polymer film extrusion processes. Singular functions are a generalization of the idea of eigenfunctions and were first applied to process identification and control problems by Gay (1989) and Gay and Ray (1988, 1995). In order to capture the spatio-temporal nature of distributed parameter systems, a convenient form of representation has been the medium of partial differential equations (PDE). If one possesses knowledge of the exact

*Corresponding author. Tel.: 001 608 263 4732; fax: 001 608 262 0832; e-mail: [email protected].

equations, boundary conditions and parameters that constitute the conventional PDE model, then a variety of methods exist for DPS control system design. In the absence of such exact model information, one has to resort to process identification methods. Hence, the challenge lies in identifying an unknown DPS from straightforward input/output experiments and then designing a controller that would take into account the spatially varying nature of the system. This has been the motivating theme for earlier research work in our group (Gay and Ray, 1988, 1995; Chakravarti and Ray, 1997). Some significant progress in this regard was made by Gay and Ray (1988, 1995) who showed that using an integral equation model facilitated viewing the distributed parameter system in terms of an input/output framework. Use of empirically determined singular functions helped account for the spatial nature of the problem. 1.1. Boundary Identification In the previous work (Gay and Ray, 1995), spatially varying inputs were used. In the simplest situation, by varying the spatial form of the inputs, appropriate steadystate experiments sufficed to identify the steady-state singular functions. Dynamic spatially varying inputs can also

0009-2509/99/$ — see front matter  1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 9 8 ) 0 0 5 0 0 - 4

1182

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

be used to provide more information and to identify the non steady-state singular functions (Gay, 1989; Gay and Ray, 1988). In this work, boundary inputs are considered. For this case, since steady-state experiments yield limited information, it becomes imperative to incorporate dynamic information. This necessitates reworking the methodology that had been developed for spatial inputs. The primary objective of this paper is to develop a framework where, by suitable design of dynamic experiments, the singular functions of an unknown DPS can be determined. The approach adopted has been to first address the general case of spatially varying inputs and then specialize for the case of boundary inputs. The class of DPS addressed include those linear systems for which a Green’s function exists. It is also assumed that the distributed parameter systems exhibit first order dependence with respect to time. No detailed knowledge of the underlying PDEs is needed. In our examples here, the models merely serve as surrogates to provide “data” for the actual distributed parameter processes. 1.2. Karhunen—Loe` ve expansion In developing the dynamic identification framework, use of the Karhunen—Loe`ve (K-L) expansion technique is also explored (Sirovich 1987; Fukunaga, 1990). The computational approach is described in Appendix C. K-L decomposition is primarily used in the analysis of turbulence or spatiotemporal patterns. From the spatiotemporal pattern, a symmetric spatial correlation operator with spatially dependent eigenfunctions (or topos) and a symmetric temporal correlation operator with timedependent eigenfunctions (or chronos) are constructed. The decomposition is biorthogonal in the sense that both the topos and chronos form complete orthogonal bases for the respective spaces. Where necessary, we use the terminology boundary topos to reflect the fact that the topos are identified from the transient profiles arising in response to an input at the boundary. We also make the following distinction between the singular functions and boundary topos. While both sets of basis functions can be computed from empirical data, only the singular functions are fundamental properties of the DPS. The boundary topos are purely empirical; however, they can still be used to good effect in the identification and control of distributed parameter systems with boundary inputs. Applications of the K-L procedure can be found in a variety of areas such as fluid mechanics, pattern-recognition, statistics, signal analysis, meteorology, data compression and catalysis. The K—L decomposition is also referred to by other names such as proper orthogonal decomposition (Berkooz et al., 1993; Graham et al., 1993; Lumley, 1970; Sirovich, 1987), biorthogonal decomposition (Aubry et al., 1991) or method of empirical orthogonal functions (Lorenz, 1956). In chemical engineering,

Graham et al. (1993) and Chen et al. (1993) have applied the K-L expansion technique for analysis of spatiotemporal patterns on catalytic surfaces. Chen and Chang (1992), who presented an identification and control methodology for distributed systems with unknown nonlinear dynamics, used the K-L decomposition to identify the topos from the output. More recently, Park and Cho (1996) have employed the K-L decomposition to obtain a set of basis functions that could be used in a Galerkin procedure for obtaining a good low-dimensional model of the DPS. Zhou et al. (1996) have demonstrated the use of the K-L expansion in conjunction with feedforward neural networks for purposes of developing an empirical model for a fixed-bed reactor. Also, Rigopoulos et al. (1997) have employed the K-L expansion in identifying the significant features of the profiles in continuous sheet forming processes. 1.3. Boundary control There is mathematical literature that deals with the underlying theory of problems involving boundary control of distributed parameter systems (e.g. Fattorini, 1968; Zolesio, 1991). Many of the applications are simple examples such as the classic pure heat conduction problem (Olmstead and Schmitendorf, 1997; Olmstead, 1980); however, Olivei (1974) designed a boundary modal control scheme for regulating the temperature of the melt during crystal-growth in crucibles. Also, Hanczyc and Palazoglu (1996) have applied sliding mode control theory to a class of distributed parameter systems where the application considered was the control of an isothermal tubular reactor. There is, however, a dearth of literature that provides a systematic framework for the identification of an unknown DPS and subsequent control using only boundary inputs. The current effort tries to fill this void by showing how one can use the empirically identified singular functions as the basis for a pseudo-modal boundary control scheme to achieve the desired control objective for even an unknown or partially known DPS. It is important to realize that the range of control problems addressed with the use of boundary inputs can be significantly different from the ones with spatial inputs. 1.4. Organization of the paper The organization of this paper is as follows. The PDE and integral equation models for distributed parameter systems are presented in Section 2. Also included is a brief discussion of the singular value theory of integral equations. The properties of the singular functions and issues relating to their identification when using spatial or boundary inputs are also elucidated. Section 3 describes the dynamic framework for the identification of singular functions of unknown DPS. Examples in

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

Section 4 illustrate the performance of the identification methodology when applied to linear self-adjoint and non-self-adjoint distributed parameter systems with boundary inputs. Section 5 introduces the concept of boundary topos. Section 6 demonstrates how the singular functions or boundary topos can be used in the development of a pseudo-modal feedback boundary control scheme. Examples of boundary control applications are presented in Section 7. The important contributions of the paper are summarized in Section 8.

2. Theory 2.1. Model forms In order to demonstrate the theory in exact fashion, we will choose a linear system for illustration. However, it should be noted that some of the underlying ideas in the identification and control methodology have also been shown to be applicable to nonlinear processes (Gay and Ray, 1988; Chakravarti et al., 1995). We consider distributed parameter systems which can be represented by linear, second order, parabolic partial differential equations with constant coefficients: *y *y *y !a #l #by(z, t)"0 *t *x *z

(1)

*y !L y"0. X *t

while a large value of Peclet number as compared to the Nusselt number would indicate that the system is highly non-self-adjoint. Different systems are studied by varying Pe, Nu and the boundary conditions. The PDE model for a DPS with boundary inputs can always be recast into an integral equation of the following type:



y(z, t)"

(2)

The equation is defined on the dimensionless spatial interval 0)z)1 for time t*0. Homogeneous initial conditions (y(z, 0)"0) are assumed without loss of generality. The input or control term is applied at the boundary and thus appears in the boundary conditions (described by operators B and B ):   B y(z, t)"f (t), (3a)  B y(z, t)"0. (3b)  The subscripts 0 and 1 refer to z"0 and z"1 respectively. Notice that the boundary input is applied only at the entrance, z"0, as in most flow problems. It is assumed that a'0. Such a PDE model can be used to describe a one-dimensional heat transfer in a thin metal plate, in which case y(z, t) refers to temperature. f (t) could correspond to heat flux or temperature at z"0 as per the boundary condition. a corresponds to the thermal diffusivity (k/oC ) scaled by the characteristic length (for N example, plate length) ¸, l to velocity and b to the external heat transfer coefficient (h/oC ). The dimensionN less quantities Nusselt number and Peclet number are defined as Nu"b/a and Pe"l/a, respectively. Pe"0 would indicate that the DPS is symmetric or self-adjoint

R

g(z, t!q) f (q) dq, (4)  where g(z, t), the boundary Green’s function, can be determined for any unknown DPS by use of an impulse or a step input. The step input, f (t)"1, (usually easier to implement) leads to a corresponding step response, S(z, t). From this, the impulse response may be determined as follows: *S(z, t) g(z, t)" . *t

(5)

Having identified g(z, t), the response of the system to any arbitrary time-varying input can be determined using equation Eq. (4). Eq. (4) can be derived directly from the PDE formulation or alternately obtained as a simplification of the following integral equation model for spatially distributed inputs:



 R

k(z, f, t!q)u(f, q) dqdf (6)   through careful use of the Dirac delta function (Butkovskiy, 1982). Such an integral equation formulation (Gay and Ray, 1995) is more general than the PDE model since the distributed input and output functions, u(z, t) and y(z, t), only need to be L (i.e. square integrable) on the region +z, t: 0)z)1, t*0,. The kernel k(z, f, t!q) is also L. The integral equation can also be expressed in compact operator notation as: y"Ku. The Laplace transform of Eq. (6) yields the following transfer function representation for the DPS:

y(z, t)"

and in operator notation as

1183





k(z, f, s)u(f, s) df, (7)  where s, a complex parameter, denotes the transform domain. Eq. (7) is a Fredholm integral equation of the first kind, with an L kernel k(z, f, s). The kernel is also referred to as the Green’s kernel transfer function and can be interpreted as the distributed analog of the transfer function for lumped systems. y(z, s)"

2.2. Singular value theory of integral equations For more details on the material presented in this section, the reader is referred to the dissertation of Gay (1989). The L kernel k(z, f, s) generates the linear compact integral operator K so that Eq. (7) may be written

1184

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

as: y(z, s)"Ku(f, s). Typically the kernel is unsymmetric and complex. The adjoint integral operator K* is characterized by the adjoint kernel k*(z, f, s)"k(f,z,s), where the overbar denotes the complex conjugate. The operators K and K* are related through the L inner product as shown: (Ku, y)"(u, K*y)

(8)

for any two L functions u and y. The simplest definition of inner product with a weighting function of unity is used. Now K"K* implies that K is self-adjoint or symmetric or Hermitian. If KK*"K*K, then K is referred to as a normal operator. Obviously, the class of symmetric operators is a subset of the class of normal operators. Both symmetric and normal operators possess a convergent bilinear expansion for the kernel in terms of a complete set of orthonormal eigenfunctions and their corresponding eigenvalues. For a symmetric operator, the eigenvalues are real and distinct, while for a normal operator they may be complex. The mathematical theory dealing with the properties of symmetric and normal operators is well established and can be found in most texts on functional analysis (see, e.g. Naylor and Sell, 1982; Stakgold, 1979; Ramkrishna and Amundsen, 1985). 2.2.1. Singular functions and singular values In order to identify an unknown DPS, a theory is needed which can extend the useful properties of symmetric operators (distinct eigenvalues, orthonormal eigenfunctions, convergent expansions) to the more general case of nonnormal operators. Such a theory does exist in the mathematical literature and is based on the singular value analysis of the kernel of the compact integral operator K. This theory was first given by E. Schmidt in 1907 and extended to real general kernels by Smithies in 1937. The books by Smithies (1958) and Cochran (1972) provide a detailed account of this theory. The key results are indicated below. The more general integral equation model [see Eq. (7)] for a linear DPS assumes only that the kernel of K is L. For this more general case of a nonnormal operator (and so an arbitrary kernel), a set of eigenfunctions may or may not exist and even if it does, it may not be orthonormal in the sense of the usual L inner product. Further the eigenfunctions may fail to provide a convergent expansion. However from the unsymmetric kernel of the nonnormal operator, two Hermitian kernels may be constructed as follows:

eigenfunctions. If the sequence of functions +w (z, s), deG notes the eigenfunctions of the kernel kk*(z, f, s) and +p (s), the corresponding eigenvalues, then by definition: G



p (s)w (z, s)" G G



kk*(z, f, s)w (f, s) df. (11) G  A corresponding sequence of functions +v (z, s), can be G defined in terms of the adjoint kernel k(f, z, s) and the eigenfunctions and eigenvalues of Eq. (9):



p (s)v (z, s)" G G



(12) k(f, z, s)w (f, s) df. G  It can be easily shown that the functions +v (z, s), are G actually the eigenfunctions of the kernel k*k(z, f, s) with +p (s), being the corresponding eigenvalues. For comG pactness, we use the operator notation: pw "KK*w , (13a) G G G p v "K*w , (13b) G G G 1 1 K*Kv " K*KK*w " pK*w "pv . (13c) G p G p G G G G G G Substituting for K*w in Eq. (13a) from Eq. (13b) yields G the analog of Eq. (12) which defines +w (z, s), in terms of G the kernel k(z, f, s) and +v (z, s),: G  p (s)w (z, s)" k(z, f, s)v (f, s) df. (14) G G G  The elements of the sequences +w (z, s), and +v (z, s), are G G known as the left and right singular functions of Eq. (7), respectively, and the corresponding +p (s), are referred to G as the singular values. The singular functions satisfy the following orthonormality relations:



 

 w (z, s)w (z, s) dz"d , (15a) G H GH   (15b) v (z, s)v (z, s) dz"d , GH G H  where the Kronecker delta is defined in the usual manner as:



1 if i"j, d " GH 0 if iOj.

(9)

By convention, the singular values are always positive or zero and are arranged in descending order of magnitude:

k(m,z, s)k(m, f, s) dm. (10)  The symmetry implies that each of the kernels possesses a set of distinct eigenvalues and associated orthonormal

p (s)*p (s)* 2 *0   The left and right singular functions form orthonormal bases for the ranges of the operators K and K* respectively. Of course, if K is self-adjoint (e.g. is the diffusive operator in the heat conduction problem with Dirichlet

 

kk*(z, f, s)"

k*k(z, f, s)"



 

k(z, m, s)k(f,m,s) dm,

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

boundary conditions), then the left and right singular functions are identical and are simply the eigenfunctions themselves. Also they are independent of the complex parameter s. The eigenvalues of K or K* are then j (s)"p (s). G G For the case of K being a normal operator, there again exists a simple relationship between its eigenfunctions

(z, s) and eigenvalues j (s) and the singular functions G G and singular values: w (z, s)" (z, s), (16a) G G v (z, s)" (z, s) sgn j (s), (16b) G G G p (s)""j (s)", (16c) G G where the symbol sgn denotes the signum function and is defined as follows: j sgn j" , (jO0), sgn 0"0. "j" " ) " denotes magnitude of the complex number. Smithies (1958) provides a rigorous proof of the above equivalence. However, for the general case of a non-normal operator, there do not seem to exist in the mathematical literature results that analytically link the eigenfunctions and singular functions for the corresponding kernel. As mentioned previously, this could be because there may not exist a set of eigenfunctions. However, such a general L kernel will always possess a complete set of orthonormal singular functions with distinct, real singular values which provide a bilinear expansion for the kernel that converges in the mean-square sense. Thus singular value analysis of non-self-adjoint systems is the natural generalization of eigenvalue analysis of self-adjoint systems. Unfortunately, the task of deriving analytic expressions for the singular functions is cumbersome which makes numerical computation of singular functions a more viable alternative. The bilinear expansion for the kernel k(z, f, s) which can be expressed as  k(z, f, s)" w (z, s)p (s)v (f, s) (17) G G G G is essentially the infinite-dimensional analog of the matrix singular value decomposition K"WRV* discussed in Appendix A. It is indeed possible to construct other bilinear expansions. For example, one could use the eigenfunctions and adjoint eigenfunctions of K (if they exist). However, by virtue of the Approximation ¹heorem due to Schmidt (Stewart, 1993), for a given finite truncation order N the bilinear expansion in terms of the singular functions provides the best mean square approximation to the kernel. A proof of this can also be found in Butkovskiy’s book (Butkovskiy, 1969).

1185

Substituting Eq. (17) in Eq. (7), one obtains



  y(z, s)" w (z, s)p (s) v (f, s)u(f, s) df. (18) G G G  G Thus the output y(z, s) can be expanded in terms of the left singular functions +w (z, s),. Similarly, the input G u(z, s) has a convergent expansion in terms of the right singular functions (using the corresponding bilinear expansion for the adjoint kernel):



  u(z, s)" v (z, s)p (s) w (f, s)y(f, s) df. G G G  G

(19)

2.3. Steady-state singular functions As noted above, for self-adjoint systems the singular functions are frequency-independent while for nonself-adjoint systems they are in principle dependent on frequency. However, even for the more general non-selfadjoint problem, the frequency dependence is often not pronounced (cf. Gay and Ray, 1995) so that the steadystate left and right singular functions, i.e. +w (z, 0), G & +v (z, 0),, provide a good basis for expanding the G output and input, respectively, even when y and u are time varying. Thus, for the case of spatial inputs, the steady-state (s"0) form of Eq. (18):



  (20) y(z, 0)" w (z, 0)p (0) v (f,0)u(f, 0) df G G G  G provided the basis for computing the steady-state singular functions from input/output experiments. 2.4. Spatial versus boundary input identification problems When considering the problem of identification of singular functions using either boundary or spatial inputs, it is important to realize that a distributed parameter system intrinsically has both spatial and temporal characteristics. For the case of spatial inputs (Gay and Ray, 1995), steady-state (s"0) information often suffices. By adjusting the spatial profiles of the inputs (e.g. zones of control action or number of discrete actuators along the spatial domain), one can identify the singular functions or good approximations from steady-state experiments. Dynamic information can also be useful for spatial control (e.g. Gay, 1989; Gay and Ray, 1988). However for the case of boundary inputs, only dynamic information is useful for the identification of singular functions.

3. Identification of singular functions In this section we shall describe how one may identify the steady-state singular functions from both spatial inputs

1186

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

and boundary inputs. Also the use of the Karhunen— Loe`ve basis functions in the computations will be discussed. 3.1. Spatial inputs Eq. (6) may be represented by the following operator notation: y"K u, RX

(21)

where the subscript tz is introduced to explicitly represent the spatiotemporal nature of the operator. Similarly, Eq. (7) can also be expressed in operator notation as y(z, s)"K u(z, s), QX

(22)

where the subscript sz in K denotes the Laplace transQX form of K . RX The objective is to obtain the singular functions of the spatial operator alone. In other words, we seek to isolate the spatial component, say K , of the spatiotemporal X operator K or K and then compute its singular RX QX functions. By projecting onto a suitable set of N orthonormal basis functions, +s (z),, using the Galerkin criterion it is G possible to reduce the spatiotemporal DPS represented by Eq. (7) to the following temporal lumped system: b(s)"K(s)a(s),

(23)

where vectors a and b represent the inputs and outputs, respectively, and are expressed as follows:





a (s)  a (s) a(s)"  $

b (s)  b (s) and b(s)"  $

a (s) ,

b (s) ,

with

 

a (s)" G b (s)" G





u(z, s)s (z) dz G

 y(z, s)s (z) dz G 

The N;N matrix K corresponds to the finite dimensional projection of the spatiotemporal operator K onto the +s (z), basis. Issues related to the choice of QX G basis functions are discussed in the subsequent section. For purposes of discussion here, it suffices to say that the basis functions s (z) satisfy the following orthonormality G condition:



 s (z)s (z) dz"d . G H GH 

The elements of the matrix K are



 s (z)K s (z) dz G QX H  It follows that for a self-adjoint DPS, K will be a symmetric matrix. Correspondingly, for a non-self-adjoint DPS, K will be unsymmetric. Steady-state identification corresponds to s"0 in Eq. (23). Thus, K(0) is essentially the steady-state gain matrix for the reduced system. ¹he question one might ask is: ¼hat information does K(0) provide with respect to the distributed parameter system? To answer this question, we consider the example of a parabolic DPS, with a spatially distributed input, represented by the following equation: K " GH

*y "L y#u(z, t). X *t

(24a)

Homogeneous initial conditions: y(z, 0)"0

(24b)

and homogeneous boundary conditions (described by operators B and B ):   B y(z, t)"0, B y(z, t)"0 (24c)   are assumed without loss of generality. The equation is defined on the dimensionless spatial interval 0)z)1 for time t*0. The subscripts 0 and 1 thus refer to z"0 and z"1, respectively. Taking the Laplace transform of Eq. (24a)) yields: L y(z, s)"u(z, s), (25) QX where L "(s!L ). Notice that the differential operQX X ator L in Eq. (25) is the inverse of the integral operator QX K in Eq. (22). Both operators are spatiotemporal. At QX steady-state L "!L indicating that the spatioQX X temporal operator reduces to the underlying spatial operator. This reasoning also suggests that K(0) is essentially the finite-dimensional projection of the spatial operator K . Consequently, the singular vectors of K(0) when X projected back using the set of spatial basis functions +s (z), yield the steady-state singular functions of K (or G QX K ), i.e. the singular functions of K . RX X ¹he next question is: How does one use dynamic information to determine the singular functions of the spatial operator K or the steady-state singular functions of the X spatiotemporal operator K ? QX Eq. (23) can be inverted into the time domain by use of the convolution theorem:



R (26) b(t)" K(t!q)a(q) dq.  Thus, the finite-dimensional projection K(t), of the spatiotemporal operator K onto the +s (z), basis, is the RX G Green’s function of the reduced lumped system.

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

Expressing y(z,t) and u(z,t) in terms of the orthogonal basis functions +s (z),: G , y(z, t)" b (t)s (z), (27a) G G G , u(z, t)" a (t)s (z), (27b) G G G and substituting in Eq. (24a) yields: , , , db s (z) G" b (t) L s (z)# a (t)s (z). (28) G X G G G G dt G G G Applying the Galerkin condition to determine the b and G a , which is essentially taking the inner product of both G sides with s (z), and using the orthonormality condition H of the basis functions +s (z),, one obtains: G db , H" ¸ b #a , (29) HG G H dt G where ¸ " s (z) L s (z) dz or in more succinct vecHG X G  H tor—matrix notation as db "Lb#a, dt

(30)

where the N;N matrix L corresponds to the projection of the spatial operator L . The solution can be written as X R b(t)"eLtb # eLR\Oa(q) dq, (31)   where the subscript 0 denotes the initial condition. For the case of homogeneous (i.e. zero) initial conditions, as in this situation, the solution reduces to





b(t)"

R

eLR\Oa(q) dq.

(32)  Eq. (26) and (32) are input-output models for the same reduced lumped system. The former arises from the integral equation model for the DPS [Eq. (6)] while the latter from the differential equation model [Eq. (24)]. Hence, one obtains: eLt"K(t)

(33a)

or ln K(t) L" , t

(33b)

which directly relates the projection of the spatiotemporal operator K to the projection of the spatial operRX ator L . Note that L "K\, which in the projected X X X space implies that the steady-state value of K(t), K satisQQ fies the following relation: K "lim K(t)"L\, QQ R

(34)

1187

which can be verified by taking the Laplace transform of Eq. (33a) and setting s"0. Some of the terms in Eq. (33) involve functions of matrices and this subject is briefly elaborated on in Appendix B. Thus, in principle one can use dynamic information to obtain the singular functions of K which constitutes the X spatial component of the spatiotemporal operator K . RX There are a number of multivariable identification techniques that can be used to help determine K(t) in Eq. (26) or equivalently K(s) in Eq. (23). Here we use the following procedure. Let us use a form of impulse for the input: u(z, t)"u (z)d(t), (35)  where u (z) is some spatial input function. For this input,  the integral equation (6) reduces to





k(z, f, t)u (f) df (36)   while the PDE representation of the DPS as described by Eq. (24), a non-homogeneous differential equation with homogeneous initial conditions, can be equivalently expressed in terms of a homogeneous differential equation with non-homogeneous initial conditions. The boundary conditions remain the same: y(z, t)"

*y "L y, X *t

(37a)

y(z, 0)"u (z), (37b)  B y(z, t)"0, B y(z, t)"0. (37c)   In terms of the slab heating problem, this idea corresponds to heating the slab momentarily, at t"0 to achieve y(z, 0)"u (z), and then observing the dynamics  at different spatial locations, i.e. the spatiotemporal evolution of the system. For an input of the form described by Eq. (35), the projection space model represented by Eq. (26) reduces to b(t)"K(t)b , (38)  where the elements of the vector b are obtained by  projecting u (z) onto the +s (z), basis. Thus, the meas G ured output for this experiment approximates the Green’s function of the reduced system scaled by the initial condition. The solution to the corresponding PDE, as represented by Eq. (31), reduces to b(t)"eLtb . (39)  In practice, the output measurements, y(z, t) and consequently b(t), are not continuous but discrete in time. Hence, practical implementation necessitates obtaining the following discrete-time analog to Eq. (39): b(*t)"eL Rb ,  b((k#1)*t)"eL Rb(k*t), k"1,2,2,M.

(40a) (40b)

1188

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

The sampling intervals *t (assumed uniform) are chosen to be small enough to approximate the continuous behavior as closely as possible. For convenience, we will use the subscript k instead of k*t to denote the kth time instant, i.e. b "b(k*t). In this discrete time framework, I the quantity e* R can be interpreted as a map, P, which transforms b to b . Thus, we have I I> ln P eL R"P or L" . (41) *t The problem at hand is to get an estimate of P from knowledge of b for k"1, 2,2, M, assuming M samples I taken at intervals of *t. If matrices X and X are defined   such that their columns correspond to b and b reI> I spectively for k"1,2,2,M, then X "P X . (42)   Since the measurements are taken until steady state is achieved, b "b . Hence, the Mth column of X can + +>  be set equal to the (M!1)th column. Both X and X are   of dimension M;N and thus non-square. Hence, Eq. (42) cannot be solved exactly for P. However, it is possible to obtain a pseudo-inverse P : P "(X X2)(X X2)\ (43)     which is a good solution in the sense of Frobenius norm of the residual matrix, [X !PX ] (Golub and van   Loan, 1989). Thus, one can obtain P , which is a ‘good’ approximation of the map P. L can subsequently be calculated through use of Eq. (41). These results can be used to determine either the singular functions or the eigenfunctions (when they are different) as one wishes. Determination of steady-state singular functions: Singular Value Decomposition of the matrix L yields the coefficients, in the +s (z), basis, of the singular functions G of the operator L , or equivalently K\ . The singular X X vectors of L when projected back using the +s (z), basis G yield the singular functions of L . Since L "K\, the X X X left and right singular functions of L are the right and X left singular functions of K . The reciprocals of the X singular values of L correspond to the singular values X p of K . G X Determination of eigenfunctions and adjoint eigenfunctions: The methodology simplifies if one seeks to determine the eigenfunctions and adjoint eigenfunctions of L rather than the singular functions. The operation of X taking the logarithm of matrix P as in Eq. (41) becomes redundant since the eigenvectors of P and L are identical. The eigenvalues are related quite simply as: ln(jP ) jL" . *t

(44)

The column eigenvectors of P when projected back using the +s (z), basis yield the eigenfunctions of L while the G X row eigenvectors yield the adjoint eigenfunctions.

3.2. Some choices for basis functions In carrying out the computations in the last section, one must choose a suitable set of orthogonal basis functions, s (z). In Gay and Ray (1995), B-splines were usually G selected; however any complete set of orthogonal basis functions will suffice. Hence when choosing a basis, the primary questions to be posed are 1. Can the important information be captured using as few basis functions as possible? 2. Can the basis functions be chosen keeping in mind the physical aspects of each particular DPS, yet be universal in the application sense? These considerations suggest the use of the K—L expansion, a technique used by dynamicists in the analysis of complex spatiotemporal patterns. In the current context, the spatiotemporal pattern corresponds to the output of the DPS, namely y(z, t). A useful feature of the K-L decomposition is that the appropriate set of orthonormal basis functions can be computed from the experimental data itself and these basis functions are optimal for representing that data set. From the output y(z, t), a symmetric spatial correlation operator S with kernel S(z, z) is constructed as follows:



1 2 S(z, z)" y(z, t) y(z, t) dt. (45) ¹  The eigenfunctions of the symmetric spatial correlation operator S constitute the desired set of orthonormal spatial basis functions t (z), so that St "j t or in G G G G more elaborate notation:









S(z, z)t (z) dz"j t (z). (46) G G G  Consequently, +t (z), are referred to as spatially depenG dent eigenfunctions or topos. Orthonormality implies:

t (z)t (z) dz"d . G H GH  The eigenvalues j refer to the energies of the respective G modes. The modes are arranged in descending order of their energies: j *j *2*0.   The corresponding time varying amplitudes, d (t), which G can be obtained by projecting the output, y(z, t), onto the topos:



d (t)" G



y(z, t)t (z) dz (47) G  are essentially the eigenfunctions of the symmetric temporal correlation operator T with kernel

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

q(t, t)"y(z, t)y(z, t) dz and hence satisfy the following  relation analogous to the one for topos:



1 2 q(t, t) d (t) dt"j d (t). (48) G G G ¹  The temporal correlation operator T will have the same set of eigenvalues j as the symmetric correlation operG ator S. The amplitudes are also referred to as timedependent eigenfunctions or chronos and satisfy the following orthogonality condition:



2




2



d (t) dt d (49) G GH   making the K-L expansion time-space symmetric (Aubry et al., 1991). The output or spatiotemporal pattern, y(z, t), can generally be expressed in terms of many sets of basis functions that provide convergence when used in an infinite series expansion: d (t) d (t) dt" G H

 y(z, t)" b (t)s (z). (50) G G G However, from a practical standpoint, one would like to have a truncated series expansion with as few terms as possible, say N. For such an N-term expansion, , y (z, t)" b (t)s (z) (51) , G G G the K-L expansion provides the optimal basis in the least squares sense for the data set being analysed (Butkovskiy, 1969; Fukunaga, 1990), i.e. the residual R defined as , 1 2  (y(z, t)!y (z, t)) dz dt (52) R " , , ¹   is minimized for the choice: b (t)"d (t) and s (z)"t (z). G G G G It should be stressed that the basis functions obtained using the K-L procedure are optimal only for the data set analysed and thus only for the particular input used. The K-L expansion decomposes the output (in response to a particular input) orthogonally with respect to space and time as follows:








, , y(z, t)" d (t)t (z) " l c (t)t (z) , (53) G G G G G G G where c (t) refers to the normalized version of the chronos G and l ("(j ) provides a measure of the energy conG G tained in each mode. The introduction of these two new quantities is motivated by the use of a computational framework based on the matrix SVD for numerical implementation of the K-L decomposition (see Appendix C). The identified topos will depend on the choice of input in addition to the intrinsic nature of the DPS. On the other hand, the singular functions (and eigenfunctions)

1189

are properties of the DPS and independent of the input. Thus the K-L expansion is most useful as a computational tool in calculating the eigenfunctions or singular functions of the DPS. Chen and Chang (1992) also developed a framework for the identification of eigenfunctions and adjoint eigenfunctions of an unknown DPS using the K-L decomposition. Their contribution differs from ours in that, (i) the emphasis here is on identification of singular functions as opposed to eigenfunctions and adjoint eigenfunctions as in their work; and (ii) the identification scheme developed here provides reasonable results with only a small number of sensors, whereas their work required a very large number of sensors. 3.3. Boundary inputs The classical eigenfunction—eigenvalue problem is characterized by a non-homogeneous differential equation with homogeneous boundary conditions. For the DPS described by Eq. (24), the associated eigenfunction—eigenvalue problem is L (z)"j (z), n"1,2,3,2 , (54a) X L L L where j are the eigenvalues and (z) are the corresL L ponding eigenfunctions. The boundary conditions are B y"0, (54b)  B y"0. (54c)  In contrast, the boundary control problem is posed in terms of a homogeneous differential equation: *y(z, t) "L y(z, t) X *t

(55a)

with non-homogeneous boundary conditions: B y"u(t) and B y"0 (55b)   Again, the equation is defined on the dimensionless spatial interval 0)z)1 for time t*0. The subscripts 0 and 1 refer to z"0 and z"1, respectively. In order to determine the boundary Green’s function, we will use u(t)"1 to provide the data for analysis. Thus, in order to develop a methodology for obtaining the eigenfunctions or more generally, the singular functions, it is necessary to decompose the boundary input problem into two constituent problems: 1. ¹he steady-state problem: L y(z)"0; B yN "1, B yN "0 (56) X   which captures the inhomogeneity in the boundary conditions. The overbar is used to denote steady-state. 2. ¹he transient problem which contains the pertinent dynamic information for determination of the

1190

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

eigenfunctions or the singular functions. Letting y(z, t)"y(z)#½(z, t), yields: *½(z, t) "L ½(z, t); B ½"0, X  *t

B ½"0 

(57)

with initial conditions ½(z,0)"!y(z). The equivalent integral equation representation for the DPS of Eq. (57) is



½(z, t)"



k(z, f, t)[!y(z)] dz,

(58)

 where k(z, f, t) is the Green’s function of the DPS. Notice that Eq. (58) is identical in form to Eq. (36). Consequently, the methodology for identification of the singular functions of L or K exactly follows the approach X X described in Section 3.1 using spatial inputs. We provide here a brief summary of the main steps. 1. The dynamic profiles, y(z,t), and the steady-state profile y(z), in response to a unit step input at the boundary are obtained. 2. ½(z, t) is obtained by subtracting the steady-state contribution from the dynamics, ½(z, t)"y(z, t)!y(z). 3. K-L decomposition is performed on ½(z, t) to obtain a set of topos. These form the basis functions for projection onto a coefficient space (chronos). 4. From the chronos, a good approximation, P , to the map P is constructed [see Eq. (40), (42) and (43)]. The projection, L, of the spatial operator L onto the X topos basis is computed through use of Eq. (41). 5. The singular vectors of L when projected back using the topos yield the desired singular functions. Similarly an expansion of the eigenvectors of L provides the eigenfunctions of the DPS.

4. Examples of identification with boundary inputs The distributed parameter systems chosen as examples belong to the class of second-order linear parabolic partial differential equations with constant coefficients as shown in Eq. (1). Neumann boundary conditions are employed, so that *y (0, t)"u(t), *z

(59a)

*y (1, t)"0. *z

(59b)

The equation is defined on the dimensionless spatial interval 0)z)1 for time t*0. The subscripts 0 and 1 thus refer to z"0 and z"1, respectively. In simple terms if the DPS represents a metal plate or a tubular reactor, we seek to identify the singular functions of the DPS from the dynamic information obtained by adjust-

Fig. 1. Heating of a metal slab from one end.

ing the heat flux at one end (z"0) of the metal plate (see Fig. 1) or the temperature/concentration of the feedstream to the reactor. There are also other factors which warrant consideration in this identification methodology: 1. Sampling time, duration and number of samples. When implementing the identification procedure in real time, samples are taken over the spatial domain at uniform time intervals. Obviously, the sampling time can alter the identification results significantly. Since the primary emphasis is on understanding the spatial structure of the problem, the approach adopted was to make the sampling time small enough so as to obtain a ‘good’ approximation of the continuous-time behavior. The general rule-of-thumb with respect to sampling duration was to ensure that both short-time and long-time behavior was included in the data set. Thus, in these examples, our procedure was to take 100 samples and a ‘small enough’ sampling time of 0.5 s. 2. Number of sensors for output measurement. Generally, it was assumed that a relatively small number of sensors (five or six) were available for output measurement. To provide a benchmark for assessment of the identification results with a few sensors, simulations are also performed with a large number of sensors (about 21). For these examples, the sensors are assumed to be equi-spaced. 3. Comparison with spatial identification results. All results of the dynamic boundary identification scheme are compared to the steady-state identification results of Gay (1989). For self-adjoint problems, the singular functions are independent of time and are the eigenfunctions which are known analytically. For non-selfadjoint systems, left and right steady-state singular functions are determined to a high degree of precision by employing 50 zones of spatial input and 51 sensors for output measurement. These singular functions will be referred to as ‘actual’ singular functions. Henceforth, it is also to be assumed that even when not explicitly mentioned, ‘singular functions’ actually refers to the steady-state singular functions. 4. Dimension of reduced temporal system. A key step in the dynamic identification framework is projecting the spatiotemporal DPS onto a finite set of orthonormal basis functions to obtain a reduced temporal system of finite order, say N. The choice of N is crucial. For

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

example, a very small value of N, would neglect important information. Choosing an arbitrarily large value of N has disadvantages from an implementation standpoint as the problem can often become ill-conditioned. In such a situation, the condition number, i provides a stable quantifiable measure of the conditioning of the reduced order system. Since the topos, obtained from a K-L decomposition of the output, are used as the spatial basis functions, the conditioning of the reduced order system can be assessed by inspecting the ratio of the K-L singular values, i"l /l , in  , Eq. (53). The examples below will show that in general a value of i around 100 provides a reasonable compromise. 4.1. Example 1: Self-adjoint system with Neumann boundary conditions The first example corresponds to a self-adjoint system defined by Eq. (1) with Nu"20, Pe"0 and Neumann boundary conditions [Eq. (59)]. Physically, this corresponds to heating a metal strip at one end. The left and right singular functions are identical to the eigenfunctions for this model problem which are known to be +1, (2 cos(i!1)nz for i*2,. Fig. 2 shows that the boundary identification scheme yields nearly exact results when using 21 sensors for output measurement and the first 6 topos to construct the reduced system as indicated by the solid lines. Notice that only the first three pairs of singular functions are presented. For higher singular functions ('3) the identification results were not so precise. The singular values are nearly identical to the analytically determined values (Table 1). The more stringent test case of five sensors at z"0.0, 0.0.25, 0.5, 0.75, 1.0 and using four topos (l /l "115.64)   to construct the reduced system also yields comparable results. As illustrated in Fig. 2, the left singular functions, which constitute an orthonormal basis for the output space, are in slightly better agreement as compared to the

1191

right singular functions, which constitute an orthonormal basis for the input space. However, both the identified singular functions are reasonable approximations of the eigenfunctions. The singular values also are in quite reasonable agreement with exact values. With five sensors a maximum of five topos can be identified through a K-L decomposition of the dynamic output profiles translated by the steady-state profile. However, only three pairs of singular functions may be identifiable with a reasonable level of accuracy in this example. Since l /l "2090.61, using five topos would   yield an extremely ill-conditioned reduced order system. However, the fact that l /l "32.8 would suggest that   using three topos would neglect valuable information. Fig. 3 clearly depicts the deteriorated identification results with the use of three topos. This observation is independent of the number of sensors used, as the results with 21 and 5 sensors demonstrate (see also Table 1). 4.2. Example 2: Non-self-adjoint system with Neumann boundary conditions This example corresponds to a non-self-adjoint system defined by Eq. (1) with Nu"20, Pe"10 and Neumann boundary conditions [Eq. (59)]. The Peclet number is essentially a measure of the convective contribution or in

Table 1 Singular values of DPS of Exmaple 1 using boundary inputs. Comparison of exact singular values with the singular values obtained by means of the boundary identification scheme using 21 sensors, 6 and 3 topos; 5 sensors, 4 and 3 topos p G

Actual

Boundary identification 21 Sensors

p  p  p 

0.05000 0.03348 0.01681

5 Sensors

6 Topos

3 Topos

4 Topos

3 Topos

0.04990 0.03348 0.01681

0.04874 0.02966 0.01083

0.05161 0.03184 0.01810

0.05289 0.02916 0.01065

Fig. 2. Identified right and left singular functions [l (z) and u (z)] for G G DPS of Example 1 using boundary inputs. Solid and dashed lines correspond to the identified singular functions using 21 sensors, 6 topos and 5 sensors, 4 topos, respectively.

1192

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

Fig. 4. Empirically determined eigenfunctions, (z), for Example 2 usG ing boundary inputs.

Fig. 3. Identified right and left singular functions (l (z) nad w (z)) for G G DPS of Exmaple 1 using boundary inputs and 3 topos. Solid lines correspond to actual eigenfunctions, long and short dashes to indentified singular functions using 21 sensors, 3 topos and 5 sensors, 3 topos, respectively.

mathematical terms, a measure of the non-self-adjointness. A complete set of eigenfunctions, though not orthogonal, does exist for this system [7]. The current identification framework can be used to compute the eigenfunctions and eigenvalues for such a non-self-adjoint system by projecting the eigenvectors of L onto the +s (z), basis, as indicated in the discussions at the end of G Sections 3.1 and 3.3. The identified eigenvalues of the PDE system (1) are j "20.0178, j "54.5140 and   j "88.3046 which compare favorably with the actual  values of 20.0, 54.8696 and 84.4784, respectively. Fig. 4 depicts the non-orthogonal nature of the identified eigenfunctions. Use of such a basis leads to slow and often non-monotonic convergence to desired output profiles (Gay and Ray, 1995). A purpose of this example was to illustrate the usefulness of singular functions as compared to eigenfunctions in such situations. Fig. 5 represents the identified left and right singular functions using 21 sensors and 6 topos to construct the reduced order temporal system. The solid lines correspond to the ‘actual’ singular functions computed using the steady-state SVD-based method of Gay and Ray (1995). The dashed lines represent the results of the present work using boundary inputs. It is possible to identify the left singular functions, which correspond to the out-

Fig. 5. Identified right and left singular functions (l (z) and w (z)) for G G DPS of Example 2 using boundary inputs and 21 sensors, 6 topos. Solid lines correspond to ‘actual’ singular functions using spatial inputs (Gay and Ray, 1995) and dashed lines to the identified singular functions using boundary inputs and 21 sensors, 6 topos.

put space, and the singular values (Table 2) almost exactly. However, the identified right singular functions (which correspond to the input space), while a reasonable approximation for the ‘actual’ right singular functions, display some oscillations. Even with only six sensors for output measurement, it is possible to identify a reasonable set of singular functions for this non-self-adjoint system using boundary

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204 Table 2 Singular values of DPS of Exmaple 2 using boundary inputs. Comparison of ‘actual’ singular values obtained using spatial inputs (Gay and Ray, 1995) with singular values from boundary inputs using 21 sensors, 6 topos and 6 sensors, 3 topos p G

p  p  p 

Spatial identification (actual)

Boundary identification 21 Sensors, 6 topos

6 Sensors, 3 topos

0.07542 0.02265 0.01152

0.07420 0.02259 0.01149

0.07809 0.02273 0.01218

Fig. 6. Identified right and left singular functions [l (z) and u (z)] for G G DPS of Example 2 using boundary inputs and 6 sensors, 3 topos. Solid lines correspond to ‘actual’ singular functions and dashed lines to identified singular functions with boundary inputs using 6 sensors, 3 topos.

inputs as shown in Fig. 6. Only three topos are used to construct the reduced-order temporal system (l /l "40.55, l /l "196.70). Table 2 provides a direct     comparison of the singular values obtained with six sensors and three topos to the spatial identification results of Gay and Ray (1995) as well as the boundary identification results using 21 sensors and 6 topos. The agreement with exact values in all cases is quite good.

5. Boundary topos Eq. (4) provides a general representation for linear DPS with boundary inputs. The unit step response S(z, t),

1193

and consequently the boundary Green’s function g(z, t) [see Eq. (5)], can be ascertained experimentally through the use of a unit step input at the boundary. Application of the K-L decomposition (Section 3.2) to the boundary Green’s function, g(z,t), or equivalently the step response function, S(z, t)"R g(z, t) dt, would yield a set of spatial  basis functions (topos) that are an optimal basis in the least squares sense for representation of dynamic output profiles that arise in response to boundary step inputs. These functions will be referred to as boundary topos in this work. We use the adjective boundary to explicitly refer to the form of input used. The boundary topos, although purely empirical in nature, can be identified systematically for unknown DPS and subsequently put to effective use in feedback control applications. Since for linear autonomous systems all possible spatial profiles are time integrations of the boundary Green’s function Eq. (4), one could argue that those boundary topos are ‘optimal’ spatial basis functions for boundary inputs. When using K-L decomposition for analysis of spatiotemporal patterns or complex dynamics arising in nonlinear DPS, it is standard practice to first subtract the time-averaged part of the data and apply the K-L decomposition to the time-varying part. This is because the K-L modes are sought to explain the dynamics that constitute the primary fluctuations from some mean asymptotic behavior. In the current context, we are dealing with a linear DPS with a boundary input and hence the issue of complex dynamics is not pertinent. When developing the framework for identification of singular functions using boundary inputs, the first step in Section 3.3 was to transform the DPS characterized by nonhomogeneous boundary conditions to one with homogeneous boundary conditions. Notice that the standard practice of subtracting the time-averaged part for nonlinear DPS, when applied to linear DPS with boundary inputs indirectly achieves the same effect, i.e. it transforms the problem to one with homogeneous boundary conditions. For linear DPS, subtracting the time-averaged profile is tantamount to a simple translation of the spatial profiles by redefining the zero profile to correspond to the time-averaged profiles. This would suggest that one could also obtain an equivalent set of boundary topos by subtracting the final steady-state profiles from the transient profiles. Thus in order to obtain a set of boundary topos, the K-L decomposition may be applied to any one of the following: 1. The actual step response, S(z, t). 2. The step response translated by the time-averaged profile, S(z, t)!S(z), where S(z)"1/¹2S(z, t) dt,  ¹ being the sampling duration. 3. The actual step response translated by the final steady-state profile, S(z, t)!S(z), where S(z)" S(z, t P R).

1194

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

In the procedure for determining the singular functions of an unknown DPS with boundary inputs, the spatial basis functions (or topos) obtained from K-L decomposition on ½(z, t)"y(z, t)!y(z) were used to project the spatiotemporal DPS and obtain a reduced temporal system (see Section 3.3 for details). Closer perusal of the list above would indicate that performing the K-L decomposition on +S(z, t)!S(z), essentially yields the same set of basis functions, which we refer to as boundary topos. (For the sake of clarity, we will employ the notation u (z), G to denote the boundary topos as opposed to t (z) for the G topos. We will continue to use the symbol l to denote the G corresponding K-L singular value.) Computing the boundary topos is easier than computing the singular functions since there is no need to go through the steps of determining spatial basis functions and using the condition number to obtain a well-conditioned reduced temporal system. Consequently, the issues related to ill-posedness of problems becomes irrelevant. Presented now are the boundary topos for the two examples studied in Section 4. 5.1. Example 1: self-adjoint DPS with Neumann boundary conditions Fig. 7 shows the first three identified boundary topos, u (z) (obtained from a K-L decomposition of the output G

Table 3 K-L singular values of the DPS of Example 1 using boundary inputs. Comparison of results with 21 and 5 sensors l G

21 Sensors

5 Sensors

l  l  l  l  l 

1.3489;10\ 1.7627;10\ 4.3660;10\ 1.0969;10\ 2.2422;10\

1.3244;10\ 1.6041;10\ 4.0397;10\ 1.1453;10\ 6.3351;10\

translated by the steady-state profile in response to a step input at the boundary). The results with five sensors compare well to the results with 21 sensors. Table 3 compares the K-L singular values. Notice that using five sensors does not yield an accurate estimate of l .  However, this is inconsequential because one rarely uses the fifth mode since l is about three orders of magni tude less than l . Also, use of a step boundary input at  z"1 instead of z"0 leads to the identification of boundary topos u* (z), which are exact mirror images of G u (z) as illustrated in Fig. 7. This, of course, is due to the G self-adjoint nature of the DPS. 5.2. Example 2: Non-self-adjoint DPS with Neumann boundary conditions The identification results with 6 and 21 sensors for the non-self-adjoint system of Example 2 are presented in Fig. 8 and Table 4. With 6 sensors the first two modes are identified almost exactly. The identification results deteriorate for the higher modes although one still obtains reasonable approximations for the third, fourth and perhaps even the fifth mode.

6. Pseudo-modal feedback control scheme

Fig. 7. Identified boundary topos reflecting the symmetric nature of the DPS of Example 1. u (z) corresponds to the usual case of the input at G z"0 while u* (z) corresponds to the case of the input at z"1. Solid G lines and dashed lines correspond to identified boundary topos using 21 and 5 sensors, respectively.

The discussion in Sections 3 and 5 would suggest that there are two candidate sets of basis functions — steadystate left singular functions, +w (z, 0),, and boundary G topos, u (z), that can be employed in a pseudo-modal G feedback control scheme for representation of the output profiles. It will be demonstrated that reasonable results can be obtained using either the left singular functions or the boundary topos. A schematic diagram of the pseudo-modal feedback control scheme is provided in Fig. 9. For distributed parameter systems, the setpoint corresponds to a spatial profile as opposed to a particular value for lumped parameter systems. To keep the discussion general we will allow for the setpoint to be either a static profile or a dynamically varying profile. Therefore, the setpoint is represented as y (z, t), which reduces to y (z) for the QN QN case of the setpoint being a static profile. The objective in

1195

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

the control problems studied below is either achieving a particular setpoint from a homogeneous state or maintaining an existing profile in the presence of disturbances. e(z, t) is the deviation of the output profile, y(z, t), from

the setpoint, y (z, t), and is computed as QN e(z, t)"y (z, t)!y(z, t). (60) QN The pseudo-modal error vector is then computed by projecting the deviation profile onto the boundary topos or alternatively, the left singular functions. For purposes of discussion here, we use the boundary topos. Thus,



e (t)" G







(y (z, t)!y(z, t))u (z) dz" QN G

 

e(z, t)u (z) dz. G

(61) The modal controller computes the control coefficients a (t) corresponding to each of the modes. Any suitable G form of control law can be used. For example, with proportional control, a (t)"Ke (t). (62) G G The net control action, which constitutes the boundary input, can be computed as follows: , 1 u(t)" a (t), (63) l G G G where use of the inverse of the K-L singular value can be justified as follows. The simplified integral equation representation for the DPS with boundary inputs is



y(z, t)" Fig. 8. Identified boundary topos for non-self-adjoint DPS of example 2 with Neumann boundary conditions. Solid lines and dashed lines correspond to identified boundary topos using 21 and 6 sensors, respectively.

Table 4 K-L singular values of the DPS of example 2 using boundary inputs. Comparison of results with 21 and 6 sensors l G

21 Sensors

6 Sensors

l  l  l  l  l  l 

8.1689;10\ 1.0034;10\ 2.0813;10\ 4.8340;10\ 1.0537;10\ 1.8538;10\

8.1731;10\ 1.0029;10\ 2.0157;10\ 4.1553;10\ 8.6578;10\ 1.4839;10\



g(z, t!q)u(q) dq. (64)  Expanding the output, y(z, t), in terms of the boundary topos, u (z), and using the optimal spatiotemporal exG pansion for the boundary Green’s function (see Section 5) yields:



 b (t)u (z)" u (z)l c (t!q)u(q) dq. (65) G G G G G  G G Taking inner products of both sides of the above equation with u (z) and using the orthonormal property of the H boundary topos, one obtains,



b (t)"l H H



c (t!q)u(q) dq (66) H  which would suggest use of 1/l as a scaling factor in Eq. G (63).

Fig. 9. Boundary pseudo-modal feedback control scheme.

1196

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

Eq. (61) would seem to indicate that in order to compute the modal error vector e(t), one needs to reconstruct the entire spatial profile at each time instant from the available measurements and then perform the integration. However in practice this is accomplished quite easily through the use of weighting matrices. Assuming N measurement locations, it is possible to express both the setpoint profile as well as the actual output in terms of a set of N linearly independent trial basis functions, +q (z),, as follows: I , y(z, t)" q (z)b (t), (67a) I I I , y (z, t)" q (z)b (t). (67b) QN I IQN I Corresponding to the set of N trial basis functions are a set of N linear functionals, +q*( ) ),, which are used to I determine the trial basis coefficients b (t) and b (t) based I IQN on the output measurements: b (t)"q* (y(z, t)), (68a) I I b (t)"q* (y (z,t)). (68b) I QN IQN The deviation profile, e(z, t), can also be expressed in terms of the trial basis functions as follows: , e(z, t)" q (z)[b (t)!b (t)]. I IQN I I Similarly for the boundary topos we have

(69)

, u (z)" q (z)u , (70) G H HG I where u is the coefficient of the jth trial basis function HG q (z) in the series approximation to u (z). Substituting Eq. H G (69) and (70) in the expression for the modal error vector shown in Eq. (61), one obtains:



 , , q (z)(b (t)!b (t)) q (z)u dz I IQN I H HG  I H which can be simplified to e (t)" G




(71)

of output trial basis coefficients, b(t), compare it to the vector of setpoint trial basis coefficients, b (t), and then QN scale the difference by the quantity U2Q. This is completely equivalent to performing the integration in Eq. (61). The performance of the pseudo-modal control scheme is now studied through a few examples.

7. Examples of boundary control Three example applications are considered. The first two correspond to linear distributed parameter systems. The last example is of a nonlinear tubular chemical reactor. In all cases, only a finite number of sensors (five to six) are assumed to be available for measurement of the output. 7.1. Example 1: pure diffusion/conduction problem with Dirichlet boundary conditions The first example corresponds to a self-adjoint system with Nu"0, Pe"0 and Dirichlet boundary conditions. Mathematically, it can be represented as follows: *y *y "a , *t *z

(74a)

y(0, t)"u(t),

(74b)

y(1, t)"0.

(74c)

Homogeneous initial conditions are assumed without loss of generality. Five pointwise output measurements at z"0.0, 0.25, 0.50, 0.75, 1.0 are considered. Physically, this corresponds to controlling the temperature profile in a stationary metal slab (Pe"0) by varying the temperature, u(t), at one end, z"0. Also, no external heat transfer is assumed (Nu"0). Using the approach outlined in Section 5, a set of boundary topos is identified as shown in Fig. 10. The



 , , e (t)" u2 q (z) q dz (b (t)!b (t)). (72) H I IQN I G GH  H I In vector—matrix notation, the above equation can be expressed as e(t)"U2Q(b (t)!b(t)), (73) QN where U is the matrix of trial basis coefficients of the boundary topos and Q is the Gram matrix constructed from the known trial basis functions. The inner products,  q (z)q (z) dz, can be exactly computed through the use H  I of quadrature rules (Schumaker, 1981). In this work, the choice of cubic B-splines (order 4) proved to be an adequate trial basis. Thus, when actually performing the control calculations, the practice is to compute the vector

Fig. 10. Identified boundary topos for the DPS of Example 1 with pure diffusion and Dirichlet boundary conditions.

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

corresponding K-L singular values are l "1.4078,  l "0.2057 and l "0.0269. One could further compute   the actual left and right singular functions, which in the case of this model problem are identical to the eigenfunctions, known to be +(2 sin inz for i*1,. However, for purposes of this example, the boundary topos are chosen as the basis for the pseudo-modal feedback control scheme. Also, only the most dominant mode is used in computing the control action. Thus the summation in Eq. (63) reduces to a single term indicating that the controller essentially tries to minimize e , the error be tween the projections of the output profile and the specified setpoint profile onto the first boundary topos [see Eq. (61)]. Fig. 11 shows how a PI-controller (K"0.5, 1/q "50), ' based on one mode, exhibits some initial overshoot and finally achieves the desired linear setpoint profile. Clearly for such a diffusion process, the only profiles that can be achieved closely at steady-state are linear profiles like the setpoint in Fig. 11. Fig. 12 illustrates the performance of the pseudo-modal controller (K"0.5, 1/q "50) for a different setpoint ' profile, one that cannot be achieved exactly at steadystate. Also shown are the profiles achieved using two lumped parameter single-input single-output (LPS SISO) controllers. For both cases, the input is the temperature at the boundary, u(t), while the output is either the temperature at the mid-point, y(0.5, t), or at z"0.25, i.e. y(0.25, t). Notice that the two LPS controllers (K"1, 1/q "50) minimize the error exactly at z"0.25 or ' z"0.5. On the other hand, the pseudo-modal controller takes into account the entire spatial information and tries to make the deviation profile, e(z), orthonormal to the first boundary topos. In this particular instance, a LPS SISO controller with output at z"0.35, would perform

Fig. 11. Pseudo-modal feedback control for setpoint change from y(z)"0 to y(z)"1!z. 1 mode of control action with K"0.5 and 1/q "50 is used. The setpoint is exactly achieved at t"1.2. Solid line ' corresponds to the setpoint profile. Dashed lines correspond to the transients.

1197

almost exactly like the DPS controller. This will not be the case in general. Since the qualitative nature of the profiles attainable for this pure-diffusion process at steady-state are fairly limited, a study was undertaken to see how the system responds to time-varying inputs. In other words, the performance of the pseudo-modal controller was tested for the case of a time-varying setpoint. The setpoint was obtained as the output of the system described in Eq. (74) for the case when the temperature at the boundary was varied in a sinusoidal fashion, u(t)"1#sin(10t). The desired profiles are no longer just linear as shown by the solid lines in Fig. 13a. For static setpoints, the primary role of integral action is to eliminate steady-state offset. This does not apply in general for time-varying setpoints (depends on the frequency of variation). Improper choice of q could easily ' cause the system to go unstable. Hence only proportional control is used. Snapshots of the setpoint profiles at different instants of time and the controlled profiles using 1 mode of control action with a proportional gain of 10 are shown in Fig. 13a. The time-history at the midpoint is also shown in Fig. 13b. With one mode of control action and a gain, K"100, the pseudo-modal controller is able to achieve the desired periodic setpoint. 7.2. Example 2: Non-self-adjoint problem with Neumann boundary conditions This example corresponds to a non-self-adjoint system with Nu"20, Pe"10 and Neumann boundary conditions: *y *y *y ! #10 #20y(z, t)"0, *t *z *z

(75a)

Fig. 12. Comparison of DPS and LPS SISO controllers for Example 1. 1 mode of control action with K"0.5 and 1/q "50 is used for the  DPS cotnroller. The two LPS SISO controllers (K"1, 1/q "50) use ' as setpoint y (z"0.5) and y (z"0.25). QN QN

1198

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

Fig. 14. Identified left singular functions for the non-self-adjoint DPS of Example 2 with Neumann boundary conditions.

Fig. 13. Boundary control to time-varying setpoint profile with 1 mode and K"100. (a) The spatial profiles at different instants and (b) history at the midpoint. Solid lines correspond to the setpoint and the dashes to the controller output. Fig. 15. Pseudo-modal feedback control for setpoint change from y(z)"0. 1 mode of control action with K"0.5 and 1/q "20 is used. ' The setpoint is achieved at t"0.6. Solid line corresponds to the setpoint profile. Dashed lines correspond to the transients.

*y (0, t)"u(t), *z

(75b)

*y (1, t)"0. *z

(75c)

Homogeneous initial conditions are assumed without loss of generality. Six pointwise output measurements at z"0.0, 0.2, 0.4, 0.6, 0.8, 1.0 are considered. This DPS was studied in Section 4 (see Example 2). Physically it could correspond to an isothermal tubular reactor with a first order reaction where the concentration profile is controlled by means of the inlet flux or alternatively, a moving heated slab with Newton cooling along its length and heating at z"0. The pseudo-modal controller for this example uses the identified left singular functions as shown in Fig.14. The associated singular values are p "0.07809, p "   0.02273 and p "0.01218.  Fig. 15 shows how the pseudo-modal controller, using 1 mode of control action gradually achieves the desired

setpoint profile. The controller gain is set to 0.5 and an integral time constant 1/q "20 is used. ' The setpoint in Fig.16 is a linear profile, one that cannot be exactly achieved using a boundary input. Notice that a LPS SISO controller (K"1, 1/q "50), with ' the output at z"0.5, exactly achieves the setpoint at z"0.5. On the other hand, the DPS controller drives the system to a steady-state profile so that the deviation from the setpoint is orthogonal to the first left singular function. 7.3. Example 3: nonlinear tubular chemical reactor with Danckwerts boundary conditions The last example is of a non-isothermal tubular reactor in which a first-order, irreversible, exothermic reaction takes place. This is a nonlinear distributed parameter system. The source of the nonlinearity is the exponential

1199

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204 Table 5 Dimensionless parameters for tubular reactor model Parameter

Value

Definition

Pe F Pe K ¸e a b d c g U m G g G

5 5 1 0.875 0.5 25 13 1.0 0.85 1.05

Heat Pe´clet number Mass Pe´clet number Lewis number Damko¨hler number Heat of reaction Activation energy Heat transfer coefficient Wall temperature Inlet concentration Inlet temperature (control variable)

Fig. 16. Comparison of DPS and LPS SISO controllers for Example 2. 1 mode of control action with K"0.5 and 1/q "20 is used for the ' DPS controller. The LPS SISO controller (K"1, 1/q "50) uses as ' setpoint y (z"0.5). QN

dependence of the reaction rate on temperature. Hence in the current study we focus only on the temperature behavior of the system. This particular model has been taken from the work of Alvarez et al. (1981). Similar reactor models have been used by Gay (1989) and Chen and Chang (1992) for investigative studies of identification and control of nonlinear DPS. For more details on the model, the reader is referred to the article of Alvarez (1981). Here the coupled mass and energy balances are directly given in dimensionless form: *m *m *m " !Pe !ameB\E, K *z *t *z

(76a)

*g 1 *g *g " !Pe #abmeB\E#c(g !g) F *z U *t ¸e *z

(76b)

as well as the Danckwerts boundary conditions: *m (0, t)"Pe (m!m ), K G *z

(77a)

*g (0, t)"Pe (g!g ), F G *z

(77b)

*m (1, t)"0.0, *z

(77c)

*g (1, t)"0.0. *z

(77d)

The above model of the reactor can be viewed as representing either an empty tubular chemical reactor or a packed bed reactor represented by a pseudo-homogeneous model. Table 5 contains the definitions and values of the dimensionless parameters (Alvarez et al. 1981). m refers to the dimensionless concentration and g to the dimensionless temperature. The reference temperature is 500 K.

Fig. 17. Input/output structure for the tubular reactor.

The input—output structure for this tubular reactor is shown in Fig. 17. It is assumed that there are five sensors located at z"0.0, 0.25, 0.5, 0.75, 1.0 for measurement of output temperature. Concentration is not measured. In his work, Gay (1989) used as input 3 independent zones of wall cooling (spatially distributed input). In the current work, a single zone of fixed cooling is assumed. The control variable is the feed temperature (boundary input). Qualitatively speaking, the two control variables will achieve different effects. Windes (1986) whose dissertation constituted an extensive study on the modelling and control of packed bed reactors, noted that the wall temperature influences the temperature throughout the reactor and is a good choice when trying to control the maximum reactor temperature or the overall yield at the reactor exit. Wall temperature has little effect on the location of the hot spot. On the other hand, inlet temperature and flow rate can significantly affect the location of the hot spot. The identification and control results presented are based on temperature measurements alone using the full nonlinear model. Before proceeding to identify the characteristic boundary topos, the steady-state version of Eqs. (76) and (77) is solved to obtain the base temperature and concentration profiles. This requires solution of a boundary value problem. The differential equations were discretized using the method of orthogonal collocation on finite elements (Villadsen and Stewart, 1967; Villadsen and Michelsen, 1978; Carey and Finlayson, 1975). Forty elements with 3 collocation points were found to be adequate. The set of

1200

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

nonlinear algebraic equations were then solved using a Newton solver NLEQ1. An alternate approach, which yielded identical steady-state profiles was to use as initial condition, m(z, 0)"0.02 and g(z, 0)"1.1, discretize into a set of ordinary differential equations in time using orthogonal collocation on finite elements and subsequently integrate until steady-state [using DDASAC (Caracotsios and Stewart, 1985)]. The advantage of this approach is that one can obtain a set of topos as well as the steady-state profiles with a single simulation. However, these topos do not explicitly account for the particular form of input, i.e. boundary or spatially distributed. They do provide a reasonable representation of the temperature profiles. In keeping with the spirit of the current work, the boundary topos were identified from the dynamic temperature profiles in response to a scaled step change of !0.02 in the feed temperature, which corresponds to 10 K. The sampling duration was five seconds. A total of hundred samples were taken. The first three of the identified boundary topos are shown in Fig. 18. The associated K-L singular values, which provide a measure of the importance of each mode, are presented in Table 6. The particular control problem addressed is a regulatory or disturbance rejection problem. The Damko¨hler number, a, is the ratio of the diffusion and reaction

time-scales. Thus with an increase in catalyst decay, the reaction time scale would increase and the Damko¨hler number would be expected to decrease. Considered here is the case of a decreasing from 0.875 to 0.8. In the absence of control, the hot spot is lowered and shifted significantly downstream, as shown in Fig. 19. Consequently, the conversion is dramatically reduced as manifested by the increase in reactant concentration at the exit of the reactor. A pseudo-modal controller, based on the dominant boundary topos alone, with a proportional gain of 10 (no integral action was needed) is able to maintain a temperature profile much closer to the original temperature profile. Although the DPS controller performs reasonably well in the above disturbance rejection problem, the spatial properties of the DPS controller are not highlighted. For example, one could obtain comparable performance from a simple LPS SISO proportional controller with judicious choice of the output location. Choosing the temperature at a location close to the hot spot as the

Fig. 18. Identified boundary topos for the nonlinear tubular chemical reactor.

Table 6 Identified K-L singular values for tubular reactor Mode i

K-L singular value l G

1 2 3 4 5

2.106;10\ 5.651;10\ 1.465;10\ 1.069;10\ 1.788;10\

Fig. 19. Regulatory temperature control for the nonlinear tubular reactor of Example 3. The disturbance arises due to the onset of catalyst decay. a decreases to 0.8 from 0.875. The DPS controller (K"10) is based on the dominant mode. The solid lines correspond to the setpoint profiles, long dashes to the controlled profiles and short dashes to the profiles in the absence of control.

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

output would enable the LPS SISO controller to successfully reject this disturbance. However, using a temperature measurement downstream (z"0.5) could lead to inverse response and poor controller performance. It is conceivable that the DPS controller would react better to spatially inhomogeneous disturbances. Another yardstick for the DPS controller would be its ability to stabilize complex spatiotemporal patterns. A boundary controller based on the ideas presented in this work is shown to achieve good results when applied to a reaction-diffusion system characterized by complex dynamical behavior (Chakravarti and Marek, 1995).

8. Summary Singular functions are fundamental properties that account for the spatial nature of distributed parameter systems. In principle, the singular functions are frequency dependent. From a practical standpoint, Gay and Ray (1995) found the steady-state singular functions to be adequate. A simple analysis in Section 3.1 showed that the steady-state singular functions, which correspond to the spatiotemporal operator K of the DPS, are equivalent QX to the singular functions of its spatial component K . X Previous work (Gay and Ray, 1995) showed that for the case of spatially distributed inputs, it was possible to design simple steady-state or dynamic input/output experiments to identify the steady-state right and left singular functions. Here we present another framework for the identification of singular functions of unknown distributed parameter systems based on dynamic input/output experiments. An important step here is the choice of basis functions used in projecting the spatiotemporal DPS to obtain a reduced temporal system. Section 3.2 motivated the use of the topos obtained by K-L decomposition of the dynamic output. Of course, the primary motivation for development of the dynamic identification scheme was its applicability to DPS with boundary inputs, for which case steady-state experiments yield little information. As was shown in Section 4, the dynamic identification scheme yields good estimates of the desired steady-state singular functions with only a small number of sensors for output measurement. It would also be helpful to notice that from a feedback control perspective, with boundary inputs, only the identification results concerning the left (output) singular functions are of significance. The right (input) singular functions are not used. This observation does add some value to the dynamic identification framework since estimation errors were much greater for the right singular functions. In any case, the condition number provided a convenient means of circumventing some of the difficulties associated with ill-posedness. For the examples considered in this work, a condition number of around 100

1201

seemed to provide a reasonable tradeoff between illconditioning and neglecting information. The current work also introduced the concept of boundary topos, which are an optimal basis for representation of the boundary Green’s function, g(z, t) or equivalently the step response S(z, t). The boundary topos are computed through the use of the K-L decomposition. Thus the K-L procedure has been employed in two distinct situations. Sections 6 and 7 illustrated how the empirically identified left (output) singular functions or boundary topos can be used as a basis for a pseudo-modal feedback control scheme. The range of static profiles achievable with boundary inputs for linear parabolic DPS was fairly limited. Thus, the case of time-varying setpoints was also considered.

Appendix A. Matrix singular value decomposition Singular value decomposition (SVD) of an arbitrary complex matrix K3Cn;m yields K"W RV*,

(A1)

where W and V are unitary matrices of order n;n and m;m, respectively, such that R"W*KV"diag(p ,2,p ), p"min[n, m],  N where p *p *2*p *0.   N R is referred to as the singular value matrix of K and the diagonal entries p are the singular values. W and V are G the left and right singular matrices respectively. Associated with each singular value p in R are the ith column G of W, denoted by w , and the ith column of V, denoted by G v , which are termed the ith left and right singular vectors G respectively. The singular values and vectors satisfy the following properties for i"1 to p: KK*w "pw , (A2a) G G G K*w "p v , (A2b) G G G K*K€ "pv , (A2c) G G G K€ "p w . (A2d) G G G The singular values are defined to be the positive square roots of the eigenvalues of the Hermitian matrices KK* and K*K and are hence unique. The left singular vectors w are the eigenvectors of KK* while the right singular G vectors v are the eigenvectors of K*K. Each pair of G singular vectors is unique with respect to each other only up to a scaling factor eHF of unit modulus. For example, the pair of singular vectors [w , v ] is completely equivaG G lent to [!w ,!v ]. In principle, any number of equivaG G lent pairs can be constructed by varying h.

1202

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

The K-L decomposition can be implemented through use of the SVD algorithm, in which case K corresponds to the dynamic output and could be thought of as a spatiotemporal data matrix. The n rows represent the spatial coordinate while the m columns correspond to the time coordinate. The left singular vectors w are then the G spatially dependent vectors or topos and the right singular vectors v are the time-dependent amplitudes or chroG nos. Though not explicitly demonstrated in this work, the matrix SVD can also be used to obtain the singular functions from the steady-state Green’s function (Gay and Ray, 1995). In this case, the n rows correspond to the steady-state output measurements while the m columns correspond to the m independent inputs. Standard routines in LINPACK (1989) or LAPACK (1992) can be used for numerical computation of the matrix SVD.

Appendix B. Functions of matrices The book by Golub and Van Loan (1989) provides a comprehensive discussion on the subject of functions of matrices. Several methods exist for approximating the function of a matrix, say h(A). One approach involves use of the Jordan canonical form. For the case of repeated eigenvalues, it is possible to express A3Cn;n as follows: A"M J M\, where J is referred to as the Jordan canonical form of A and is block diagonal. The columns of M correspond to the generalized eigenvectors of A. For the case of distinct eigenvalues, J reduces to a diagonal matrix. The diagonal entries correspond to the eigenvalues. Further the columns of M are now the eigenvectors of A. In such cases, M is referred to as the modal matrix. Based on the above Jordan characterization of a matrix, it is possible to express the function h(A) as follows: h(A)"M h(J) M\. Note that the columns of M are the eigenvectors of A and also h(A). For the case of exponential, logarithmic, sine and cosine functions, it is also possible to construct Taylor series approximations.

Appendix C. Computational framework for the K-L decomposition Given a spatiotemporal pattern or data set y(z, t), the K-L decomposition can be performed to obtain the topos, chronos and associated singular values (or the energies of the modes) through use of the stable matrix SVD algorithm described in Appendix A. However, there

could be some variations with respect to implementation based on how one views the distributed system. We consider the general case of output measurements being available only at a few finite locations over space but at regular (short enough) sampling intervals. The simple-minded approach would be to just take the raw data at each of the spatial locations over time and construct a ‘discrete’ spatiotemporal data matrix, Y"y(z , t ). One could then apply the SVD and obtain I J the left and right singular vectors as demonstrated in Appendix A. The left singular vectors correspond to the values of the topos at z , z ,2,z . This approach essen  , tially yields a discrete approximation of the ‘actual’ topos (continuous over the spatial domain, 0)z)1). The discreteness is not as crucial in the time domain since in practice by sampling fast enough, one can approximate the continuous-time behavior. The more rigorous approach would involve accounting for the spatial nature by reconstructing the entire profile from the few finite measurements. Application of the SVD would then directly yield the continuous topos. The measurements are still assumed to be discrete in time. It will be shown that this rigorous approach can be implemented in a relatively straightforward manner by the use of Gram matrices (Schumaker, 1981) and Cholesky weighting factors (Golub and Van Loan, 1989). The main ideas from the computational perspective are based on the work of Gay (1989). The problem at hand is: Given measurements of a particular output y at a finite number of spatial locations, say N, and over a time period with a small enough sampling rate, determine the topos and chronos. The truncated K-L expansion for y(z, t) is , y(z, t)+ t (z)l c (t). G G G G

(C1)

Since there are N measurement locations, we choose a set of N linearly independent trial basis functions for the output: +q (z),2,q (z),2,q (z),.  I , Corresponding to these output trial basis functions are a set of N linear functionals: +q*( ) ),2,q*( ) ),2,q*,  I , which are used to determine the coefficients b (t), I b (t)"q*(y(z, t)) I I The spatially distributed output y(z, t) can then be expressed in terms of the trial basis functions as follows: , y(z, t )+ q (z)b (t ), J I I J I

(C2)

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

where the subscript l refers to the sampling instant. The topos can also be expressed in terms of the trial basis functions: , t (z)" q (z)t , (C3) G H HG H where t is the coefficient of the jth trial basis function HG q (z) in the series approximation to the ith topos t (z). H G Substituting the series approximations for the output and topos from Eq. (C2) and (C3) into the truncated K-L expansion of Eq. (C1) yields: , , , (C4) q (z)b (t )" q (z) t l c (t ). I I J H HG G G J I H G Applying the Galerkin condition, with respect to the trial basis functions q (z), to Eq. (C4) results in the set of H N algebraic equations:







 , , q (z) q (z)b (t ) dz" q (z) q (z) F I I J F H   I H , t l c (t ) dz. (C5) HG G G J G This can be expressed in more succinct matrix notation as follows: QB"QWKC2,

(C6)

where the following notation has been used:



Q"





q (z)q (z) dz , F I

 B"+b (t ),, I J C2"+c (t ), G J W"+t ,, HG K"+diag(l ),. G By construction, the Gram matrix Q is symmetric positive-definite. The Cholesky decomposition of Q yields the Cholesky factor D (Q"D2D), an upper triangular matrix with positive elements. Thus, Eq. (C6) can be further reduced to: DB "D( "C2.

GHI Bª

GHI (ª

(C7)

BK can be ascertained from dynamic output measurements. Application of the SVD algorithm to Bª yields the left and right singular matrices, Wª and C, as well as the diagonal matrix K of singular values. The columns of the right singular matrix C correspond to the normalized chronos. The matrix of the trial basis coefficients of the topos W can be calculated from Wª as follows: W"D\WK .

(C8)

1203

The continuous topos can then be obtained through use of Eq. (C3). This completes the computational framework for the calculation of the (continuous) topos and chronos from dynamic data at a finite number of spatial locations.

References Alvarez, J., Romagnoli, J.A., & Stephanopoulos, (1981). Variable measurement structures for the control of a tubular reactor. Chem. Engng Sci., 36 (10), 1695—1712. Aubry, N., Guyonnet, R., & Lima, R. (1991). Spatiotemporal Analysis of Complex Signals: Theory and applications. J. Stat. Phys. 64, 683. Berkooz, G., Holmes, P. & Lumley, J.L. (1993). The Proper Orthogonal Decomposition in the analysis of turbulent flows. Annu. Rev. Fluid Mech. 25, 539—575. Butkovskiy, A.G. (1969). Distributed control systems, New York: American Elsevier. Butkovskiy, A.G. (1982). Green’s functions and transfer functions handbook. Chicester: Ellis Horwood. Caracotsios, M., & Stewart, W.E. (1985). Sensitivity analysis of initial value problems with mixed odes and algebraic equations, Comput. Chem. Engng., 9, 359—365. Carey, G.F., & Finlayson, B.A. (1975). Orthogonal collocation on Finite Elements. Chem. Engng. Sci., 30, 587—596. Chakravarti, S., Marek, M., & Ray, W.H. (1995). Reaction—diffusion system with Brusselator kinetics — control of quasiperiodic routes to chaos, Phys. Rev. E, 52 (3), 2407—2423. Chakravarti, S., & Ray, W.H. (1997). Boundary identification and control of distributed parameter systems using singular functions. In Proceedings of the American Control Conference, Albuquerque, NM (pp. 2223—2227). Chien-Chong Chen, & Hseuh-Chin Chang (1992). Accelerated Disturbance damping of an unknown distributed system by nonlinear feedback, A.I.Ch.E. J, 38 (9), 1461—1476. Chien-Chong Chen, Wolf, E.E., & Hseuh-Chia Chang (1993). Lowdimensional spatiotemporal thermal dynamics on nonuniform catalytic surfaces, J. Phy. Chem. 97, 1055—1064. Cochran, J.A. (1972). ¹he analysis of linear integral equations, New York: Mc-Graw-Hill. Fattorini, H.O. (1968). Boundary Control Systems, SIAM J. Control 6 (3), 349—385. Fukunaga, K. (1990). Introduction to statistical pattern recognition. Boston: Academic Press. Gay, D.H. (1989). Identification and control of distributed parameter systems by means of the singular value decomposition, PhD thesis, University of Wisconsin at Madison. Gay, D.H., & Ray, W.H. (1988). Application of singular value methods for identification and model based control of distributed parameter systems. In McAvoy, J.T. Arkun, Y. & Zzflriou, E. (Eds.), Proceedings of the IFAC ¼orkshop on Model-Based Process Control, Atlanta, GA (pp. 95—102). Oxford: Pergamon Press. Gay, D.H. & Ray, W.H. (1995). Identification and control of distributed parameter systems by means of the singular value decomposition. Chem. Engng Sci. 50 (10), 1519—1539. Golub, G.H., & van Loan, C.F. (1989). Matriz computations. Baltimore: Johns Hopkins University Press. Graham, M.D., Lane, S.L., & Luss, D. (1993). Proper Orthogonal Decomposition Analysis of Spatiotemporal Temporal Patterns, J. Phys. Chem. 97 (4), 889—894. Hanczyc, E.M., & Palazoglu, A. (1996). Use of symmetry groups in sliding mode control of nonlinear distributed parameter systems, Int. J. Control, 63 (6), 1149—1166.

1204

S. Chakravarti, W. Harmon Ray/Chemical Engineering Science 54 (1999) 1181—1204

LAPACK Users’ Guide (1992). LINPACK Users’ Guide (1989). Lorenz, E.N. (1956). Empirical orthogonal functions and statistical weather prediction. Technical report, MIT. Lumley, J.L. (1970). Stochastic tools in turbulence. New York: Academic Press. Naylor, A.W., & Sell, G.R. (1982). Linear operator theory in engineering and science. New York: Springer. Olivei, A. (1974). Boundary modal control of the temperature of the melt for crystal-growing in crucibles. Int. J. Control, 20 (1), 129—157. Olmstead, W.E. (1980). Boundary controllability of the temperature in a long rod, Int. J. Control, 31 (3), 593—600. Olmstead, W.E., & Schmitendorf, W.E. (1977). Optimal Control of the end-temperature in a semi-infinite rod. J. Appl. Math. Phys., 28, 697—706. Park, H.M., & Cho, D.H. (1996). The use of the Karhunen—Loe´ve decomposition for the modeling of distributed parameter systems, Chem. Engng Sci., 51(1), 81—98. Ramkrishna, D., & Amundson, N.R. (1985). ¸inear operator methods in chemical engineering. Englewood Cliffs. NJ: Prentice-Hall. Rigopoulos, A., Arkun, Y., & Kayihan, F. (1997). Identification of full profile disturbance models for sheet forming processes, A.I.Ch.E. J. 43(3), 727—739.

Schumaker, L.L. (1981). Spline functions: Basic theory. New York: Wiley-Interscience. Sirovich, L. (1987). Turbulence and the dynamics of coherent structures—parts I,II,III, Quart. Appl. Math. X¸» (3), 561—590. Smithies, F. (1958). Integral equations. Cambridge: University Press. Stakgold, I. (1979). Green’s functions and boundary value problems. New York: Wiley. Stewart, G.W. (1993) On the early history of singular value decomposition. SIAM Rev., 35(4), 551—566. Villadsen, J., & Michelsen, M.L. (1978). Solution of differential equation models by polynomial approximation. Englewood Cliffs, NJ: Prentice-Hall. Villadsen, J.V., & Stewart, W.E. (1967). Solution of boundary-value problems by orthogonal collocation. Chem. Engng Sci., 22, 1483—1501. Windes, L.C. (1986). Modelling and control of a packed bed reactor. PhD thesis, University of Wisconsin at Madison. Zhou, X.G., Liu, L.H., Dai, Y.C., Yuan, W.K., & Hudson, J.L. (1996). Modeling of a fixed-bed reactor using the K-L expansion and neural networks. Chem. Engng Sci., 51 (10), 2179—2188. Zolesio, J.P. (Ed.), (1991). Boundary control and boundary variation. New York: Springer.