OPTIMAL CONTROL OF A CLASS OF DISTRIBUTED-PARAMETER SYSTEMS WITH DISTRIBUTED CONTROLS LOWELL B. KOPPEL AND YEN-PING SHIH1 School of Chemical Engineering, Purdue Uniuersip, Lafayette, Ind. 4790 7
A strong minimum principle is derived for the distributed-parameter system a x / & u ( r , t ) ] with performance index P(u) = f d f
I I F [ x ( r , f ) , u(r,t)]dr.
control theory for lumped-parameter processes is This theory has resulted in analytical solutions for optimal control of linear systems subject to performance criteria such as minimum integral-square-error, minimum time, and minimum fuel ('4thans and Falb, 1966). For other lumped-parameter systems and other performance criteria, efficient computational algorithms have been devised to obtain numerical solutions. For distributed-parameter systems, comparable results are not yet available. In the present paper, we show that an important class of distributedparameter problems, in optimal control of processes of chemical engineering interest, can be solved by using the already available results for lumped-parameter problems. We first examine a particular problem in optimal control of distributed-parameter systems. This problem offers minimal computational resistance, and in fact permits a complete, analytical solution, which enables mathematical verification of previously used conjectures on feedback control, and suggests a technique for optimizing a general class of distributed-parameter control systems by using results obtained from lumpedparameter theory. A strong minimum principle is ultimately obtained by generalizing this technique. For the particular problem, the partial differential equation describing the system dynamics is PTIMAL
and the initial and boundary conditions are
= f[x(r,t),
Tubular reactors and heat exchangers
are described by these system equations. Optimal controls are found eter problem; some analytical solutions are presented.
0 relatively advanced.
+ v(f) a x / &
by
solving a related lumped-param-
heat flux profile. Both x ( r , t ) and 4 7 , t ) are normalized. I t is assumed that the wall flux can be varied with both position and time. Equation 2 states that the initial temperature distribution is arbitrary, while Equation 3 states that the temperature of the fluid entering the exchanger is constant a t the desired final value. As a performance index, choose the following:
P(u)
=
1"1'
[p2x2(r,t )
+ u2(r,t ) ]dr dt
(4)
where p is a constant weighting factor. Necessary Conditions
Shih (1967) has shown that the necessary condition for an unconstrained u ( r , t ) to provide the minimum of P(u) is that it yields a stationary point of the Hamiltonian
(5) where the adjoint variable p ( r , t ) satisfies
p(1, t ) = 0
(7)
p ( r , t,) = 0
(8)
The minimizing u is found by differentiation of H to be
u(r, t ) = -p(r, t )
(9)
where we have assumed that u is not constrained.
A physical system which is described by such equations (after appropriate normalization) is the tubular, plug-flow heat exchanger with constant properties, constant velocity, no axial mixing of energy or mass, perfect radial mixing, and no wall capacitance (Koppei et al., 1968). Here t represents normalized time and r represents normalized axial distance from the exchanger entrance. Thus, t = 1 a t one residence time, and r = 1 a t the exchanger exit. The state variable x ( r , t ) represents deviation of the fluid temperature from the desired final temperature profile, and the control variable u ( r , t ) represents deviation of the wall heat flux from the final Present address, Cheng Kung University, Tainan, Taiwan, China. 414
I&EC FUNDAMENTALS
Solution
Because of the particularly simple form of the system and adjoint equations (Equations 1 and 6), they are easily solved by the two-dimensional Laplace transform method (Voeker and Doetsch, 1950). The solution is: xo(r
cosh [p(l - r ) ] [p(1 - r t)]
- t ) cosh
+
(r, t ) eR1
p(r, t ) =
f p x ( r , t ) tan11 [ p ( ~- r ) ]
(Y, t ) E R ~
4 ~ . L ( Y , t ) tanh
(I,
lo
= -
[p(t, - t ) ]
t ) tR2
(7,t ) U(Y,
Feedback Control
(11)
eR3
t)
where R I , R2, and R3 are the regions in r and t space shown o n 1, region Rz disappears; only R1 and R3 Figure 1. If t i in Equations 10 and 11 apply.
>
Illustrative Responses and Controls
Figures 2 and 3 show temperature piofiles and control profiles for t , = 1 and p 2 = 5. The optimal control reduces the temperature much more quickly than does the simple step control u = 0, for lvhich the response is shown in dashed lines. Figure 4 presents the optimal exit temperature x(1, t ) for t, = 1 and various values of p2. The case for p2 = 0 corresponds to the step control u = 0. As expected, higher penalties cause faster response, but require more control. Figures 5 to 7 present comparable results for t, = 0.5.
I n a previous publication on optimal control of linear, distributed-parameter systems, with a quadratic performance index (Koppel et al., 1968) a linear relation was conjectured between the optimal state and adjoint variables. In the present case, the conjecture takes the form d r , t ) = K(y, t)x(r, 0
(12)
This conjecture is applied to solve the system and adjoint equations because an analytical solution such as Equations 10 and 11 is generally not available. We follow the typical procedure by differentiating Equation 12 with respect to Y and t , substituting the results into the adjoint equation (Equation 6), and combining the result with the system equation (Equation 1) and Equation 9, to obtain
I.o
0.8
I 0.6
4
P(r,tl -ucr,t)
0.4
t
0.i 0 0
I
r
0
Figure 1. Characteristic lines for wall flux control problem
0
44
0.2
Figure 3.
r
Optimal control for
1.o
0.8
0.6
ti
3
1
1.0
08
0.6
XtlJ) 0.4
0.2
0
0.2 Figure 2.
0.4
r
0.6
Optimal profiles for t f u(r,t) = 0
OB
1.0
0
>1
t Figure
4.
Optimal exit temperature for t , VOL. 7
3
1
NO. 3 A U G U S T 1 9 6 8
415
0.2
0 Figure 5. Optimal profiles for
---
tj
= 0.5
0.4
Figure 7. Optimal exit temperature fort, = 0.5
u(r,t) = 0
1.0
Comparison of Equations 9, 11, and 17 shows that Equation 18 provides the optimal feedback control. For the present simple case this verifies a feedback conjecture similar to that used in previous work (Koppel et al., 1968). However, in that previous work, the control was a function of time only, and the conjecture was of the form
0.8
0.6
u(t) =
-p(r$0.4
sb'
K(r, t ) x ( r , t ) d r
An analytical solution is not available for the case with u a
UKt)
function of time only. The results of the present paper cannot be used as a proof of the above feedback relation, which must still be regarded as a conjecture. However, the fact that a linear feedback conjecture can be mathematically verified for the case with u a function of both position and time is encouraging.
0.2
Comparison with Dynamic Programming
0
0.4
0.2
0
Figure 6.
r
0.6
0.8
I .o
Optimal control for t , = 0.5
Sirazetdinov (1965) presents a technique for dynamic programming solution of the problem we have solved above. A modified version with finite tf is presented below. Define
Since xo(r) is arbitrary, this can be satisfied only if K ( r , t ) satisfies
dK
at
+ dK - K 2 + p2 = 0
Also, Equations 12, 7, and 8 require the following conditions on K: K(1, t ) = 0 (1 5)
K(r, t j ) = 0
(16)
The solution to Equations 14 to 16 is found by straightforward application of the method of characteristics:
K(r, t )
=
i
lr
tanh tanh
[lr
(1
G
(tj
- r)1
- t)]
(r, t ) E R I
( r , t ) ERZ,RS
(17)
Equations 12 and 9 combine to the feedback control u(r, t ) = - K ( r , t ) x ( r , t ) 416
I&EC FUNDAMENTALS
Then the optimal control must satisfy the following form of the Hamilton-Jacobi equation:
(18)
We conjecture the form of V to be
Then,
Using Equations 1 and 3, integrating by parts the term involving d r l b r , imposing the condition K(1, t ) = 0, and combining ivith the expression for T V in Equation 19, yield
dV dt
I
I
i
t
bK
bh'
Kxu
'I
+ -2 u2
dr
(23)
Minimizing with respect to u requires u = -Kx and then imposing the condition that the minimum must be zero yields Equation 14. Thus, the analytical solution obtained in Equations 10 and 11 verifies the conjectured form of the minimum performance functional in Equation 21. Finite-Heater Approximation
For tf
at
K,(r) = p tanh [k(1 - r ) ]
Consider the system with dynamics
(24)
Therefore, a physically implementable approximation to the optimal control is suggested in Figure 8. The wall flux is provided by a sequence of discrete heaters, each operated by the feedback scheme shown. The temperature at the exit of each section is used to regulate the heater over that particular section. Such approximations are necessary because the exact feedback control in Equation 18 requires infinitely many sensors and controllers. Shih (1967) has shown for the case with u a function of time only that a four-sensor approximation provides a response essentially identical to the optimal response. Similar results can be anticipated for finite approximations to Equation 18 such as that suggested in Figure 8.
For optimal control of the lumped-parameter system
dt
= u(t)
x ( 0 ) = xo with performance index
it is well known (Athans and Falb, 1966) that the optimal feedback control is u(t) =
-K(t)x(t)
dX
-
at
+ v ( t ) -dr dX
= f[x(r,
where x is a column vector of n state variables, u a column vector of m control variables, f a column vector of n functions, and v ( t ) 2 0 a given scalar function. A physical example of such a system is control of a tubular chemical reactor. Typical normalized mass and energy balances may take the form
dc -
at
+ v ( t ) bcbr = --R(c,
p
tanh
[p(tf
- t)]
T)
-
(27)
x(r, 0 ) =
(28)
Explicit expressions for the optimal response and control are
The similarity between these results and those obtained in Equations 10 and 11 for the analogous distributed-parameter system leads to the generalization presented below. This leads to a strong minimum principle for a n important class of distributed-parameter processes.
(32)
is specified for all r such that - h ( t f ) h(t) =
with gain
K(t) =
t ) , u(r, t ) ]
where c ( r , t ) is the dimensionless local concentration of the reacting species, u is the dimensionless velocity of the fluid in the exchanger, T(r,t ) is the dimensionless local temperature, R is the dimensionless local rate of decomposition expressed as a function of concentration and temperature, k is the dimensionless heat released per unit amount of reaction, and q is the dimensionless heat flux applied a t the reactor wall to control the concentration and temperature. I n general, there may be many reacting species, and a mass balance must be written for each one, thereby increasing the number, n, of components in the state vector. We assume that
Comparison with Lumped-Parameter System
-
Finite heater configuration
Generalization
3 1, the optimal gain is stationary
dX
Figure 8.
6
r
6 1, where
s,'
u(7)dt
(33)
but that u(r, t ) is defined only on 0 6 r 6 1. Physically, this means that the process lies on 0 6 r 6 1, but that some of the fluid which will enter the process during the operating interval 0 6 t 6 tf lies a t negative r a t time zero. I n other words, Equation 32 is a specification of the initial condition on the fluid initially in the process, as well as on the fluid which, though not initially in the process, does enter the process during the operating interval. I n the sequel, we for convenience generally regard u ( t ) as the velocity of an incompressible fluid moving in the r-direction through the tubular process, and X ( T , t ) as the state of the moving fluid, although Equation 30 and the results obtained are not restricted to such physical interpretation. Then, h ( t ) is the total distance traversed by the fluid a t time t, since the control process was started. VOL. 7
NO. 3 A U G U S T 1 9 6 8
417
where F is a scalar function. By virtue of Eqdation 35, Pl(u; r,) can be minimized by solving a lumped-parameter problem only. Kow define a performance index for the overall distributed-parameter process as
P(u) =
J
1
pl(u; ro)dro
(38)
-h(tO
f
I
Since for each ro, Pl(u; ro) involves values of x and u not involved in P1 for any other ro, and since each Pl(u; r o ) is minimized, P(u) will necessarily be minimized because it is merely the sum of all the independent values of Pl(u; ro). The performance index of Equation 38 may be written in a more recognizable form by interchanging the order of integration over r and t. As ro is varied, the appropriate values of to(ro) and t,(ro) change. The lowest value of T, of interest is -h(t,), since fluid initially a t lower values of r, does not reach the process entrance (r = 0) during the control interval 0 6 t 6 t,. Figure 9 shows the critical characteristic lines, plots of time us. position for fluid particles as they travel through the process. Points on the base line, t = 0, correspond to the various values of ro of interest. Figure 9, a, shows the case for h(t,) 1, while Figure 9, 6, shows the case for h(t,) 6 1. These diagrams make computation of t,(r,) and t f ( r o ) more convenient. For example, for 1 - h(t,) 6 yo 6 1, these figures show that a particle initially on the base line reaches r = 1 before t,, and further that the time required to reach r = 1 is given by h-I(l - ro). Continuing in this manner we obtain [ & ( i ois ) defined to be zero for particles initially in the process and t,(r,) is defined to be t , for particles finally in the process] :
>
0 I-h(t$
r
I
(b) h(tJ)4 I Figure
9.
Characteristic lines
The path
r ( t ) = r,
+ h(t)
(39)
is a characteristic of Equation 30. I t is the positional behavior of a fluid particle, initially a t ro, as a function of time. Consider the state variable behavior seen by an observer moving with the fluid particle, X ( Y , h, t ) . Then, using the chain rule
+
(34) Therefore, Equations 30 and 34 show that, for such an observer, the dynamics are dx -
dt
(TO
+ h, t )
= fIx(r0
+ h, t ) ,
+ h, t ) l
~ ( r o
These relations apply to either case shown in Figure 9. We next plot to(ro) and t,(ro) as functions of r, for the two cases shown in Figure 9. The results are shown in Figure 10. The integration in Equation 38 is a double integral over ro and t , with integration over t first, and then ro (Figure 10). This can be decomposed into three integrals over the separate regions where different expressions apply for t,(ro) and t,(ro), the limits on the inner integral. Thus, for h(t,) 3 1
(35)
which is equivalent to the lumped-parameter system dx
-
dt
u [ro = f ( x , u)
>
provided t to(r0),where to(ro) is the time a t which a fluid particle, initially at ro, first enters the process-Le., reaches r = 0. For ro 0, to(ro) is defined to be zero. Similarly, t f ( r o ) is defined as the time at which a fluid particle, initially a t ro, first leaves the process-Le., reaches r = 1. For ro 6 1 - h(t,), t,(ro) is defined to be t,, the time at which the control process is terminated; such particles do not reach r = 1 during the control process. Now, consider the fluid initially a t ro, and define for this fluid the performance index to be minimized as:
+ h ( t ) , t l ) dt dr,
(40)
and primes on R,, R,, and R, convert this to the result for h(t,) 6 1. I n either case, Figure 10 shows that when the order of integration is reversed, the result is
>
+ h ( t ) , t ] , ~ [ r +, h ( t ) , t l ] d t 418
I&EC FUNDAMENTALS
(37)
Finally, changing the variable of integration to r = ro yields
P(u) = c d t
J1 F[x(r,
t ) , u(r, t ) ] d r
+ h(t) (41)
the usual form of performance criterion for the distributedparameter system of Equation 30.
Case
11. h ( t f ) 6 1
Minimum Principle
Thus, the following strong result is obtained : Let x*(t - t o ; tf - to, x,) and u*(t - t o ; tf - to, x,) be the transient and control which provide the absolute minimum of
Pl(u) =
i?[x(t),
u(t)ldt
for the system t
with
x(to) = x,
0
Then, the transient and control which provide the absolute minimum of
..l o 4 b) h(t+)4 I Figure 10.
P(u) = 1 " d t
Extrance and exit times
u(t)l dt
F[x(r, t ) , U(Y, t)]dr
for the system
Thus, the computational procedure is to solve the lumpedparameter problem involving minimization of
PI(u)= J:F[x(t),
1'
with
(42) are as follows:
for the system of Equation 36 for every initial value x, which appears in Equation 32. Computational methods available for lumped-parameter systems may be used. This will yield a n optimal transient and control
x* = x*(t
- to;
t,
Case
I. h(t,)
1
X(Y, t ) =
- to, x,) (43)
u* = u*(t
- to;
tf
-
to, x,)
-
with ( t , to) and x,, the total control time and initial state vector, appearing as parameters. T o construct the solution to the distributed-parameter problem, a t each point (Y, t ) we calculate the pertinent initial position Yo
=
Y
- h(t)
U(Y,
t) =
(44)
Then, the initial value xo in x* and u* should be replaced by X,[Y h ( t ) ] which is obtained from the specification of the initial condition in Equation 32. The initial and final times in x* and u* should be replaced by values obtained from Figure 10. For each point (r, t ) , the appropriate replacement values are : Case I. h(t,) 1
Case
11.
-
X(Y, t ) =
>
(y,
(T,
t)eRb, Rc t)&
U(Y,
t) =
(45) VOL. 7
NO. 3 A U G U S T 1 9 6 8
419
+
An extension of this result follows immediately: Let x*(t; to, t,, y o , x,) and u * ( t ; to, t,, y o , x,) be the transient and control which provide the absolute minimum of
PI(^ = f h x ( t ) , u ( t ) ,yo
fW,40,yo
+ h(t), t l
fb,u l F[x, u ]
x(t3 = xo Then the control which provides the absolute minimum of
dt
l1
=
v ( t ) - = f[x(r, t ) , U(Y, t ) , Y,
bY
Y
6 1
t]
= u
2- (P22 + u') 2
dx dt
P(u) =
+
6
-=f
F [ x ( Y t, ) , U(Y, t ) , Y, tldr
bX
+ 1 - t,
The solution to the lumped-parameter problem with
for the system
bX bt
t
+ 1 - t,
Example 1. Consider again the system studied at the beginning of the paper, with
with
P(u) =
o < r < t t 6 Y 6 t
+ -
+ h ( t ) , tldt
for the system
dx - = dt
u*{r; t , - t Y, xo(r - t ) ) ; l l * { t ; tf, X,(Y - t ) ] ; u*(t; t 1 Y, Xo(Y - t ) l ;
Jot'Fdt
in this case is (Athans and Falb, 1966)
with
X(Y, 0) = Xo(Y); -h(t,) h(t,)
6
Y
6 1
6 1
is
U(Y,
We now specify, instead of Equations 2 and 3, an arbitrary 1. Thus, the entrance entrance function xo(r), for -tf 6 Y temperature behavior is considered to be arbitrary. Then, applying Equations 48 and 48' yields the solution to the original distributed-parameter problem:
>1 t,
cosh [ p ( l - r ) ] cosh G(t 1 Y)];
0 6 t 6 Y
+ -
Special Case For the important special case u ( t ) = 1, it follows that h ( t ) = h-l(t) = t. Then, the optimal transient and control become: Case I. t l 1
>
\ K ( Y ,t ) =
X(Y, t ) =
{
tanh G(1 P tanh G(t,
p
U(Y,
CASE11. t,
U(Y,
6
- Y)];
- -91;
-
0
6 t 6 t / - (1 - Y) (1 - Y ) 6 t 6 t ,
t,-
(1
6 t 6 t/
t ) = -K(r, t ) x ( r , t )
1
t) =
I
Case 11. t l 6 1
b(l - Y ) ] - . G(t + 1 - Y)]'
cash xoo(r
- t , cosh
\
t+l-t,1
Nomenclature
A B F(x, u)
x n constant matrix n x m constant matrix = functional, integrand of performance index = n =
VOL. 7
NO. 3
AUGUST 1968
421
= n-dimensional functional = = = =
= = = = =
= = = =
=
R1,R,, R a t , etc. = = r = t = = = = =
GREEKLETTERS
Hamiltonian gain matrix gain gain matrix for related lumped-parameter system stationary gain gain matrix defined by Equation 33 dimension of control vector dimension of state vector performance index performance index of related lumpedparameter system adjoint variable n X n positive definite constant matrix n X m positive definite constant matrix regions divided by characteristic lines normalized distance, re[O, 11 normalized time, t e [ o , t f ] or te[to, t i ] control vector minimum of performance index, Equation 19 velocity defined by Equation 19 state variable
= weighting factor in performance index
2t)
=
SUPERSCRIPTS
delta (impulse) function
*
=
T
= transpose
optimal
SUBSCRIPTS
f 0
= final time = initial time
literature Cited
Athans, M., Falb, P. L., “Optimal Control,” Chap. 9, McGrawHill. New York. 1966. Koppel, L. B., Shih, Y. P., Coughanowr, D. R., IND.ENG.CHEM. FUNDAMENTALS7, 286 (1968). Shih, Y. P., ‘‘Optimum Feedback Control of Tubular Process,” thesis for Ph.D. degree, - - Purdue University,. - West Lafayette, . Ind., August 1967. Sirazetdinov, T. K., Automation Remote Control 26, 1449 (1965). Voelker, D., Doetsch, G., “Die Zweidimensionale Laplace Transformation,” Verlag Birkhauer, Basel, 1950. RECEIVED for review August 30, 1967 ACCEPTEDApril 12, 1968 Work supported by a research grant from the National Science Foundation.
SENSITIVITY ANALYSIS IN OPTIMAL SYSTEMS BASED ON T H E MAXIMUM PRINCIPLE T. M. C H A N G ’ A N D
C . Y. WEN
Chemical Engineering Department, West Virginia University, Morgantown, W . Va
.
The sensitivity equations, necessary for determining the effects o f the first-order variations in parameters on the optimal objective function and decisions of a continuous process, are derived based on Pontryagin’s maximum principle. Methods for optimal design of parameter-sensitive systems are presented. The algorithm o f the maximum principle i s extended to problems with sensitivity constraints. The sensitivities of the optimal temperature and the maximum conversion owing to a single exothermic reaction taking place in a tubular reactor are calculated.
solution of an optimization problem is obtained based on set of numerical values of the parameters. Usually, the values of the parameters are subject to change owing either to the uncertainties in the experimental evaluation or to variations of the operating and surrounding conditions. Therefore, it is desirable to know how the optimal solution changes if some of the parameters are changed. T o answer this question, it is necessary to make a sensitivity study on an optimal system. So far, the sensitivity analysis has been made in a few classes of optimization problems such as the linear programming problem (Garvin, 1960; Hadley, 1962) and the convex quadratic programming problem (Boot, 1963, 1964). Since the maximum principle has been shown to be powerful in optimizing continuous processes, a sensitivity study on problems solved by the maximum principle is reported here. If the sensitivity analysis shows that the optimal performance of a system is strongly dependent on the selected parameter values which are subject to variation, the actual system perHE
Ta
1
422
Present address, Goodyear Tire & Rubber Co., Akron, Ohio l&EC FUNDAMENTALS
formance may substantially deviate from the specifications. To ensure better system performance over a range of parameter values, the parametric sensitivity must be taken into consideration in optimal system design. In other words, the decision variables of a system should he determined on the basis of not only optimality but also sensitivity of the system. I t is the purpose of this paper to derive sensitivity equations and show how they may be used to develop methods for the design of parameter-sensitive systems. Deflnition of Sensitivity
T h e sensitivity of x 1 to the variation of Let xi = f(wl). near w1 = al is defined as the ratio of the percentage change in .xl to the percentage change in w 1 a t w l = 2 j 1 (Bode, 1945). Let 4 x 1 = x 1 - 31 and 4 w = w1 where z = f ( @ l ) . T h e sensitivity of x1 to the variation of w 1 a t ~1 is expressed as
wl
where 4x1/4w1 is called the sensitivity coefficient.
