- Email: [email protected]
Controller design for parabolic distributed parameter systems using finite integral transform techniques
J Proc Cont. Vol 6, No 6, p p 359-366, 1996 Copyright ,g 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0959-1524/96 $15.00 + 0.00
ELSEVIER
0959-1524(95)00036-4
Controller design for parabolic distributed parameter systems using finite integral transform techniques Masatoshi Yoshida and Shigeru Matsumoto* Department of Biochemistry and Engineering, Tohoku University, Aramaki, aza Aoba, Aoba-ku, Sendai 980, Japan A procedure is presented for designing a control system for distributed parameter systems of parabolic type based on the reduced-order decoupled state-space model obtained by a finite integral-transform technique. A Kalman filter is used as an observer to estimate the state variables, and state feedback control is performed. The method was applied to a one-dimensional heat conduction process and a moving bed adsorber. Both state estimation and control performances were satisfactory in spite of the model and parameter uncertainties. Following this controller design approach, the searching algorithms for the optimal sensors' and the optimal actuators' allocation problems were solved. These algorithms were applied to a one-dimensional heat conduction process in order to confirm the state estimator and controller performance. The fastest state estimation could be achieved by assigning the sensors at the optimal locations and the desired state distribution was realized with a few actuators located at the optimal positions. Copyright qg 1996 Elsevier Science Ltd
Keywords: distributed parameter system; finite integral transform; heat conduction process A large n u m b e r o f chemical processes involving heat and/or mass transfer or chemical reaction are modelled by partial differential equations o f parabolic type. Hence the design o f an effective controller for such systems is of considerable importance in industrial practice. F o r chemical processes, cases are often met where the dynamic behaviour is well u n d e r s t o o d and the parameters are substantially known. In such cases, it is expected that better control p e r f o r m a n c e can be achieved by making the best use o f knowledge on the process to be controlled. F r o m this point o f view, we have p r o p o s e d a procedure for designing a control system based on the physical model o f distributed parameter system (DPS) w-.
for0
y = 4hD/pCpd 5 = l/pCp In Equation (l), z denotes the dimensionless distance based on the length l, and Q(z, t) the distribution o f heat generation rate. The b o u n d a r y conditions are assumed to be the third kind:
3T -& -+ 3T
Consider a one-dimensional heat c o n d u c t i o n system which involves heat generation and heat removal to the surroundings. If the thermal properties are assumed to be constant, the d y n a m i c behaviour is expressed by:
at
c~ 32T l 2 07,2
7(T-
T=) +
6Q(z,t
at z = 0
Bi~T= a t z = l
(2)
(3)
where Bi~ and Bin are Biot numbers at z -- 0 and z -- 1 respectively. N o w , consider the control problem to keep the temperature T(z, t) at a desired profile 7"*(z) by manipulating the heater outputs Q(z, t). The above model can be transformed to a lumped parameter system (LPS) by use of the finite Fourier transform technique 3. The Fourier transform o f T(z, t) is defined over the interval z -- 0-1 as follows:
Algorithm Jar controller design
-
BiiT = BiIT~
--+Bi~T= 9z
Theory
3T
1, t > 0 , where:
(1)
*To whom correspondence should be addressed Paper presented at the 5th International Symposium on Process Systems Engineering, Kyongju, Korea, 30 May-3 June 1994
7{T(z,t)} = Ti(t ) = I~X(~,, z ) T ( z , t)dz
359
(4)
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
360
where X(~, z) is an eigenfunction which satisfies the homogeneous boundary condition.
(11)
y = Ccx where:
X(~i, z) = ~ cos ~,.z + Bi I sin ~,.z
(5) y = [T(z I , t),..., T ( Z j , t) ..... T ( z m, t)] x
The eigenvalues ~ are the positive roots of:
C¢ =[Oji],
tan ~i = ~, ( Bi 1 + Bi2) /(~ f - BilBi2)
(6)
The inverse formula is Y-'{T,.(t)} - T ( z , t ) = ~, X ( ¢ " z ) ~ ( t ) i=l N ( ~ )
(7)
where N(~) is the norm of eigenfunction defined by
N(~,) =
( ~ + Bi~) 1 -t ~ + Bi~ ÷ Bi,
(8) Applying the finite Fourier transform to Equations (1)-(3), the following ordinary differential equations are obtained: d~ dt -
o¢ 2 ~i + Y
-
-
+ 5Q~ + diT~
(9)
for i = 1,2 ..... o% where:
[Bil~i + BizX(~i,1)],
d i -- [ 0 ¢ / l 2 4- ~ / / ~ 2 ] .
Q, = F { Q ( z , t ) } In the above equations, T, are independent of each other with respect to eigen-mode i, and hence fully decoupled. Neglecting_modes of_eigenvalues higher than n, and representing T,(t) and Q,(t), respectively, as state variables and manipulated variables, we have the following state-space expression: = AcX +Bcu 4-dcw
(10)
where
x=[~,~ u = [
.....
T, . . . .
T.] T
,Q2 . . . . .
W = T~
A c = diag {a l, a 2 ..... a i .... a,} B c = 5I d~ = [d l, d2,..., d i .... dn] v a, = -(o¢~/l 2 + Y) for i = 1..... The cut-off eigen-mode n should be selected by taking into account the time constants v~(= -1/ai) of the eigen mode, which decrease according to the order of the eigen mode, and the possible minimum sampling interval which is subject to the limitation of the hardwareL If measurements are made at m points z = Zl,...,Zm, and if the observed temperatures are regarded as output variables, the following output equation is obtained by use of the Fourier-inverse transform:
( J = l ..... m ; i = l ..... n)
(pj, = X ( ~ i, zj )/N(~i ) Based on the state-space model Equations (10) and (11), the feedback control system can be constructed as shown in Figure 1. Since both matrices Ac and Bc are diagonal, t_here is no coupling between the manipulated variable Qj and state variable T,. except when j -- i. Therefore, the controller in Figure 1 can be designed for a SISO system. As the state variables cannot be measured directly, they must be estimated by an observer. In this work, a Kalman filter will be adopted for the reason of the compensation of the model uncertainty and the observation noise. In practice, this control system must be realized in a discrete system. If the discretization is done on the basis of a sampling time At, the state-space model becomes: x(k + 1) = Ax(k) + Bu(k) + w y (k)
(12)
y(k) = Cx(k) + VN(k )
(13)
In the above equations, the system noise WN and the measurement noise VN are introduced instead of the disturbances. These noises are assumed to be independent and white Gaussian noise, and their covariances are denoted by Q and R respectively. The estimation algorithm by the current-type Kalman filter is given by: x(k + l) = Ax(k) + Bu(k)
(14)
x(k + 1) = x(k) + K(y(k + 1) - Cx(k + 1))
(15)
K -- APCT[R + CPCT] 1
(16)
where P is a covariance matrix of estimation error and the solution of the following Riccati equation: p = Q + APA T - APCT[R + CPCT]-~CPA T
(17)
Searching algorithm Jor optimal sensor allocation The estimation performance is affected by the allocation of the observation points. It is preferable that sensors are located at which the state estimation is carried out as precisely and fast as possible. Such sensor allocation is called optimal sensor allocation. If we want to achieve the smallest estimation error, the following objective function should be minimized:
J =~e(k)Ve(k) k=O
where:
(18)
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
361
Manipulated variables +
Des i red temp. T(z,t)
×l
I el
_
>
Gclcsl
÷
i
Heater )ower
_1
I
x21
52
Plant
1 ~ _ ~ UCntS) ]
HI Estimated temp. /% T(z,t)
4--
2 -4 • ®
-4 Figure 1
I
L
T(z l,t)
I
T(znv t )
Kalman filter
Tn Meas. temp.
~ Estimated
state variables
Structure of the control system
e(k) = x ( k ) - x(k)
Searching algorithm f o r optimal actuator allocation
In case of a Kalman filter, the above objective function can be approximated by the covariance matrix of estimation error with initial error vector (refer to Appendix A): Y = e(0)Tpe(0)
(19)
In order to search the optimal allocation of m sensors, we define the sensor allocation vector zs as: ZS ---- [/1'/2 . . . . . Zm]T
(20)
In order to keep the temperature at a specified target profile T*(z) by using a limited number of heaters, the heaters must be allocated at suitable locations so that the best control performance can be achieved. Consider the allocation problem of r heaters with the width of 20" as shown in Figure 2. We define here the actuators' allocation vector z, and heater output vector q as follows: Z a = [Z1,Z2,...,Zj .....
Zr]T, (24)
q = [ql,q2 ..... G ..... q,]T
if sensor positions are changed, the coefficient matrices C of the output equation and P also change. Therefore, the objective function depends on the sensor allocation vector• We use a conjugate gradient method to search the minimum value of the objective function. The gradient vector is defined as follows: r
aJ_lOJ c)zs
oJ
aJ
L&, ' Oz~ '"" Ozj
M az m
OJ - e(0)T_~ - ~ = e(0)
where
(22)
dz I
OP = S ~_i sT _ K _ _ ~ p s T o~CT Oz/ 0z/ - SP ~ KT
u = Hq
H = [~',j]
(i = 1..... n;
(,7
~i, z) dz
(23)
We consider the following objective function to be minimized: I
J = Io(Ts(z) - T * ( z ) ) 2 d z
(26)
where T~(z) is the steady-state temperature realized by r heaters. Denoting the Fourier transform of T l z ) and T*(z), respectively by T~,i and T*, the above function can be rewritten as follows (refer to Appendix C): J = f01T * ( z ) 2 d z - 2 x
where S=A-KC
j = 1..... r) (25)
(21)
The elements of the gradient vector are expressed in the form of a Lyapunov equation as shown in Equation (23) (refer to Appendix B), and then they can be calculated numerically.
~z/
If each heater is assumed to have a uniform heat generation rate, the manipulated variable vector u is related to q by use of a Fourier transform.
*T N x s + x ~ N x
s
(27)
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
362
L.
1 T(z,t)
I ~ d
'7" 7" I I
I I
I I
I I
I I
I I
0
z=O
Figure2 X*
=
D
z1
--~
I I
I I
! I
I I
~
..........
'"hi I
2o
zj
t
~
I I
.................
I
t
Heater
I
z=l
Zr
Schematic diagram of one-dimensional heat conduction system [,F~,..
Experimental
-, T .,T~]
a p p a r a t u s and p r o c e d u r e
x~ = - A~I(BcHq + d0v) Simple heat conduction process
N = diag{N(~l)~,N(~2)-' ..... N(~,)-'}
Figure 3 shows the experimental apparatus. The process to be controlled is a metal rod made of stainless steel 0.48 m long and 0.3 m in diameter. The temperature profile along the rod length is controlled by manipulating 20 band heaters (with a maximum power of 50 W) wound around the rod. The temperature was measured at 11 points by copper-constant thermocouples. A personal computer was used for execution of the estimation and control algorithm with a sampling time of 20 s. The PI control law was used for the controller in the proposed control scheme. The values of physical parameters in the mathematical model, shown in Table 1, are quoted from the literature except for the heat transfer coefficients, which were determined by preliminary experiments.
7"~,i is an element of vector x~ which is the steady-state solution of Equation (10). Since the coefficient matrix H is a function of the actuators' allocation vector Za, the value of J from Equation (27) also depends on both Za and q. In this paper, the gradient method is used to find the minimum value of the objective function. The gradient vector of J can be formulated as shown below:
~...~
aJ -a- q- - ~
[~JaJ aJ ~lql" ' o~q2 ' ' ' ' ' aq j
(28)
,...,
'''''
aJ C~qr
(29)
Moving bed adsorber
The elements of these gradient vectors are derived by differentiation of Equation (27) with respect to z a and q as follows:
OJ = 2(x,T_xT)NAcIB ~ o~I -~zjq
0zj
Figure 4 shows the experimental apparatus of the moving bed adsorber, which consists of a stainless steel column of 1.2 m high, with outside diameter 0.03 m and inside diameter 0.025 m. Adsorbent particles were fed in
Table1
(30)
cgJ = 2(x,T_xT)NAEIBcH &l cgqj Oqj
(m2/s) ( W / ( m . K)) hD (W/(m 2 " K)) h~ (W/(m 2. K)) h: (W/(m 2 . K)) 7 ( s 1)
(31)
Power24Ch Controller]
D,I 2t
Thermocouple
i
4.44 × 10 6 16.3 40.3 239 232 1.3 × 10 3
GPIB
1,1 ,91
Fgil
I
[ t~gi~
Datalogger
Figure3
Parameters for the experimental process
Experimental apparatus of simple heat conduction process
GPIB
Co~~
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
a Adsorber
. Manometer
b Hopper
Digiufl-multimeter
c Rotary-valve
Personal computer
d Blower
363
:Thermocouple
e Feed Band heater f Preheater a Power controller
g Seal gas
m
00000 00000 00000
D
d
Figure 4
Experimental apparatus of moving-bed adsorber
to the top of the adsorber from a hopper and preheated air was introduced from the bottom. Activated carbon was used as an adsorbent and adsorption of the moisture in the air occurred in the adsorber. Although the thermal dynamics of the system was complicated, it was modelled by the following simple linear partial differential equation, and the control system was designed in accordance with the procedure mentioned above.
DT at
a O:T /: &2
fl aT 1 Oz
Results and discussion
y(T- T=)+ag+aq
(32)
where C c = {l -- e ) C p p p p
{(1
+
eCpgpg}S + CPoPoSo
e)Zp +ek.g}S+k.oS o
G¢ =
Ce Ce
y-
10 s. The cut off eigen-mode was set at 7. The thermal properties in the mathematical model were quoted from the literature, and heat transfer coefficients were obtained from the step response of the temperature profile in the adsorber. The values are shown in Table 2. The control experiments were carried out for the various numbers of observation points in order to test the state estimator performance.
hm/ Ce
Simple heat conduction process Figure 5 shows an example of controller performance, in the form of the transient states from a certain initial state to the target distribution indicated by a broken line. The measurements were carried out at : = 0, 0.5 and 1. Solid lines indicate the estimated temperature profiles. Fairly good estimation was achieved only with a few points of observation and the temperature distribution coincides with the target distribution after 2400 s. Table 2
Thermal properties and mass flow rate
1 Ce
Control of axial temperature distribution was performed by manipulating 20 100 W band heaters wound around the adsorber. A personal computer was used for the measurements and the regulation of the power controller, and the time interval for control was set at
Air Activated carbon Stainless steel
2
('p
p
If"
( W / ( m . K})
(J/(kg- KI)
(kg/nf)
/kg/s)
0.03 1.5 21.5
1017 75~ 461
109 ~50 7700
7.27 × 10 5.33 × 10 4
Void fraction t Heat transfer coefficient h (W/(m ~ - K))
0.35 2(1
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
364
1.0
140 '
'
110
' Stainle's steel t = 2400.
•
i
0.8
L
i
Error- e(t)Te(t)
L
~¢_. 0.6
Z = [0.1,0.9] T Z = [0.4,0.6]T
0
L0.4
L.LJ
1-
0.2
50
I
0
200
0 •
200
•
f
,
0_2
°
0.4
5
I'
0 6
z
Figure
i
~(o)r~(o)
t = 120~N~
~ so
i
(a)
400
•
O.
1
(b)
Controller performance for the simple heat conduction
process
IOOC
BOO
Time [s]
400. 300
[-3
600
~- I 0 0 ~
Oi 0
1.0
03 C3
I
I
I
I
0.2
0.4
I
I
I
0.6
I
0.8
1.0
zF] l I
/ \
,
/
i
/
i
/
I
Figure 7
Effect o f sensor location on the state estimator performance: (a) transition o f estimation error; (b) target profi]e
~ - i
250 ~o
r I
(-xJ
-~I
, /,,/4#~ I
0
i.
0.2
~
0.4
I
;; .4 /i
0.6
/
0.8
1.0
.~
150
1St SENSORPOSITION o:iniUal
Figure 6
point
o:opUmal point
Convergence to optimal sensor locations for two sensors
5O
Figure 6 shows the result of the optimal sensor allocation with two sensors searched for the transient state toward the stead-state profile as shown in Figure 7(b). The white circles in the figure indicate the initial positions, the black points are the optimal points, and the solid lines are the trajectories of searching. The result shows that the same optimal position can be found despite different initial positions. The dashed lines are contour plots of the objective function of Equation (19). Figure 7(a) indicates the trajectories of estimation error for various pairs of observation points during the transient state toward the steady-state profile as shown in Figure 7(b). In this case, the optimal sensor allowcation was obtained as z~ = [0.27, 0.73] x. Obviously, the fastest speed of estimation was obtained when the sensors were located at the optimal allocation. An optimal allocation problem was solved for the case that the temperature distribution of stainless steel plate 48 cm long was controlled to a desired profile by means of several plate-heaters 2.5 cm wide with a maximum generation rate of 65 W. Figure 8 shows the steady-state temperature distributions attained by various numbers of heaters located at the optimal locations. The dashed line indicates the target profile. It is found that the target profile can be realized with five actuators.
Moving bed adsorber The control result is shown in Figure 9, where the closed response of the temperature distribution is indicated from the initial state. In this case, the state estimation
0
~
0.0
i
I
0.2
I
I
I
0.4
I
0.6
I
I
0.8
|
1.0
zH Figure 8 Steady state temperature profile achieved by the actuators which were optimally located was carried out by using 11 temperature measurements at the marked positions in the figure (the positions were not at optimal allocation). In Figure 9, it can be seen that the temperature distribution was estimated well except for the initial period and controlled satisfactorily within an error of + 5°C after 2500 s. Figure 10 shows the heater-output distributions corresponding to Figure 9. Zero heater-output around the middle of the column is considered to be the effect of disturbance due to the heat of adsorption.
Conclusion A new controller design algorithm for parabolic DPS is proposed and applied to a simple heat conduction process and a moving bed adsorber. The results show that the state distribution can be estimated accurately by using only a few sensors, and controller performance is also satisfactory. Furthermore, methods for searching the optimal sensors' and actuators' allocation are presented, and it is demonstrated that the state estimation and control performances are improved by assigning sensors and actuators at the optimal locations. The proposed control scheme is basically similar to the conventional modal control structure 4.5. However, it
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
170~
V
V
V
V
aoL
V
V
V
V
Since S is a stable matrix, e(k) converges to zero as k increases. The objective function J of Equation (18) can be rewritten using the initial error e(0) as follows:
V °
t:25oo
J : e(0)Tp'e(0)
gOF
olo~
"
----0.4
0.6 z
P' = I + STp'S
o
0.8
100
. -
.
estimated distribution, o measurements) .
.
.
.
.
.
.
_j
Now, P, the solution of the Riccati equation Equation (17), can be rearranged in the form of a Lyapunov equation: p = sPSr+
~" so
0.0
(A.3)
.0
[-]
Figure 9 Controller performance for the moving-bed adsorber ( - - - , target distribution, - - - - ,
is a solution of the following Lyapunov
j t =os__ ~
0.2
(A.2)
• where P' equation:
10 0.0
365
Q + KRK v
(A.4)
P' is related to P as shown below by calculating the trace of the above equation.
0.2
0.4
0.6 z
0.8
1.0
[-]
traceP = trace{V' + SV'S T + S:V'S T2 + ...} = trace{V'(l + STS + STeS2 + ..3} = trace{V'P' }
(A.5)
Figure 10 Distributions of heater outputs during the control in Figure 9 where: is essentially different in the formulation of state-space representation, where it has several advantages compared with the conventional model control approach. The most distinctive one is that the state variables can be estimated directly from a few observations using a Kalman filter and the whole state distributions are obtained numerically from the data of state variables. Furthermore, boundary conditions are involved in the state equation and consequently this approach can be also extended to the controller design of the boundary control problem.
V'
= Q
+ KRK r
In this study, the covariance matrices are assumed to be diagonal matrices with constant coefficients such as Q = 0.1I, R = 0.1I. Therefore, V' can be assumed as the coefficient matrix. Consequently, P' in Equation (A.2) may be replaced by P, as the form of the objective function is not affected.
Appendix B: Derivation of Equation (23) Differentiating Equation (15) with respect to sensor location z/
References 1 Matsumoto, S. and Yoshida. M. Kagaku Kogaku Ronhunshu, 1987, 13, 145; English translation in Int. Chem. Eng. 1989, 29, 158 2 Yoshida, M., Matsumoto, S. et al. Kagaku Kogaku Ronbunshu, 1988, 14, 147; English translation in Int. Chem. Eng. 1990, 30, 142 3 Ozisik, M. N. 'Heat Conduction', John Wiley, New York, 1980 4 Lausterer, G K. and Eitelberg, E, in 'Distributed Parameter Control Systems', (Ed. S. G. Tzafestas), Pergamon Press, Oxford, 1981 5 Ray, W. H. and Lainitois, D. G. 'Distributed Parameter Systems Identification. Estimation and Control', Marcel Dekker, New York, 1976
3P = A 0P A - A 0z/ &j - AP
c3CT
3P C r [ R + C P C T ] I C P A 8:j [R + CPC T]
_ APCT 0[R + CPC T] i CPA APC T [ R + C P C
T] I OK?pA
c): I 3t' - A P C T [ R + C P C T] I C ~ A
Appendix A: Derivation of Equation (19) When a steady state Kalman filter is used as a state estimator, the dynamics of state estimation error e can be expressed as follows: e(k + 1) = Se(k)
I CPA
o3z/
(A.1)
(B.1) Since [R + CPCT]q[R + CPC T] = 1, the following relation can be derived.
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
366
3[R + C P C v ]-:
Wf'X({*'z)X(;'z)dzl = ~ J0
O,
for i * j for i = j
N({,), 3[R + C P C v ] = - [ R + C P C v ]-:
• [R + CPC v ]q
Oz/
(C.4)
Then the third term o f the right-hand-side in Equation (C.2) can be simplified as:
ac
= - [R + C P C v ] -I oqz7 P C v [R + C P C s l-'
s,i i=l N(:,)
- [R + cPCT ]-Ic--~ cT [R +cPCT] -1
'
*=l
l
z) 2 dz '
~2
¢lzj
s,i
- [R + C P C T ] - 1 C P - ~ [ R
i=l
+ c P C T ] -1
OZ,
(B.2) R e a r r a n g i n g E q u a t i o n (B.1) with E q u a t i o n (B.2) and the gain matrix K in E q u a t i o n (16), E q u a t i o n (23) is
If eigenvalue modes greater than n are negligible, T~.~ and ~ become elements o f the state variable of x s and x* respectively. Therefore, E q u a t i o n (27) can be derived by substituting Equations (C.3), (C.5) into Equation (C,2) and by using these state variables.
derived.
Appendix C: Derivation of Equation (27)
Nomenclature
By using the finite Fourier transform, the steady state
Bi Cp
temperature Ts(z ) can be expressed as:
T~(z):~X(~i'z)~, i=l N ( { , )
(C.1) '
Substituting the above equation into E q u a t i o n (26), we have:
j:
i ~ T , ( z ) Z d z - 2 i i T , ( z ) ~ " X(~,,z) ~,.dz
i=1 N(~i) -
X(~i'z) ~ i X(¢,)
'
dz
'
(C.2)
The second term o f the right-hand side in E q u a t i o n (C.2) can be rewritten as follows:
,'=1
N(~i)
d g h l N P Q q S T T~ t u W w X x y z a t A p
Biot number , hl/)~ specific heat (J/kg . K) diameter (m) heat generation rate (W/m) heat transfer coefficient (W/mz • K) length (m) norm of eigenfunction covariance matrix of estimation error heat generation rate (W/m~) heat of adsorption (W/m) cross-sectional area (m:) temperature (K) surrounding temperature (K) time (s) vector of manipulated variable mass flow rate (kg/s) disturbance eigenfunction vector of state variable vector of output variable dimensionless distance thermal diffusivity (mZ/s) void fraction thermal conductivity (W/m • K) eigenvalue density (kg/m3)
Superscripts
'
Fourier transform
= 2~ ~ J o T * ( z ) X ( ~ i , z ) d z = 2
~,:, ~
5 (c.3)
The kernel function is an o r t h o g o n a l function expressed as follows:
Subscripts g i o p
gas phase mode of eigenvalue outer casing particle phase
ELSEVIER
0959-1524(95)00036-4
Controller design for parabolic distributed parameter systems using finite integral transform techniques Masatoshi Yoshida and Shigeru Matsumoto* Department of Biochemistry and Engineering, Tohoku University, Aramaki, aza Aoba, Aoba-ku, Sendai 980, Japan A procedure is presented for designing a control system for distributed parameter systems of parabolic type based on the reduced-order decoupled state-space model obtained by a finite integral-transform technique. A Kalman filter is used as an observer to estimate the state variables, and state feedback control is performed. The method was applied to a one-dimensional heat conduction process and a moving bed adsorber. Both state estimation and control performances were satisfactory in spite of the model and parameter uncertainties. Following this controller design approach, the searching algorithms for the optimal sensors' and the optimal actuators' allocation problems were solved. These algorithms were applied to a one-dimensional heat conduction process in order to confirm the state estimator and controller performance. The fastest state estimation could be achieved by assigning the sensors at the optimal locations and the desired state distribution was realized with a few actuators located at the optimal positions. Copyright qg 1996 Elsevier Science Ltd
Keywords: distributed parameter system; finite integral transform; heat conduction process A large n u m b e r o f chemical processes involving heat and/or mass transfer or chemical reaction are modelled by partial differential equations o f parabolic type. Hence the design o f an effective controller for such systems is of considerable importance in industrial practice. F o r chemical processes, cases are often met where the dynamic behaviour is well u n d e r s t o o d and the parameters are substantially known. In such cases, it is expected that better control p e r f o r m a n c e can be achieved by making the best use o f knowledge on the process to be controlled. F r o m this point o f view, we have p r o p o s e d a procedure for designing a control system based on the physical model o f distributed parameter system (DPS) w-.
for0
y = 4hD/pCpd 5 = l/pCp In Equation (l), z denotes the dimensionless distance based on the length l, and Q(z, t) the distribution o f heat generation rate. The b o u n d a r y conditions are assumed to be the third kind:
3T -& -+ 3T
Consider a one-dimensional heat c o n d u c t i o n system which involves heat generation and heat removal to the surroundings. If the thermal properties are assumed to be constant, the d y n a m i c behaviour is expressed by:
at
c~ 32T l 2 07,2
7(T-
T=) +
6Q(z,t
at z = 0
Bi~T= a t z = l
(2)
(3)
where Bi~ and Bin are Biot numbers at z -- 0 and z -- 1 respectively. N o w , consider the control problem to keep the temperature T(z, t) at a desired profile 7"*(z) by manipulating the heater outputs Q(z, t). The above model can be transformed to a lumped parameter system (LPS) by use of the finite Fourier transform technique 3. The Fourier transform o f T(z, t) is defined over the interval z -- 0-1 as follows:
Algorithm Jar controller design
-
BiiT = BiIT~
--+Bi~T= 9z
Theory
3T
1, t > 0 , where:
(1)
*To whom correspondence should be addressed Paper presented at the 5th International Symposium on Process Systems Engineering, Kyongju, Korea, 30 May-3 June 1994
7{T(z,t)} = Ti(t ) = I~X(~,, z ) T ( z , t)dz
359
(4)
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
360
where X(~, z) is an eigenfunction which satisfies the homogeneous boundary condition.
(11)
y = Ccx where:
X(~i, z) = ~ cos ~,.z + Bi I sin ~,.z
(5) y = [T(z I , t),..., T ( Z j , t) ..... T ( z m, t)] x
The eigenvalues ~ are the positive roots of:
C¢ =[Oji],
tan ~i = ~, ( Bi 1 + Bi2) /(~ f - BilBi2)
(6)
The inverse formula is Y-'{T,.(t)} - T ( z , t ) = ~, X ( ¢ " z ) ~ ( t ) i=l N ( ~ )
(7)
where N(~) is the norm of eigenfunction defined by
N(~,) =
( ~ + Bi~) 1 -t ~ + Bi~ ÷ Bi,
(8) Applying the finite Fourier transform to Equations (1)-(3), the following ordinary differential equations are obtained: d~ dt -
o¢ 2 ~i + Y
-
-
+ 5Q~ + diT~
(9)
for i = 1,2 ..... o% where:
[Bil~i + BizX(~i,1)],
d i -- [ 0 ¢ / l 2 4- ~ / / ~ 2 ] .
Q, = F { Q ( z , t ) } In the above equations, T, are independent of each other with respect to eigen-mode i, and hence fully decoupled. Neglecting_modes of_eigenvalues higher than n, and representing T,(t) and Q,(t), respectively, as state variables and manipulated variables, we have the following state-space expression: = AcX +Bcu 4-dcw
(10)
where
x=[~,~ u = [
.....
T, . . . .
T.] T
,Q2 . . . . .
W = T~
A c = diag {a l, a 2 ..... a i .... a,} B c = 5I d~ = [d l, d2,..., d i .... dn] v a, = -(o¢~/l 2 + Y) for i = 1..... The cut-off eigen-mode n should be selected by taking into account the time constants v~(= -1/ai) of the eigen mode, which decrease according to the order of the eigen mode, and the possible minimum sampling interval which is subject to the limitation of the hardwareL If measurements are made at m points z = Zl,...,Zm, and if the observed temperatures are regarded as output variables, the following output equation is obtained by use of the Fourier-inverse transform:
( J = l ..... m ; i = l ..... n)
(pj, = X ( ~ i, zj )/N(~i ) Based on the state-space model Equations (10) and (11), the feedback control system can be constructed as shown in Figure 1. Since both matrices Ac and Bc are diagonal, t_here is no coupling between the manipulated variable Qj and state variable T,. except when j -- i. Therefore, the controller in Figure 1 can be designed for a SISO system. As the state variables cannot be measured directly, they must be estimated by an observer. In this work, a Kalman filter will be adopted for the reason of the compensation of the model uncertainty and the observation noise. In practice, this control system must be realized in a discrete system. If the discretization is done on the basis of a sampling time At, the state-space model becomes: x(k + 1) = Ax(k) + Bu(k) + w y (k)
(12)
y(k) = Cx(k) + VN(k )
(13)
In the above equations, the system noise WN and the measurement noise VN are introduced instead of the disturbances. These noises are assumed to be independent and white Gaussian noise, and their covariances are denoted by Q and R respectively. The estimation algorithm by the current-type Kalman filter is given by: x(k + l) = Ax(k) + Bu(k)
(14)
x(k + 1) = x(k) + K(y(k + 1) - Cx(k + 1))
(15)
K -- APCT[R + CPCT] 1
(16)
where P is a covariance matrix of estimation error and the solution of the following Riccati equation: p = Q + APA T - APCT[R + CPCT]-~CPA T
(17)
Searching algorithm Jor optimal sensor allocation The estimation performance is affected by the allocation of the observation points. It is preferable that sensors are located at which the state estimation is carried out as precisely and fast as possible. Such sensor allocation is called optimal sensor allocation. If we want to achieve the smallest estimation error, the following objective function should be minimized:
J =~e(k)Ve(k) k=O
where:
(18)
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
361
Manipulated variables +
Des i red temp. T(z,t)
×l
I el
_
>
Gclcsl
÷
i
Heater )ower
_1
I
x21
52
Plant
1 ~ _ ~ UCntS) ]
HI Estimated temp. /% T(z,t)
4--
2 -4 • ®
-4 Figure 1
I
L
T(z l,t)
I
T(znv t )
Kalman filter
Tn Meas. temp.
~ Estimated
state variables
Structure of the control system
e(k) = x ( k ) - x(k)
Searching algorithm f o r optimal actuator allocation
In case of a Kalman filter, the above objective function can be approximated by the covariance matrix of estimation error with initial error vector (refer to Appendix A): Y = e(0)Tpe(0)
(19)
In order to search the optimal allocation of m sensors, we define the sensor allocation vector zs as: ZS ---- [/1'/2 . . . . . Zm]T
(20)
In order to keep the temperature at a specified target profile T*(z) by using a limited number of heaters, the heaters must be allocated at suitable locations so that the best control performance can be achieved. Consider the allocation problem of r heaters with the width of 20" as shown in Figure 2. We define here the actuators' allocation vector z, and heater output vector q as follows: Z a = [Z1,Z2,...,Zj .....
Zr]T, (24)
q = [ql,q2 ..... G ..... q,]T
if sensor positions are changed, the coefficient matrices C of the output equation and P also change. Therefore, the objective function depends on the sensor allocation vector• We use a conjugate gradient method to search the minimum value of the objective function. The gradient vector is defined as follows: r
aJ_lOJ c)zs
oJ
aJ
L&, ' Oz~ '"" Ozj
M az m
OJ - e(0)T_~ - ~ = e(0)
where
(22)
dz I
OP = S ~_i sT _ K _ _ ~ p s T o~CT Oz/ 0z/ - SP ~ KT
u = Hq
H = [~',j]
(i = 1..... n;
(,7
~i, z) dz
(23)
We consider the following objective function to be minimized: I
J = Io(Ts(z) - T * ( z ) ) 2 d z
(26)
where T~(z) is the steady-state temperature realized by r heaters. Denoting the Fourier transform of T l z ) and T*(z), respectively by T~,i and T*, the above function can be rewritten as follows (refer to Appendix C): J = f01T * ( z ) 2 d z - 2 x
where S=A-KC
j = 1..... r) (25)
(21)
The elements of the gradient vector are expressed in the form of a Lyapunov equation as shown in Equation (23) (refer to Appendix B), and then they can be calculated numerically.
~z/
If each heater is assumed to have a uniform heat generation rate, the manipulated variable vector u is related to q by use of a Fourier transform.
*T N x s + x ~ N x
s
(27)
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
362
L.
1 T(z,t)
I ~ d
'7" 7" I I
I I
I I
I I
I I
I I
0
z=O
Figure2 X*
=
D
z1
--~
I I
I I
! I
I I
~
..........
'"hi I
2o
zj
t
~
I I
.................
I
t
Heater
I
z=l
Zr
Schematic diagram of one-dimensional heat conduction system [,F~,..
Experimental
-, T .,T~]
a p p a r a t u s and p r o c e d u r e
x~ = - A~I(BcHq + d0v) Simple heat conduction process
N = diag{N(~l)~,N(~2)-' ..... N(~,)-'}
Figure 3 shows the experimental apparatus. The process to be controlled is a metal rod made of stainless steel 0.48 m long and 0.3 m in diameter. The temperature profile along the rod length is controlled by manipulating 20 band heaters (with a maximum power of 50 W) wound around the rod. The temperature was measured at 11 points by copper-constant thermocouples. A personal computer was used for execution of the estimation and control algorithm with a sampling time of 20 s. The PI control law was used for the controller in the proposed control scheme. The values of physical parameters in the mathematical model, shown in Table 1, are quoted from the literature except for the heat transfer coefficients, which were determined by preliminary experiments.
7"~,i is an element of vector x~ which is the steady-state solution of Equation (10). Since the coefficient matrix H is a function of the actuators' allocation vector Za, the value of J from Equation (27) also depends on both Za and q. In this paper, the gradient method is used to find the minimum value of the objective function. The gradient vector of J can be formulated as shown below:
~...~
aJ -a- q- - ~
[~JaJ aJ ~lql" ' o~q2 ' ' ' ' ' aq j
(28)
,...,
'''''
aJ C~qr
(29)
Moving bed adsorber
The elements of these gradient vectors are derived by differentiation of Equation (27) with respect to z a and q as follows:
OJ = 2(x,T_xT)NAcIB ~ o~I -~zjq
0zj
Figure 4 shows the experimental apparatus of the moving bed adsorber, which consists of a stainless steel column of 1.2 m high, with outside diameter 0.03 m and inside diameter 0.025 m. Adsorbent particles were fed in
Table1
(30)
cgJ = 2(x,T_xT)NAEIBcH &l cgqj Oqj
(m2/s) ( W / ( m . K)) hD (W/(m 2 " K)) h~ (W/(m 2. K)) h: (W/(m 2 . K)) 7 ( s 1)
(31)
Power24Ch Controller]
D,I 2t
Thermocouple
i
4.44 × 10 6 16.3 40.3 239 232 1.3 × 10 3
GPIB
1,1 ,91
Fgil
I
[ t~gi~
Datalogger
Figure3
Parameters for the experimental process
Experimental apparatus of simple heat conduction process
GPIB
Co~~
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
a Adsorber
. Manometer
b Hopper
Digiufl-multimeter
c Rotary-valve
Personal computer
d Blower
363
:Thermocouple
e Feed Band heater f Preheater a Power controller
g Seal gas
m
00000 00000 00000
D
d
Figure 4
Experimental apparatus of moving-bed adsorber
to the top of the adsorber from a hopper and preheated air was introduced from the bottom. Activated carbon was used as an adsorbent and adsorption of the moisture in the air occurred in the adsorber. Although the thermal dynamics of the system was complicated, it was modelled by the following simple linear partial differential equation, and the control system was designed in accordance with the procedure mentioned above.
DT at
a O:T /: &2
fl aT 1 Oz
Results and discussion
y(T- T=)+ag+aq
(32)
where C c = {l -- e ) C p p p p
{(1
+
eCpgpg}S + CPoPoSo
e)Zp +ek.g}S+k.oS o
G¢ =
Ce Ce
y-
10 s. The cut off eigen-mode was set at 7. The thermal properties in the mathematical model were quoted from the literature, and heat transfer coefficients were obtained from the step response of the temperature profile in the adsorber. The values are shown in Table 2. The control experiments were carried out for the various numbers of observation points in order to test the state estimator performance.
hm/ Ce
Simple heat conduction process Figure 5 shows an example of controller performance, in the form of the transient states from a certain initial state to the target distribution indicated by a broken line. The measurements were carried out at : = 0, 0.5 and 1. Solid lines indicate the estimated temperature profiles. Fairly good estimation was achieved only with a few points of observation and the temperature distribution coincides with the target distribution after 2400 s. Table 2
Thermal properties and mass flow rate
1 Ce
Control of axial temperature distribution was performed by manipulating 20 100 W band heaters wound around the adsorber. A personal computer was used for the measurements and the regulation of the power controller, and the time interval for control was set at
Air Activated carbon Stainless steel
2
('p
p
If"
( W / ( m . K})
(J/(kg- KI)
(kg/nf)
/kg/s)
0.03 1.5 21.5
1017 75~ 461
109 ~50 7700
7.27 × 10 5.33 × 10 4
Void fraction t Heat transfer coefficient h (W/(m ~ - K))
0.35 2(1
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
364
1.0
140 '
'
110
' Stainle's steel t = 2400.
•
i
0.8
L
i
Error- e(t)Te(t)
L
~¢_. 0.6
Z = [0.1,0.9] T Z = [0.4,0.6]T
0
L0.4
L.LJ
1-
0.2
50
I
0
200
0 •
200
•
f
,
0_2
°
0.4
5
I'
0 6
z
Figure
i
~(o)r~(o)
t = 120~N~
~ so
i
(a)
400
•
O.
1
(b)
Controller performance for the simple heat conduction
process
IOOC
BOO
Time [s]
400. 300
[-3
600
~- I 0 0 ~
Oi 0
1.0
03 C3
I
I
I
I
0.2
0.4
I
I
I
0.6
I
0.8
1.0
zF] l I
/ \
,
/
i
/
i
/
I
Figure 7
Effect o f sensor location on the state estimator performance: (a) transition o f estimation error; (b) target profi]e
~ - i
250 ~o
r I
(-xJ
-~I
, /,,/4#~ I
0
i.
0.2
~
0.4
I
;; .4 /i
0.6
/
0.8
1.0
.~
150
1St SENSORPOSITION o:iniUal
Figure 6
point
o:opUmal point
Convergence to optimal sensor locations for two sensors
5O
Figure 6 shows the result of the optimal sensor allocation with two sensors searched for the transient state toward the stead-state profile as shown in Figure 7(b). The white circles in the figure indicate the initial positions, the black points are the optimal points, and the solid lines are the trajectories of searching. The result shows that the same optimal position can be found despite different initial positions. The dashed lines are contour plots of the objective function of Equation (19). Figure 7(a) indicates the trajectories of estimation error for various pairs of observation points during the transient state toward the steady-state profile as shown in Figure 7(b). In this case, the optimal sensor allowcation was obtained as z~ = [0.27, 0.73] x. Obviously, the fastest speed of estimation was obtained when the sensors were located at the optimal allocation. An optimal allocation problem was solved for the case that the temperature distribution of stainless steel plate 48 cm long was controlled to a desired profile by means of several plate-heaters 2.5 cm wide with a maximum generation rate of 65 W. Figure 8 shows the steady-state temperature distributions attained by various numbers of heaters located at the optimal locations. The dashed line indicates the target profile. It is found that the target profile can be realized with five actuators.
Moving bed adsorber The control result is shown in Figure 9, where the closed response of the temperature distribution is indicated from the initial state. In this case, the state estimation
0
~
0.0
i
I
0.2
I
I
I
0.4
I
0.6
I
I
0.8
|
1.0
zH Figure 8 Steady state temperature profile achieved by the actuators which were optimally located was carried out by using 11 temperature measurements at the marked positions in the figure (the positions were not at optimal allocation). In Figure 9, it can be seen that the temperature distribution was estimated well except for the initial period and controlled satisfactorily within an error of + 5°C after 2500 s. Figure 10 shows the heater-output distributions corresponding to Figure 9. Zero heater-output around the middle of the column is considered to be the effect of disturbance due to the heat of adsorption.
Conclusion A new controller design algorithm for parabolic DPS is proposed and applied to a simple heat conduction process and a moving bed adsorber. The results show that the state distribution can be estimated accurately by using only a few sensors, and controller performance is also satisfactory. Furthermore, methods for searching the optimal sensors' and actuators' allocation are presented, and it is demonstrated that the state estimation and control performances are improved by assigning sensors and actuators at the optimal locations. The proposed control scheme is basically similar to the conventional modal control structure 4.5. However, it
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
170~
V
V
V
V
aoL
V
V
V
V
Since S is a stable matrix, e(k) converges to zero as k increases. The objective function J of Equation (18) can be rewritten using the initial error e(0) as follows:
V °
t:25oo
J : e(0)Tp'e(0)
gOF
olo~
"
----0.4
0.6 z
P' = I + STp'S
o
0.8
100
. -
.
estimated distribution, o measurements) .
.
.
.
.
.
.
_j
Now, P, the solution of the Riccati equation Equation (17), can be rearranged in the form of a Lyapunov equation: p = sPSr+
~" so
0.0
(A.3)
.0
[-]
Figure 9 Controller performance for the moving-bed adsorber ( - - - , target distribution, - - - - ,
is a solution of the following Lyapunov
j t =os__ ~
0.2
(A.2)
• where P' equation:
10 0.0
365
Q + KRK v
(A.4)
P' is related to P as shown below by calculating the trace of the above equation.
0.2
0.4
0.6 z
0.8
1.0
[-]
traceP = trace{V' + SV'S T + S:V'S T2 + ...} = trace{V'(l + STS + STeS2 + ..3} = trace{V'P' }
(A.5)
Figure 10 Distributions of heater outputs during the control in Figure 9 where: is essentially different in the formulation of state-space representation, where it has several advantages compared with the conventional model control approach. The most distinctive one is that the state variables can be estimated directly from a few observations using a Kalman filter and the whole state distributions are obtained numerically from the data of state variables. Furthermore, boundary conditions are involved in the state equation and consequently this approach can be also extended to the controller design of the boundary control problem.
V'
= Q
+ KRK r
In this study, the covariance matrices are assumed to be diagonal matrices with constant coefficients such as Q = 0.1I, R = 0.1I. Therefore, V' can be assumed as the coefficient matrix. Consequently, P' in Equation (A.2) may be replaced by P, as the form of the objective function is not affected.
Appendix B: Derivation of Equation (23) Differentiating Equation (15) with respect to sensor location z/
References 1 Matsumoto, S. and Yoshida. M. Kagaku Kogaku Ronhunshu, 1987, 13, 145; English translation in Int. Chem. Eng. 1989, 29, 158 2 Yoshida, M., Matsumoto, S. et al. Kagaku Kogaku Ronbunshu, 1988, 14, 147; English translation in Int. Chem. Eng. 1990, 30, 142 3 Ozisik, M. N. 'Heat Conduction', John Wiley, New York, 1980 4 Lausterer, G K. and Eitelberg, E, in 'Distributed Parameter Control Systems', (Ed. S. G. Tzafestas), Pergamon Press, Oxford, 1981 5 Ray, W. H. and Lainitois, D. G. 'Distributed Parameter Systems Identification. Estimation and Control', Marcel Dekker, New York, 1976
3P = A 0P A - A 0z/ &j - AP
c3CT
3P C r [ R + C P C T ] I C P A 8:j [R + CPC T]
_ APCT 0[R + CPC T] i CPA APC T [ R + C P C
T] I OK?pA
c): I 3t' - A P C T [ R + C P C T] I C ~ A
Appendix A: Derivation of Equation (19) When a steady state Kalman filter is used as a state estimator, the dynamics of state estimation error e can be expressed as follows: e(k + 1) = Se(k)
I CPA
o3z/
(A.1)
(B.1) Since [R + CPCT]q[R + CPC T] = 1, the following relation can be derived.
Controller design for parabolic distributed parameter systems: M. Yoshida and S. Matsumoto
366
3[R + C P C v ]-:
Wf'X({*'z)X(;'z)dzl = ~ J0
O,
for i * j for i = j
N({,), 3[R + C P C v ] = - [ R + C P C v ]-:
• [R + CPC v ]q
Oz/
(C.4)
Then the third term o f the right-hand-side in Equation (C.2) can be simplified as:
ac
= - [R + C P C v ] -I oqz7 P C v [R + C P C s l-'
s,i i=l N(:,)
- [R + cPCT ]-Ic--~ cT [R +cPCT] -1
'
*=l
l
z) 2 dz '
~2
¢lzj
s,i
- [R + C P C T ] - 1 C P - ~ [ R
i=l
+ c P C T ] -1
OZ,
(B.2) R e a r r a n g i n g E q u a t i o n (B.1) with E q u a t i o n (B.2) and the gain matrix K in E q u a t i o n (16), E q u a t i o n (23) is
If eigenvalue modes greater than n are negligible, T~.~ and ~ become elements o f the state variable of x s and x* respectively. Therefore, E q u a t i o n (27) can be derived by substituting Equations (C.3), (C.5) into Equation (C,2) and by using these state variables.
derived.
Appendix C: Derivation of Equation (27)
Nomenclature
By using the finite Fourier transform, the steady state
Bi Cp
temperature Ts(z ) can be expressed as:
T~(z):~X(~i'z)~, i=l N ( { , )
(C.1) '
Substituting the above equation into E q u a t i o n (26), we have:
j:
i ~ T , ( z ) Z d z - 2 i i T , ( z ) ~ " X(~,,z) ~,.dz
i=1 N(~i) -
X(~i'z) ~ i X(¢,)
'
dz
'
(C.2)
The second term o f the right-hand side in E q u a t i o n (C.2) can be rewritten as follows:
,'=1
N(~i)
d g h l N P Q q S T T~ t u W w X x y z a t A p
Biot number , hl/)~ specific heat (J/kg . K) diameter (m) heat generation rate (W/m) heat transfer coefficient (W/mz • K) length (m) norm of eigenfunction covariance matrix of estimation error heat generation rate (W/m~) heat of adsorption (W/m) cross-sectional area (m:) temperature (K) surrounding temperature (K) time (s) vector of manipulated variable mass flow rate (kg/s) disturbance eigenfunction vector of state variable vector of output variable dimensionless distance thermal diffusivity (mZ/s) void fraction thermal conductivity (W/m • K) eigenvalue density (kg/m3)
Superscripts
'
Fourier transform
= 2~ ~ J o T * ( z ) X ( ~ i , z ) d z = 2
~,:, ~
5 (c.3)
The kernel function is an o r t h o g o n a l function expressed as follows:
Subscripts g i o p
gas phase mode of eigenvalue outer casing particle phase