which contained about 50% solids, was fed to a crystal puripure p-dichlorobenzene was fication column and 99.5+% produced. Beer concentrates have been made by freeze concentration of beer in a pilot plant. A beer containing, for example, 4.9 volume % alcohol was cooled in the chiller to about 26’ F. to form a slurry containing 40 weight Yo ice crystals. This slurry was processed in a crystal purification column which yielded a high purity water product and a beer concentrate (mother liquor) containing about 6 weight 7 0 alcohol. T h e beer concentrate was reconstituted by adding water equal to the volume of water removed in the freeze concentration. An analysis of the original and reconstituted beer is presented in Table I. A taste panel judged the reconstituted sample equal to the original beer. A coffee extract which contained 28 weight % dissolved coffee solids was cooled to form a slurry containing 33 weight yo ice crystals. 1his slurry was processed in a crystal purification column which yielded a water product and a coffee concentrate (mother liquor) which contained 37 weight yG dissolved coffee solids. Laboratory evaluation of a sample of the concentrate showed that the coffee was essentially equal to the feed solution.
Table 1.
Beer Analysis
Color, ‘/*-inch cell, Lovibond, Series 52, O Apparent extract, Tc Real extract Alcohol, weight yc Alcohol, volume 9; Original extract, yc Real deeree of attenuation. % Reducing sugars, as maltose, 7c Protein ( N X 6.25), yo Dextrins, etc., Yr Iodine reaction Starch Amylodextrin Ervthrodextrin Iron, p.p.m. Acidity, as lactic acid, % pH value Diacetyl, p.p.m. Isohumulones, p.p.m. (Rigby Method 11)
Original
Reconstituted
2.2 2.63 4.32 3.65 4.87 11.42 62.2 1 03 D 28 3 29
2.4 2.98 4.76 3.85 4.93 12.22
Negative Negative Faint trace 0.03 0.12 4.2 0.06
Negative Negative Faint trace 0.05 0.14 4.2 0.08
16.3
15.0
61 . O
1 26 0 30
3.50
High purity p-xylene has been produced commercially in columns with cyclic solids movement. The units operate continuously with an on-line time of better than 95%. Operating labor is low with only one man per shift in multicolumn plants. Nomenclature
A . = cross sectional area. sa. ft. C i = average specific heat of crystals between T , and T,, B.t.u./’lb. O F. p. = gravitional conversion factor F = void fraction, dimensionless H , = heat of fusion of highest purity product, B.t.u./lb. R = fraction of a pound of reflux recrystallized per pound of crystals fed to column s = specific surface, sq. ft./cu. ft. T , = crystal feed temperature a t column inlet, O F. T , = melting point of highest purity product, O F. AP = pressure drop, lb./sq. ft. y = length of crystal bed, ft. p = mother liquor viscosity, lb./ft. sec. I
dB
A
= rate of flow, cu. ft./sec.
Literature Cited
(1) Arnold, P. M. (to Phillips Petroleum Co.), U. S. Patents 2,540,083, and 2,540,997 (Feb. 6, 1951); Reissue 34,038 (July 19, 1955). ( 2 ) Hachmuth, K. H. (to Phillips Petroleum Co.), Ibid., 2,764,878 (Oct. 2, 1956). (3) Harrison, R. C., Cottle, J. E., McKay, D. L., 111th Dechema Colloquium, Frankfurt, Germany, March 1964. (4) hlcKay, D. L., Goard, H. W., Chem. Eng. Progr. 61, No. 11, 99 (1965) ,-_-( 5 ) Schmidt, J. (to Phillips Petroleum Co.),U. S. Patent 2,617,274 (Nov. 11, 1952); Reissue 23,810 (March 30, 1954). (6) Thomas, R. W. (to Phillips Petroleum Co.), Ibid., 2,854,494 (Sent. 30. 19581. ( 7 j \l\jeedm’an, J.‘A. (to Phillips Petroleum Co.), Ibid., 2,747,001 (May 22, 1956). ( 8 ) Winnick, Jack, Powers, J . E., A.I.Ch.E. J . 7, No. 2, 303 (1961). ( 9 ) Woolcock, J. W., 111th Dechema Colloquium, Frankfurt, Germanv. March 1964. (IO) Yagi, Sakae, Inoue, Hakuai, Sakamoto, Hirohiko, Kagaku Kogaku 27 ( 6 ) , 415 (1963). , I
RECEIVED for review May 2, 1966 ACCEPTED October 17, 1966 Division of Industrial and Engineering Chemistry, 151st Meeting, ACS, Pittsburgh, Pa., March 1966.
PERIODIC PROCESSES: A VARIATIONAL APPROACH F, J
. M. H0 RN
processes in the chemical industry are carried out either in the steady state or batchwise. I n this paper a inore general class of processes, the periodic processes, is considered. This class exhibits some properties which may be of great practical value. T h e theory of periodic processes includes the theory of steady-state and batch processes, since these two types can be regarded as special cases of periodic processes. Formally, periodic processes are similar to some types of steady-state recycle processes, because a cycle in time OST
A N
D R
.C.
L IN
,
Rice University, Houston, Tex.
and a cycle in space lead to very similar mathematical representations. Some general results concerning the calculation and optimization of periodic processes are derived. Parallel to this, the special example of a chemical reaction process carried out in a continuously operated tank reactor is considered. The only objective in the construction of this example was to make it very simple and yet retain the possibility of demonstrating some of the interesting properties that periodic processes may VOL. 6
NO. 1
JANUARY
1967
21
The class of periodic processes contains as subclasses the batch and the steady-state processes. Simulation and optimization of periodic processes are discussed. In particular, the conditions are explored under which an optimum steady-state process can be improved by periodic control. Some of the general results are demonstrated by means of the example of a chemical reaction process.
possess. The authors d o not claim that this example has a practical application. Steady-State and Periodic Processes
T h e dynamic behavior of many chemical systems can be represented by a set of differential equations : dx 1 -= dt
fi(X1,
xz, . . . . x,;
i
= 1,2,
y1,y z ,
, , , ,
y,)
.... n
where the x’s denote state and the y’s denote control variables. As an example, consider a continuous~y operated tank reactor in which the reactions
there may be none, one? or more than one steady state belonging to a given control, a steady state may be unstable or stable, and in the latter case this steady state is approached by the system from appropriately chosen initial states. The example represented by Equations 3 has one steady state for each possible control y in the physically admissible regions (XI 2 0, x2 2 0, y > 0) and this steady state is stable. Thus if this system is operated a t constant temperature (corresponding to constant y), the point representing the state in an XI, x2 diagram will move along a trajectory leading from the initialstate point towards the steady-state point. Another mode of operation is to use periodic control functions j(t)-that is, functions satisfying the relations : yj(t
A+B (2b)
are carried out (see Figure 1). Suppose Reaction 2a is of the a t h order. Reaction 2b is of first order, and the temperature dependence of both reaction rates is given by the Arrhenius equation. Furthermore, suppose that flow rate and composition of the feed material are given and constant and that the temperature in the tank can be controlled. Then the dynamical behavior of the system is described (see Appendix 1) by:
dt
-yxp
=
7)
- ayPx1 - x1
+I
+
7) =
i In this example there are two state variables (n = 2) and one control variable ( r = 1). If the control variables are given functions of time, and if the initial values ~((0)of the state variables are known, the state of the system can be calculated as a function of time by integrating differential Equation 1 or 3. Of special interest to industry are steady-state solutions-that is, solutions of the equations which correspond to time-independent state and control variables:
. . . . x,; i
=
y1,yz,
. . . .y,)
= 0
(4)
1, 2, . . . n
For such constant control, the steady-state values of the n state variables have to satisfy n ordinary equations. Of course:
9 A
22
q
F
m A , B and C
l & E C PROCESS D E S I G N A N D DEVELOPMENT
, . .
(5)
r
=
xi(t) for any t
(6)
1, 2, . . . n
are satisfied. Figure 2 shows qualitatively the result of a computer calculation carried out for the example represented by Equations 3. T h e periodic control function is shown in Figure 2a, and the trajectory in the phase diagram belonging to it is shown in Figure 26. This trajectory approaches asymptotically a cycle, the period of which is equal to the period of the control function. This is what one intuitively lrould expect to happen for the system in question. The cycle of Figure 2 corresponds to a periodic process. One method of computing such a cycle is to simulate the approach to it as shoxvn in Figure 2. Any other method of integrating the system of differential equations (Equations 3) subject to the boundary conditions
x i ( 0 ) = xi(.)
I
for any t
In the following discussion 7 does not necessarily denote the smallest period-that is, that T / 2 or r / 3 , etc., is, also, a period is not excluded. Obviously, a constant control is a trivial case of a periodic control because it satisfies Equation 5 for any T. If the system is subjected to a periodic control, possibly the state variables, also: become periodic functions of time. If the state variables are periodic functions, and the period is the same as for the control variables, the process is called a periodic process. In this case, as \vel1 as Equation 5, the following relations xi(t
fi(X1, x z ,
= yj(t)
j = 1,2,
A-C
dx 1 -
+
i
=
1, 2, . . . . n
(7)
can be used, also, for the computation of the cycle. From Equation 6, it follows that Equation 7 is necessary for a periodic process. If the control is periodic, Equation 7 is also sufficient for periodic process, since the state a t the beginning of a period together with the control function in this p‘riod determines the state variables as functions of time throughout the period. As in the case of steady-state operation, there may be more than one periodic state for a given periodic control, and each of the corresponding periodic processes may be either stable or unstable. These problems, and also the question of the existence of a periodic process belonging to a given periodic
a
0.8
1.6
2.0
2b I .9
x2 1.8
I .7
AI
Figure 2.
Approach to limit cycle for periodic control
6 .O
6.0
----------4 .O
7th STEP M i 18.66%
Y 2.0
Y
0.5 1.0
2.0 3.0
4.0 5.0
t Figure 3.
Y 2 .o
2.0
---- -- - ----0
42)
0.5 0
LO
2.0
3.0 4.0 5.0
0.5 0
ID
t Results of hill climbing in function space for
control, which can be attacked by means of fix points theorems, are not treated in this paper. For practical applications, if a solution of differential Equations 3 subject to Condition 7 exists, the process can be operated as a periodic process. Such a mode of operation is distinguished from others because, in a periodic process, both control sequences and output sequences are repeated again and again, and this is of practical advantage. T h e term “control of a periodic process” may refer to the periodic control functions, y ( t ) , which generate the periodic behavior, or to additional control, which is used to counteract disturbances, or to keep the process in its periodic state if this periodic state is unstable. I n this paper problems concerned with the latter meaning are not discussed.
2.0
3.0 4.0 5.0
t T
= 5.0
Relations among Batch Processes, Recycle System, and Periodic Process
A batch process is defined as a process in which the initial conditions for the differential Equations 1 describing the process and the process time, T , are given or can be set. T h e solution of the equations describing a batch process corresponds, therefore, to an initial value problem. I n industry, batch processes are often repeated consecutively with identical initial conditions each time. This is a trivial case of a periodic process. An example is a tank reactor which is periodically emptied completely and refilled with fresh raw material of a given composition (no continuous input and output in this example). Then the state at the beginning of the period is determined by the composition of the raw material and is inVOL 6
NO. 1
JANUARY 1967
23
dependent of the state a t the end of the previous period. Consider now a process in which the reactor is not completely emptied a t the end of each reaction period, so that a part of the content of the reactor remains there and is mixed with fresh raw material a t the beginning of the next reaction period. In this case, the initial state of the control period depends on the final state of the previous control period. If the process is carried out as a periodic process, the initial states of consecutive control periods are equal and the relations
Ax denotes a column vector and dg/dx(.r) is the matrix with element dg,/dx,(7) at the intersection of the ith row' and j t h column. The relation between Ax**(.) and Ax*(O) can be found by linearizing the differential Equations 1 (see Appendix 2) :
A(0) Ax*(O)
Ax**(r) =
(12)
A ( t ) is an n-dimensional square matrix satisfying
are satisfied. Here the gz denote the state transformations corresponding to removing a part of the tank content and filled up with fresh raw material. The situation described by Conditions 8 will be approached naturally in many practical cases if the sequence of operations is continued long enough. In any case, if differential Equations 1 have a solution subject to Equation 8, the system can be operated as a periodic process. Equations 1 and 8 describe a general periodic process in which a continuous operation (corresponding to Equation 1) and a discrete operation (corresponding to Equation 8) are repeated cyclically. This description contains, as special cases, the previously discussed examples of a continuously operated tank reactor and the pure batch process. In the is, Equation 8 former case, g is the identity function-that becomes Equation 7 ; in the latter case, g is identically constant-that is, the initial conditions become independent of the end conditions. Equations 1 and 8 also describe some steady-state recycle processes. For example: consider a tubular reactor for which the input depends on the output due to a recycle of material or heat. Then Equation 1 describes the change of statee.g., concentrations and temperature-along the reactor ( t is now a position variable), ~((0) refers to the state at the inlet, and x i ( . ) refers to the state a t the outlet. Equation 8 describes the input-output relationship caused by the recycle. Newton-Raphson Method for Calculation of Periodic Processes
T h e difficulty in calculating periodic processes is the matching of Boundary Conditions 8. This problem can be solved by a Newton-Raphson iteration, provided sufficiently good first approximations of the state variables at t = r corresponding to periodicity are known. Suppose XI*(.), x ~ * ( T ) , . , , . . , x,*(r) are the approximate state variables at t = 7. Equation 8 then yields the initial values
xi*(0) = gi(xl*(r), x2*(7), . , . .
xn*(7));
i = 1, 2, . . . n
(9)
The system of differential Equations 1 can be integrated with these initial values; thus the end state, x * * ( r ) , can be determined. This end state will not be equal to x * ( T ) because of the approximate character of x * ( T ) : xi**(.)
- xi*(.)
= di;
i
= 1, 2,
...
TZ
(10)
If xi*(.) is changed by a small amount, A x i * ( . ) , it follows from Equation 8 that
dA = - A ( t ) df dt dX
or in index notation (summation convention)
From Equations 11 and 12 it follows that:
The objective is to choose Ax*(.) so that:
xi**(.)
+ Axi**(r) -
[xi*(.)
-I- A x $ * ( . ) ]
=
0
because then the boundary conditions are satisfied. From Equations 10, 14, and 15
+
Since second-order terms have been neglected, x * ( T ) Ax*(.) is not the exact end state corresponding to periodicity but only an improvement of the original state, x * ( T ) , if this state was a sufficiently good approximation to justify the neglect of the second order terms. The improved value can be improved again by the same technique. Thus, the scheme of the Newton-Raphson iteration is as follows: STEP 1. Choose a starting point, x*(.), which could be obtained by simulating the approach to the periodic state. Since the Newton-Raphson method converges much faster than the simulated approach, it is better to switch over to the Newton-Raphson technique at some point. STEP 2. Calculate x * ( O ) from Equation 8 and integrate the differential Equations 1 from t = 0 to t = T with this initial condition. Calculate d from Equation 10. STEP3. After step 2 the quantities djf,/dx, are known functions of time and differential Equations 1 3 can be integrated. This corresponds to finding the complete solution of a system of homogeneous first-order differential equations. STEP4. Now the correction AX*(.) can be calculated from Equation 16 and x*(.) can be improved by applying the correction to the initial choice. The corrected values can be used for a new iteration (start again at step 2), and this procedure can be repeated again and again until a significant improvement can no longer be obtained. This method has been used in the examples discussed in this paper. Optimization of Periodic Processes
The summation convention is used here, and in the following text-that is, occurrence of indices in pairs ( jin Equation 11a) indicates summation over these indices. Equation 11a can also be written in vector notation: 24
I & E C PROCESS D E S I G N A N D DEVELOPMENT
(15)
i = 1,2, , . . . n
sor
Assume that the objective function, M , is given by:
M =
1 -
m[x(t)]dt
T h a t is, M is the average over one period of an instantaneous profit, rn? which depends on the instantaneous state. For instance: in the example of the reaction process described by Equations 2 and 3 M is given by:
M
=
lr
x2(t)dt
if product B is the desired one and the total production of B is the only economically interesting quantity. T h e problem then is to find a process for which M is a maximum. For a given periodic process M has a definite value. The first step in optimization is to investigate how this quantity changes if the periodic control varies by a small amount. This problem can be solved by the methods of variational calculus (see Appendix 3) and the result is given by the following equation (summation convention) :
Here 6M is the change of the objective function due to small changes, 6 j 2 ( t ) , of the control functions. The variables h,(t), the so called adjoint variables, satisfy the equations:
ment procedure (steps 2 to 4). This is repeated until a significant improvement can no longer be obtained.
For practical calculations steps 2 and 3 are intimately related. A solution of the adjoint Equation 20a can be found in the following way. A particular solution, A p ( t ) , of Equations 20a is determined by integrating 20a for any arbitrary initial or end condition, say X*P(O) = 0
i
(224
= 1,2, . . . . n
or h,P(7) =
i
0
(22b)
= 1,2, , . . . n
Experience shows that for stable physical systems the integration of the adjoint equations in the direction of decreasing t is usually numerically more stable and the following discussions therefore refer to boundary Conditions 22b. The matrix function, A ( t ) , of Equation 13 represents a complete solution set of the homogeneous equations corresponding to 20a. Therefore, the general solution of Equations 20a is given by :
and the boundary conditions: where form : Only periodic processes are considered here. If the control of an existing periodic process is changed slightly, some time Mill elapse, in general, before the process reaches its new periodic state. This transition is not considered here. 6M refers to the difference in the objective functions between tMo periodic states. E,quation 19 makes it possible to improve a given periodic process. Suppose that in a given periodic process controly(t) is changed t o j + ( t ) by:
j
= 1,2,
,...r
and that E is so small that the second-order terms in Equation 19 can be neglected. Then the integrand in Equation 19 becomes :
and is therefore positive. Thus, the change from y to y + increases profit M and improves the process This improvement can be carried out consecutively until a point is reached a t which no further significant improvement can be obtained. This method is known as “hill climbing in function space” and has been successfully applied to many variational problems ( 7 , 2). The adaptation of this method to the optimization of periodic processes may be outlined as follows : STEP1. Choose a control, ? ( t ) , in the interval 0 < t < 7. STEP 2. Calculate a periodic process belonging to the assumed j ( t ) by the methods outlined above. are now known, a solution STEP3. Since funciions dff,/bx, of Equations 20 can be determined. STEP4. The improvement of the policy can be carried out now by means of Equation 21. (For the determination of E see the references given.) STEP 5. The improved policy is used for a new improve-
c1,
c2,
. , . .c,
are integration constants.
X(t) = XP(t)
In vector
+ cA(t)
(24)
( A , c row vectors)
Putting Equation 24 into Equation 20b yields:
The adjoint variables are then given by:
Thus, the matrix function, A ( t ) , is used in both the matching step 2 and the improvement step 3. This means a simplification in programming the iteration scheme. Moreover, the amount of calculation can be reduced considerably if in matching step 2 h ( t ) from the previous improvement step is used. This is possible because the corrections to the control functions are usually small, and, therefore h ( t ) is not changed much. The correctness of the match, of course, is not affected by a n inaccurate A ( t ) . Only the convergence of the NewtonRaphson method is affected by this. Since the improved limiting cycle is close to the one in the previous step, an excellent first guess for x ( 7 ) is known, so that convergence is not a critical problem. Necessary Conditions for Optimality
For the optimum control, from Equation 19 : (27) if the control is unrestricted. Since the so-called strong maximum principle is valid in the case in question (see Appendix 3), the much stronger statement can be made that the expression X f a must be a maximum with respect to admissible VOL. 6
NO. 1
JANUARY 1 9 6 7
25
(if there are restrictions on 1) variations of y for any t. In this paper only a weaker version of this will be needed. If the control is unrestricted, it follows from the strong principle that
If the steady-state process is considered as a special case of a periodic process, the adjoint equations are given by Equation 20a. The transformation g becomes the identity transformation, so that the boundary conditions for the adjoint variables are :
h i ( 0 ) = hi(^)
i s necessary for an optimal control corresponding to a maximum
of the objective function. Equation 28 is written for a onedimensional control. In case of a control vector of higher dimensionality Equation 28 is to be replaced by the condition of negative indefiniteness of the Hessian matrix of X J l at any t. If Equation 27 is satisfied at any point where y is unrestricted (Appendix 3),
rn
+ h J , = constant
(30) If T is unrestricted, this derivative vanishes for the optimum period. If the control is also optimal, then from Equation 29 rn
-+-
hJi
=
M
for any t
,
n
(37)
The differential Equations 20a admit a constant solution because :
%
and
a x1
brn i = l , 2 , . . . n bxi = 1,2, . . . . n -
are time-independent for a steady-state process, and the solution is given by:
(29)
within the period. A necessary condition for optimum period T can be obtained by differentiating the objective function with respect to T . This derivative is given by (Appendix 3) :
i = 1, 2, . , .
brn
hi = - - A k i dXl;
i = 1,2,
, . . .
n
The constant solution obviously satisfies the condition of Equation 37. From Equations 27 and 38 Condition 36 is obtained again. This means that the first necessary conditions for an optimum periodic process are automatically satisfied by an optimum steady-state process. No conclusions as to the possibility of improving the steady-state process by cycling can be drawn from this condition. Conditions 28 become :
(31) (39)
improvement of an Optimum Steady-State Process
Consider a process which is operated a t steady state with given time-independent control, y 3 . The objective here is to make the instantaneous profit as large as possible. If the control is unrestricted, optimum operation implies :
brn -
=
0 j = 1,2, . . .
If the steady-state process is optimal lvithin the class of steadystate processes, the second derivative of the objective function with respect to the control (one-dimensional control is assumed now) must be negative. The second derivative can be obtained by implicit differentiation. The condition which is satisfied by the optimum steady-state process is given by:
7
3Ji Application of the chain rule yields: where
(33) (41 a)
Differentiation of Equation 4 with respect toy, gives:
bfi
ax,
-+
._
- aft =o
i
= 1,2,
.... n .... r
(34)
If the matrix bf/bx is nonsingular, the derivatives of the state variables with respect to the control variables can be determined from Equation 34.
-_-
bXk 4
1
-
bf* k byy3
=
Ak%-
=
1, 2, 1. 2,
n
r
(35)
where A k z is an element of the inverse matrix of b f , / d x k .
For a n optimum steady-state process with unrestricted control
The question then arises whether the optimum steady-state process can be improved by cycling. I n order to investigate this, the steady-state process is regarded as a periodic process (it is, in fact, periodic for any arbitrarily chosen T), and the conditions of optimality in Equations 27 and 28 are applied. 26
l&EC PROCESS DESIGN A N D DEVELOPMENT
Condition 40 is different from 39. Therefore, an optimum steady-state process which satisfies 40 may violate Condition 39. If this is the case, the optimum steady state can be improved by cycling. The physical reason for the difference between the two conditions is the following: If Condition 40 is violated, a timeindependent variation of the control improves the process ; while if Condition 39 is violated, a sufficiently small impulse disturbance of y ( t ) (compare Figure 4b) improves the process. This follows from the derivation of the strong maximum principle ( 3 ) . From Equation 41 :
B 1 = 0 iff, is a sum of two functions one of which depends only on the state and the other only on the control (separable control). B B = 0 iff, is linear in the state variables, though the coefficients of the state variables may be arbitrary functions of the control.
M -18.57 % 2 .o
2.0
LOWER BOUND
0.5
0.5
.02 .04 .06 .OB .IO
0
0
t
Figure 4. Results of hill climbing in function space for 7 == 0.1
.02
.04 .06 .08 .IO
0
.02
t
.04 .06 .08 .IO
t
----_-------_ 4.0
M.18.71 Ve
M 19.1I X
Y 2.0
2 .o
------ -- -- - -
0.5
.02 .04 .06 .08 .IO
0
t
If all three terms are zero, Condition 40 becomes identical to Condition 39. Both conditions are also identical if there is only one state variable, because then Condition 36 implies t h a t B I = B2 = B3 = 0 . Discussion of a Reaction Process
Consider:the reaction process described by Equations 2 and 3. In the steady state the following equations are satisfied :
-
a2,pxI
p:1Q
-
- xi x2
=
+1=0
(424
0
(42b)
For anyy there is only one solution ( x l 1 xq) with The objective is to maximize %?-that is,
cup
2 0, xe 2 0. (43)
m = xg
by choosing y appropriately. restricted y only if
x1
This has a solution for un-
>1
(44)
otherwise, x 2 increases monotonously withy. If Equation 44 is satisfied, the objective function has a stationary maximum with respect to control -Y or, equivalently, temperature T . I n this case the conditions in Equations 36 and 40 are satisfied for optimum operation. However, by evaluating the magnitudes ' 4 , j it can be shown that Condition 39 is violated if
(45)
P < l
Therefore, if Equations 44 and 45 are simultaneously true, the optimum steady state of the system can be improved by cycling. In order to demonstrate this an example has been calculated numerically with the p,arameters:
(46) Figure 3 shows the result of hill climbing carried out for a period T = 5.0. The starting point is the optimum steadystate control (Figure 30). This control policy corresponds to a saddle point in func1:ion space. The first variation is zero and for some changes in the control (corresponding to remaining in the steady state) the objective function decreases, while for other ones (corresponding to an impulse disturbance) (Y
= 2.0
p =
0.75
a = 1
0.5
.02 .04 .06 .08
.IO
t
B3 = 0 if the instantaneous profit is a linear function of the state variables.
-2,x1Q
0.5 0
0
.02
.04 .06 .08 .IO
t
the objective function increases. In order to depart from the saddle point, the first step in the iteration scheme discussed previously has to be modified. A relatively large change of the control (the sign does not matter) is carried out in a small interval (Figure 36). After this, hill climbing proceeds as described previously. e has been chosen constant with respect to t and the step length policy adopted was as follows: If a step does not improve the objective function, e is automatically halved and the step is repeated. T o avoid having an unreasonably small e, this value is tentatively doubled after the execution of three successful improvements ( 7 ) . Figure 3, b to f, shows results for the indicated number of improvement steps. The corresponding values of the objective functions are also given in Figure 3. I n this problem upper and lower bounds for y have to be set. Otherwise the oscillations of y grow bigger and bigger from step to step. Figure 4 shows hill climbing results for a smaller period ( T = 0.1). I n Figure 5 corresponding cycles are plotted in an x 1 - x 2 diagram. Point P corresponds to the optimum steady state. The upper cycle corresponds to the periodic process belonging to the control shown in Figure 4, step 16, and the intermediate cycle belongs to an intermediate control. By the hill climbing procedure point P is transformed step by step into the upper cycle. The curve through P is the locus of steadystate points belonging to different time-independent controls y . In the case in question, there are cycles completely above this curve, because of the special nature of the problem. I n other problems the optimum cycle may lie partly below the steady-state curve. The improvement of the steady-state process leads to relatively short cycles-that is, to periods which are small compared with a characteristic time constant describing the dynamic behavior of the system. Though the system is nonlinear, such a characteristic time can be defined, for instance, by linearization about a steady state. In the problem in question this characteristic time has the order of magnitude 1 (t is dimensionless, see Appendix 1). The results of Figures 3 and 4 indicate that the periods chosen should be much shorter. If in the hill climbing procedure T is chosen large, the improvement routine automatically leads to short cycles (see Figure 3). In the calculations of Figures 4 and 5 the period has been chosen a priori relatively short, as can be seen from the fact that the cycles of Figure 5 are relatively small. The physical reasons for this are given below. '401. 6 NO. 1
JANUARY
1967
27
Figure 5. Steady-states for various temperatures (curve with point P ) and cycles corresponding to periodic control for tank reactor example Q. corresponds to infinite fast switching between temperature limits
0.2 Calculations have been carried out starting with an arbitrary periodic control. Use of the impulse policy for the first step is not necessary then. For a given period the results obtained were essentially independent of the starting policy, except when the starting policy was a nonoptimal steady-state control. Then the optimum steady state was approached. Now the above results will be made physically plausible. T h e two rate constants k l and k 2 are plotted against each other in Figure 6. (Instead of k l and k2, y and y p may be plotted.) The heavy line is the locus of points corresponding to different temperatures :
a contour line of rn. The significance of the conditions cup > 1 for a stationary optimum now becomes apparent. If cup > 1 and p < 1 (Figure 6 6 ) , the convex hull of the attainable locus contains points corresponding to a higher profit than the maximum profit on the attainable locus. Points in the convex hull can be obtained by a sufficiently fast periodic control. For instance, consider points P1 and P2 corresponding to k l l , k 2 , and kI2, k?. The coordinates of Q are then given by:
(47)
Consider now a system which is operated periodically for two values of the control variables, j 1 and j 2 , corresponding to P’and P2. During the time interval A t l controljl is chosen and during the time interval At2 control j 2 . At1 and A& satisfy:
This curve is convex downward for p > 1 and concave downward for p < 1. For the steady state the objective function rn = xp can be expressed as a function of k l and kp. T h e contour lines of this function (broken lines) are convex for a < 1 and concave for CY > 1. O n each such line rn is constant. The general direction of increasing m is indicated by the arrow. A stationary steady-state optimum corresponds to a point of contact between the locus of attainable points (heavy line) and
Figure 6. Attainable points (heavy line) and contour lines of the steady-state objective function (broken lines) in rate constant space Arrow indicates direction of increasing proflt. Q corresponds to inflnite fast switching between P1 and P*
28
I&EC PROCESS DESIGN A N D DEVELOPMENT
(49)
If period r is sufficiently small, the state at the end of the period is related to the state a t the beginning of the period by
if y ( t ) is not varied. relation :
From this by integration by parts the
(57)
h(0)6x(0) = A ( T ) ~ x ( T )
If there exists a state satisfying
can be obtained if h ( t ) is chosen such that:
T h e system can be operated arbitrarily close to this state (a quasi-steady state) by making d sufficiently small and operating according to the switching policy given by Equation 49. Since, in the example in question, functions fi are linear in the rate constants, the quasi-steady state would result from steady-state operation with rate constants kl, kL corresponding to a point on a line connecting P1and P2. Now the reason why fast switching between high and low temperatures (see Figures 3, 4, and 5) gives good results can be understood. T h e optimum quasi-steady state corresponding to infinitely fast switching between the limits of y is represented by point Q in Figure 5. Very fast switching of the temperature is, of course, not a practical proposition. For this reason, calculations have been carried out to find a n optimum heat flux control. I n these calculations, the system is described by three state variables instead of two, because the temperature now becomes a state variable. Infinitely fast switching does not bring any improvement in this case. With the methods described in the previous sections an optimum control cycle and period T can be found. The results will be published in another paper.
This linear differential equation admits n independent solutions which form the rows of the matrix A ( t ) . Equation 57 then becomes : A(O)6x(O) = A ( T ) ~ x ( T ) >
from which Equation 12 follows. Appendix 3
The identity:
is used.
For constant
T
it follows that:
Then integration by parts yields:
Appendix 1
T h e reaction rates, by :
rl
%
+A
and rz, for Reactions 1 and 2 are given
lT
h $ 6 ~dt
(60)
If h ( t ) is chosen so that Equations 20a are satisfied, the first term vanishes. If Conditions 20b are satisfied, the second term vanishes, since for periodic processes :
T h e material balances for components A and B are given by: Thus, Equation 60 becomes Equation 13. the expression :
(53)
H=Xf+m with respect to t yields
Equation 3 can be obtained by introducing in Equations 53 the dimensionless variables
Differentiation of
- -dt
if + x
g
x
(62)
m x. + x dfay j, + aax -
From this and from Equations 1 and 20a
af
dH - _- A - Q dt aY
(54)
For the optimum process
af
X - 6y
3Y
Appendix 2
sur sor
From the identity: x(t) [L
- f ( ~y ,) ] d t = 0
which is valid for any A(t), one derives:
A(t) [:6k
-
6x
af
1
dt = 0
(55)
5
0 for admissible 6y
Since y and -Q must be proportional to an admissible 6y, it follows that d H / d t = 0 and H = constant for optimum control. Consider now a variation of 7. I n this case ~ x ( T )in Equation 61 has to be replaced by: 6X(T)
+
f(T)6T
The first part refers to the variation of the end state due to a change of y or of the initial state; the second part refers to the VOL 6
NO. 1
JANUARY 1 9 6 7
29
change in period
T.
Thus, Equation 61 has to be replaced by
d
= column vector representing errors in matching
f
= defined by Equation 1
boundary conditions (see Equation IO)
sor
From the Identity 59 and Equations 20 and 64
6
g
= state transformation corresponding to discrete
H
= Hamiltonian, defined by Equation 62 = instantaneous profit
process
m
M n r
mdt = [m(T) 4- A ( T ) ~ ( T ) ] ~ T
= profit averaged over period = number of state variables
= number of control variables = time
t 6ik
= column vector representing state variables = column vector representing control variables = 1 fori = k, 0 fori # k
A
= row vector representing adjoint variables = matrix representing a complete system of solu-
X
Y If the control is not varied, or for optimum control,
x
tions of homogeneous adjoint equations = period
I
and from this, Equation 30 follows. The strong maximum principle and the weaker version of it (Equation 28) are valid because the problem in question can be reformulated such that the states a t t = 0 and t = T occur separately in the restrictions. T h e strong maximum principle has been proved by Pontryagin et al. (3) for this case. I n order to reformulate the problem a new set of n variables, xn+ 1 , x,+ 2, . . . . . . . xzn, is introduced for which the functions f are chosen zero: fn+l
= fn+z =
..., ... = f2n = 0
That is, the new variables are time-independent. 8 can then be written:
Restrictions
REACTOR EXAMPLE = defined by Equations 54
El,
E2
hi, hz
q
n, 72
R t t’ T V
x1,
Y CY
P
= = = = = = = =
= =
= x2
=
= = =
concentrations of substances A and B in tank inlet concentration of substance A activation energies frequency factors volumetric flow rate reaction rates universal gas constant dimensionless time defined by Equations 54 real time temperature reactor volume state variables defined by Equations 54 control variable defined by Equations 54 order of reaction ratio of activation energies (see Equation 54)
literature Cited Acknowledgment
The authors are indebted to Roy Jackson, Edinburgh, for many valuable suggestions. Nomenclature
A
= inverse of matrix
df/dx
B1, B z , Ba = defined by Equations 4 c
= row vector representing constants of integration
(see Equation 24)
(1) Horns, F. Trottenier, U., Chem. 2.9. Tech. 32, 382 (1960). (2) Kelley, H. J., “Methods of Gradients in Optimization Techniques,” G. Leitman, ed., Academic Press, New York, 1962. ( 3 ) Pontryagin, L. S., BolytonskiY, V. G., Gamkrelidze, R. V., Mishchenko, E. F., “Mathematical Theory of Optimal Processes,” Interscience, New York, 1962.
RECEIVED for review May 2 1966 ACCEPTED November 15, 1966 Division of Industrial and Engineering Chemistry, 151st Meeting, ACS, Pittsburgh, Pa., March 1966. \Vork performed under systems grant from the National Science Foundation (GU-1153).
PERIODIC COUNTERCURRENT PROCESSES F . J . M . H 0 R N , Department of Chemical Engineering, Rice University, Houston, Tex. The over-all stage efficiency of a periodically operated distillation column depends in a complicated way on the number of stages in the column as well as on equilibrium and transport parameters. Simple and accurate asymptotic formulas for the stage efficiency have been derived which show that the performance of a distillation column can b e improved considerably by periodic operation.
THIS reports a continuation of previously published work ( I ) on periodic extraction processes. The study of the asymptotic behavior of stage efficiencies is extended to a more general class of stages and in particular to vapor-liquid countercurrent processes. Consider a column which is operated in the steady state or in a periodic state [compare (2) and the following sections] and 30
I & E C PROCESS D E S I G N A N D DEVELOPMENT
in which two phases, called liquid phase and vapor phase, flow in opposite directions with average flow rates p and q . Both phases consist of two components, one of which is called the reference component, and its mole or mass fractions (depending on whetherp and q denote mole or mass flow rates) in the liquid and vapor phase are denoted by x and y, respectively. It is assumed that p and q can be regarded as constant along the
