Greek Letters a = Henry's law constant, molh. atm pp = density of catalyst particle, g/cm3 Literature Cited Bertrand. G.. Thomas, P., "Practical Biological Chemistry," Translated by H. A. Colwell, pp 61-67, G. Bells and Sons Ltd., London, 1920. Bizhanov, F. B., Sokolskii. D. V., Popov, N. i., Malkina. N. Ya., Khisametdinov. A. M.. Kinet. Kafal., 8 [3], 620 (1967). Brahme, P. H., Ph.D. Thesis, University of Bombay, India, 1972. Brian, P. L. T., Hales, H. B., A.LCh.E. J., 15, 419 (1969). Choudhary, V. R., Doraiswamy, L. K., Ind. Eng. Chem., Process Des. Dev.. 14,227 (1975). Furusawa, T.. Smith, J. M., A.6Ch.E. J.. 20, 86 (1974). Hofmann. H., Bill, W.. Chim. Ing. Techn., 31, 81 (1959). Kittrell, J. R.. Adv. Chem. Eng., 8 , 110-147 (1970).
Kolthoff, I. M., Belcher, R. V.. "Volumetric Analysis," Vol. 3, pp 475-488. Interscience, New York, N.Y., 1957. Law, V. J., Bailey, R. V., Chem. fng. Sci., 18, 169 (1963). Levenburg, K.. Ouart. Appl. Math., 2, 164 (1944). Mozingo, R., Org. Syn., 21, 15-17 (1941). (astergaard, K., Adv. Chem. fng., 7, 71-131 (1968). Rubin, D. I., Chem. Eng. Prog. Symp. Ser., 59 [42], 90 (1963). Ruether, J. A., Puri, P. S.,Can. J. Chem. Eng., 51, 345 (1973). Satterfield, C. N.. "Mass Transfer in Heterogeneous Catalysis," pp 92-97, 107-1 16, MIT Press, Cambridge, 1970. Schnyder. W. A,, Ph.D. Dissertation, Polytechnic Institute of Brooklyn, 1962. Schoenemann, K., Gen. Chim., 77 [4], 89 (1957). Smith, J. M., "Chemical Engineering Kinetics." pp 384-386. McGraw-Hill Inc., New York, N.Y., 1970. Wiebe, R.. Gaddy, V. S..Heins, C., lnd. Eng. Chem.. 24, 823 (1932).
Receiued for reuiew, March 5,1975 Accepted August 11,1975
Linear Direct Digital Control Algorithms for a Class of Distributed Processes Rajakkannu Mutharasan' and Rein Luus* Department of Chemical Engineering, University of Toronto, Toronto, Ontario, Canada M5S 1A4
Two procedures for deriving linear direct digital control algorithms for a class of hyperbolic distributed parameter systems are presented. Control of a flow-forced heat exchanger is chosen as an example to illustrate the procedures. The first method involves linearizing the state equations and deriving the process transfer function. The resulting algorithm responds very well for small loads. The second method is based on directly searching for the optimal algorithm parameters. The search procedure requires only nominal computer time. This procedure of design, unlike the first method, can accommodate large load and set point changes. Due to the nonlinearities of the system equations, the optimal controller parameter values vary for different load conditions. However, this method enables one to retain the linearity of the algorithm.
With the increase in the use of digital computers for controlling process systems, study of sampled-data control systems and design of direct digital control algorithms have presented interesting and challenging problems. Several authors have presented direct digital algorithms for lumped-parameter systems (Mosler et al., 1967; Moore et al., 1970; Luyben, 1973). Recently Mutharasan and Coughanowr (1974) presented nonlinear digital control algorithms for a class of hyperbolic distributed-parameter systems. Since it is often convenient to use linear DDC algorithms so that existing DDC programs need not be changed, a study was undertaken to develop linear DDC algorithms for the same class of hyperbolic systems. The class of hyperbolic systems considered in this investigation is described by the partial differential equation
where P l ( t ) and P2(t) can be load or manipulated variables, c is the state variable, and t and x are dimensionless time and distance. A typical process that belongs to the above class is a flow-forced shell-and-tube heat exchanger which is described by ac
ac
-+ m(t)-
= -Po[m(t)Ibc at ax Department of Chemical Engineering, Drexel University, Philadelphia, Pa. 19104
where
t = -00 L x = -Y
L
0
m=l+r
U
This model has been previously used in the literature in the study of continuous control of flow-forced exchangers (Koppel et al., 1970). For this process, the manipulated variable, m ( t ) , is the velocity of the heated stream. Also, the velocity dependence of the heat transfer coefficient is included in the model. The initial conditions of the process are given by c(x, 0 ) = cl(0) exp(-Pox)
(3)
m(0) = 1
The load variable is the inlet state which is given by c(0, t ) =
Cl(t)
Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 1, 1976
(4) 137
For t < 0, the inlet temperature, cl(t) is unity and the process is at steady state. At t = 0, cl(t) deviates from unity and hence causes the outlet temperature to deviate from the set point value of the loop, exp(-Po). The flow through the exchanger is then varied (and thus m ( t ) )according to a control algorithm (which is to be designed) so that the outlet temperature returns to the set point value. In an experimental closed-loop control of the exchanger process described here, Paraskos and McAvoy (1970) show that conventional Proportional-Integral controller with Ziegler-Nichols feedback settings gives poor control. Koppel et al. (1970) have used a two-point continuous controller on an experimental flow-forced exchanger. In this investigation we will consider sampled-data control algorithms for the flow-forced heat exchanger. The method adopted in this paper enables one to retain the distributed nature of the process. Since the control system is sampled-data in nature, the output of the process is measured every sampling period t , and the manipulated variable is updated every sampling period. Between sampling periods the manipulated variable is kept constant. In general the linear DDC algorithm can be expressed as
rn(kt,) =
r=O
g i e ( ( k - i)t,)
- i5 hirn((iz- i - l)t,) =O
+-, d o 5)
I
Figure 1. The control system. Ed(1, 2 )
(11)
so that
u(z)= e-Po2-n-
2
2 - 1
(5)
From the block diagram
Minimal Algorithms
where
The process model in eq 2 can be linearized about the initial steady-state to yield
+ Poe (1- b)&
+ 0 Z-(n+l) + . . .
U ( z )= Z(e-POe-SE(O,s)}
HG,(s) =
+
= e-Poz-n
where n is an integer and nt, = 1. This is by no means a restriction of the method, but a convenient case for illustration. From Figure 1,it is seen that
where gi and hi are constants and e is error. The advantage of the linear DDC algorithm in implementation is that only one subroutine is needed for calculating the manipulated variable value for all the linear loops in the DDC system. This feature is important when one is restricted by the available active memory for the DDC system.
aE a? - - = -Po? a t ax
-Po-s
(6)
where the tilde refers to the deviation state and the bar to a steady-state value. Taking the Laplace transform of eq 6 gives
By using the boundary condition
(-
S
a = Pocl(O)(b
- l)e-Po
and nt, = 1 have been introduced. The z transform of eq 15 yields
The pulse-transfer function D ( z ) can now be derived using the expression (Tou, 1959)
in the form
we can integrate eq 7 to give (19) which can be simplified to 2.1(s ) e-(s+
Po)x
(9)
Thus the outlet state can be related to the inlet state and the manipulated variable
/- Ll D(z)=
at,]
(1 - 2 - n )
(20)
The transfer function in eq 20 can be written in the form of eq 5 as This enables one to obtain a conventional control system block diagram as shown in Figure 1. To derive a minimal prototype load algorithm, let the inlet state be subjected to unit step change. We will require the outlet state to settle down to the set-point value in one sampling period. The desired response of E(1, t ) can be expressed in the 2 domain as 138
Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 1, 1976
Application of a similar technique to give set-point minimal prototype response yields the same pulse transfer function (eq 20). This is rather a unique feature of the present class of distributed parameter systems. For most lumped parameter systems, the load minimal and set-point minimal prototype algorithms differ considerably.
I
I
,
I
I
r
7
I
-- 042t
..-.
--.
2041
0
0
5a 0 3 9 0 5 0385l 03800375-
I
I
1
I
i
I
040-
i 0 393 038-
h
2 0.37-
1
2
3
4
5
6
1
0 36-
I
0.33;
i
1
0
I
I
I
7
8
---I
0 35034I
1
2
,
I
3
4
5
TIME,
,
I
6
7
1
,
t
t Figure 2. Response of the minimal algorithm (eq 14) to a load change: E(0, t ) = 0.05u(t);t , = 0.5; Po = 1.0; b = 0.6.
Figure 4. Response of the minimal algorithm (eq 14) to a load change: E(0, t ) = 0.15u(t);t , = 0.5; PO= 1.0; b = 0.6.
14
where the subscript s denotes the set-point value and t d is the time required for the integrand to reach zero. The random search procedure of Luus and Jaakola (1973) was used to obtain optimal algorithm parameters. The calculation procedure may be summarized as follows. 1. Take initial values for the algorithm parameters
TIME,
-v
q 0 4 3 -
i
042-
5 041
gi(O),i = 0 , 1,.. . , r
h p ,i
040 039
-1
2. Choose an initial range Ri over which search is to be carried out. 3. Choose sets of random numbers in the interval [-0.5, 0.51, say 100 sets; denote these sets of random numbers by
0 38
0 37
0
= 0,1,. . . , P
3
5
0
1
2
3
4
TIME,
5
6
7
9ki.
8
t
Figure 3. Response of the minimal algorithm (eq 14) to a set point change: E A l , t ) = O.lu(t);t , = 0.5; PO= 1.0; b = 0.6.
The process was simulated with the algorithm (eq 21) and was found to respond favorably only to very small load changes. T o allow large load changes, the requirement of the response of the system was relaxed; that is, the outlet state was now required to settle within two sampling periods after the load is sensed. Thus
zd(l,z) = e-Po2-n + ye-Po,-(n+l) + 0 2 - ( n + 2 ) + . . .
(22)
where y is a coefficient between 0 and 1. The resulting pulse transfer function takes the form D(z)= [(l- 7) (27 - 1)2-1 - 7 2 - 2 1
(- $)
+
(23) The system was simulated (with y = 0.5) and sample load response and set-point response are shown in Figures 2 and 3. In Figure 4 the response of the system for a load of 15% is shown. The algorithm (eq 23) makes heavy compensation initially and therefore the response is oscillatory. Since the algorithm is derived from the linearized state equation, the response of the system is different from what was imposed (eq 22). However, linearization does lead to a workable algorithm which responds favorably for small loads. An alternate method of designing linear DDC algorithms is to directly search over the parameters gi's and hi's in eq 5. We will define the performance index as the integral of error; namely
4. Use the first set of random numbers [ K = 11 and calculate algorithm parameters by t k i .
gi = gicO)+
i = 0 , 1, . . . , r
hi = hi@) + VkiRi, i = r
+ 1,r + 2 , . . . r + p - r
h,= 1 -
c hi
P-1
i=O
5. Integrate the system equation (eq 2) for a given load change by the method of characteristics and evaluate the performance index given by eq 24. 6. Repeat steps 4 and 5 for other sets of random numbers. 7. Choose the set of algorithm parameters which gives the smallest value for the performance index and denote this set gi(O)and hi(o). 8. Reduce the size of each region rij by a small amount, say 5%. 9. Go back to step 3 and repeat calculations until a specified number of iterations have been done. In the above calculation procedure, we choose the algorithm parameters to satisfy the constraint 1+ghi=O i=O
which guarantees zero steady-state error for a stable loop. In this investigation, the search was limited to two parameter algorithms, that is r = l and p = 0 in eq 5. The parameter ho is -1 by the requirement that steady-state error should be zero. Also, the parameters were searched in the region between -10 and +lo. Tables I and I1 summarized the results for the sampling periods of 0.5 and 1.0. Although each particular parameter values gives optimal response for the particular load, the response is often sluggish or unstable for other load conditions. For example, Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 1, 1976
139
Table I. Optimal Algorithm Parameters for Different Load Conditions: t, = 0 . 5 ; =~ 1.0; ~ b = 0.6 Load
0.3 5.4874 -3.6290 3.742 X
go g,
P.I.
0.2 5.3958 -3.07617 1.4992 X
0.1 7.0131 -3.8664 3.6692 X
-0.2 9.7512 -1.1546 7.9351
-0.1
9.9634 -1.8283 2.0365
X
Table 11. Optimal Algorithm Parameters for Different Load Conditions: t, = 1 . 0 ; ~=" 1.0; b
=
-0.3 9.9919 6.4361 2.0152
X
X
0.6
~~
Load
0.3 5.2348 -1.97 75 5.0018 x
go g,
P.I.
0.2 7.5318 -5.07 21 2.1337
0.1
7.6644 -2.6133 3.9183 X
X
-0.2 9.6696 0.2996 1.1115 X
-0.1 9.9849 0.66721 2.6140 X
-0.3 9.9919 6.43611 3.242 X lo-'
--I
0 50 0 48
-+ -- .
0 46
0 44
0 042
k-0 4 0 3
0.38
_----
2036 0 34
0 32
02't
0 30 0
1
2
3
4 TIME,
5
6
7
8
9
f Ltd( ~ ~ (t 1) -, &Vi,t ) I 2dt
(25)
i= 1
where c ( ~(1,t) ) is the response of the system for load condition i. In this investigation six load conditions were chosen, namely step changes of f0.1, f0.2, and f0.3. For the system parameter PO= 1.0, b = 0.6, t , = 0.5,the optimal algorithm parameter values obtained by direct search procedure described above were go = 5.4600 and gl = -3.0018. When a sampling period of 1.0 was used instead of 0.5, these parameter values were found to be go = 5.3686 and gl = -1.3478. Load response of the system with the above algorithm parameter values for different load conditions is shown in Figures 5 and 6.
Conclusions Two approaches for deriving linear direct digital control algorithms for a class of distributed parameter systems were developed. The first method is based on the linearized state equation while the second method allows one to retain the original nonlinear system equation. The pulse transfer function of the digital controller is easily derived for the first method, since the procedure is analytic. The second method involves numerically searching for the optimal linear controller parameter values with respect to the integral of squared error criterion. Use of the stochastic optimization procedure of Luus and Jaakola (1973) makes this approach simple and requires very nominal computational effort. Since the first procedure relies on the linearized equations, the algorithm responds very well for only small load and set point changes. For large changes the response is oscillatory as may be expected. Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 1, 1976
2
3
4 TIME,
the setting for a load of -0.3 from Table I gives unstable response when the load is +0.05. This is highly undesirable. Therefore, the performance index was modified as follows
140
1
t
Figure 5. Response of the two-parameter linear algorithm to a load change: Z(0, t ) = 0.3u(t): t , = 0.5; PO= 1.0; b = 0.6; go = 5.46; gl = -3.0018.
P.I. =
0
v 5
6
7
8
9
t
Figure 6. Response of the two-parameter linear algorithm to a load change: Z(0, t ) = -0.3u(t); t. = 0.5; PO = 1.0; b = 0.6: go = 5.46: gl = -3.0018.
The second procedure, unlike the first one, can accommodate large load and set point changes. Due to the nonlinearities of the system equations the optimal parameter values vary considerably for different load conditions. However, this method enables one to retain the linearity of the algorithm and to design optimal linear algorithm for a set of given load changes. The two methods can easily be applied to other systems which can be represented by eq 1. The direct search procedure can be used to obtain optimal parameters for higher order algorithms.
Nomenclature a = a constant defined by eq 16
A , = flow-area of tube-side fluid A , = heat transfer area based on the shell side b = coefficient relating the variation of heat transfer coefficient to velocity c = dimensionless state variable cp = specific heat of tube-side fluid c , = set point value e = error gi = constants in the linear algorithm (eq 5 ) hi = constants in the linear algorithm (eq 5 ) L = length of the heat exchanger m ( t) = manipulated variable P l ( t ) ,P2(t) = load or manipulated variable PO = process parameter; see eq 2 Ri = regionof search t = dimensionless time t , = sampling period T = temperature of tube-side fluid Ti, = initial steady-state inlet temperature of the tubeside fluid T, = wall temperature of the heat exchanger u ( t ) = unit step function UO = overall heat transfer coefficient based on the outside area, a t the velocity D
0 =
initial steady-state velocity of tube-side fluid
ir = deviation in velocity from steady state x = dimensionless space variable y = distance along the heat exchanger z = 2 transform variable Greek Letters y = a constant used to specify the value of outlet state a t
the second sampling instant after the load is sensed (see eq 22) p = density of tube-side fluid v k i = random numbers in the interval [-0.5,0.5] from the set k 0 = time Superscripts tilde = refers to deviation state bar = refers to steady state
Literature Cited Koppel, L. E.,Kamman, D. T., Woodward, J. L., I d . Eng. Chem., Fundam., 9, 198 (1970). Luus, R., Jaakola, T. H. I., AEhEJ., 19, 760 (1973). Luyben, W. L., “Process Modelling. Simulation, and Control for Chemical Engineers”, p 520, McGraw-Hill, New York, N.Y.. 1973. Moore. C. F., Smith, C. L., Murill, P. W., lnst. Contr. Systems, 43, 70 (1970). Mosler. H. A,, Koppel, L. E.,Coughanowr, D. R.. AEhEJ., 13, 768 (1967). Mutharasan, R., Coughanowr. D. R.. Ind. Eng. Chem.. Process Des. Dev., 13, 168 (1974). Paraskos, J. A,, McAvoy, T. J., AlChE J., 18, 754 (1970). Tou, J. T., “Digital and SampledData Control Systems”, p 501, McGraw-Hill. New York, N.Y., 1959.
Received for review March 6, 1975 Accepted August 18, 1975 This work was performed with the asssistance of a grant from the National Research Council of Canada, Grant No. A-3515. Computations were carried out with the facilities of the University of Toronto Computer Center.
Sampled-Data Proportional Control of a Flow-Forced Tubular Reactor Rajakkannu Mutharasan’ and Donald R. Coughanowr Deparfment of Chemical Engineering, Drexel University, Philadelphia. Pennsylvania 19 104
Though continuous proportional feedback control of an isothermal tubular reactor is always stable, sampleddata proportional control can make the system output to oscillate indefinitely. Through analysis and simulations it is shown that for any feedback gain, there is a magnitude of step change in inlet concentration, which can cause the system output to oscillate indefinitely.
Introduction and Previous Literature Several published papers have appeared on continuous feedback control of distributed parameter systems. Koppel (1966) considered continuous nonlinear feedback control of tubular chemical reactors and heat exchangers. Koppel et al. (1970) reported theoretical and experimental results on two-point linear control of a flow-forced heat exchanger and extended the principle to other parametrically forced distributed parameter systems. Paraskos et al. (1970) studied experimental feedforward computer control of a flowforced heat exchanger and reported that their algorithm was far superior to conventional Ziegler-Nichols settings. Implementation of their feedforward algorithm involved solution of finite difference equations, which is rather elaborate for an industrial application. Seinfeld et al. (1970) compared continuous proportional feedback, feedforward, and optimal control of a flow-forced isothermal tubular reactor. Their paper contains useful results on offset and stability of the system under proportional control. They reported that the system is stable, irrespective of the value of the proportional gain. Oscillations in outlet concentration increased as the proportional gain is increased; however, there is an upper limit on the gain because of the physical requirement that the velocity should be greater than zero. Several continuous optimal control algorithms, for open loop and closed loop, configurations have been reported for the flow-forced reactor-heat exchanger system (Vermeychuk et al., 1973; Seinfeld et al., 1968, 1970). Review of the existing literature indicates that very few authors have worked on feedback sampled-data control of
distributed-parameter systems (Palas, 1970; Hassan et al., 1970). This topic is of great importance in an industrial environment because of an exponential increase in the use of computers as a control element. Several authors have developed methods of design of direct digital control algorithms for lumped-parameter systems (Mosler et al., 1967; Moore et al., 1970; Cox et al., 1966; Dahlin, 1968), but similar work for distributed systems is lacking. This investigation helps fill part of the gap that exists in the literature. In this paper we present the results of a study on the nature of feedback sampled-data proportional control of a flow-forced tubular reactor. Ultimate loop gain of a sampled-data control of a first-order system is lower than that of a continuous control system (Mosler et al., 1967). One conclusion we may draw is that sampling operations cause the lumped parameter system to be less stable, but a stable and responsive proportional controller can still be designed. In the present distributed parameter system, we show that sampling the output with manipulation proportional to error yields a system which sustains oscillations indefinitely (at least theoretically), thus leading to the conclusion that proportional sampled-data control is unsatisfactory for the flow-forced tubular reactor. T h e Process The process under consideration is a flow-forced tubular reactor in which an nth-order isothermal reaction, A B, is taking place. Under the assumptions of (1) plug flow of fluid, (2) uniform physical properties, and (3) negligible axial dispersion, the dynamic behavior of the tubular reac-
-
Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 1, 1976
141
