ADAPTIVE OPTIMIZATION
OF CONTINUOUS PROCESSES G. E. P. BOX’AND JARAYAJAN CHANMUGAMZ Unioersity of Wisconsin, Madison, Wis., and Princeton University, Princeton, N . J .
Frequently the best operating conditions, x, for a chemical
Conventional Control for Optimization
process depend upon the level of some uncontrollable and
A serious difficulty in attempts to maintain optimum conditions arises when important variables are not measurable while the process is running. When this happens, control in its usual sense cannot be applied on these variables. Examples of such variables are: catalyst activity, mixing patterns, ambient external conditons, flow and temperature profiles, permeability of packed columns, abrasion of particles in fluidized beds, and fouling of heat exchangers. Changes in the le,vels of these variables cause changes in the response surface relative to the coordinate space defined by the (controllable) process variables. This in turn usually results in a shift in the optimum conditions.
unmeasureable variable, c, such as catalyst activity which varies over time in some unknown way.
In these circum-
stances, a method of adaptive optimization is required in which the operating conditions, x, are continuously modified
so as to be best for the level, c, existing at that time. Starting from the simple concept of evolutionary operation, an analogous method is developed.
The influence
of the
character of the noise and o f the process dynamics is discussed.
A detailed analysis of the effects of process dy-
namics in the case of consecutive reaction systems i s given.
A practical adaptive optimizer, which is the basis for the apparatus at the University o f Wisconsin, is described.
are referred to as responses in this independent variables as process variables. For instance, in a chemical process typical responses would be yield, conversion, and purity, and could include more complex quantities such as cost per pound of product and profit. Typical process variables are temperature, concentration, pressure, and flow rate. A plant optimum refers to a particular state of operation that leads to levels of responses which are in some sense jointly best ( 4 ) . I n simple cases it may be sufficient to consider maximization or minimization of a single response. Thus, optimum conditions may be defined as those giving highest yield, lowest cost of operation per pound of product, or maximum profitability. The functional relationship between the response of interest on the one hand and the process variables on the other hand is called the process function. Often the true nature of the functional relationship between the response and the process variables is unknown. I n such cases, empirical graduation of the process function such as by “response surface” methods ( 3 ) can provide information concerning the local nature of the relationship. EPENDENT VARIABLES
Darticle and
1 Present address, Statistics Department, University of Wisconsin, Madison, Wis. 2 Present address, Engineering Department, International Finance Gorp., Washington, D. C.
2
I&EC FUNDAMENTALS
As an example, let the response, 7, of a plant be a function only of the flow rate, x (process variable), and the catalyst activity, c (uncontrollable variable), as shown in Figure 1. The problem is to adjust x so that maximum response is attained at any level of catalyst activity, c. If the process is tvorking at flow rate x = X Z , the state of the catalyst activity being c = CZ, response 7 would be at its maximum value. If the catalyst activity now changes to c = c1 or c = c3, response 7 is no longer at its maximum. If there is some way of measuring catalyst activity c while the process is running and if the functional relationship 9 = f ( x , c) is known, then it is possible to adjust the level of the controllable process variable, x , so as to maximize 7 again. If c = c1, then x must be adjusted to the new level, x = x i ; if c = 63, x must be adjusted to x = x 3 to obtain optimum response. This “table look-up’’ method is essentially the basis of many current schemes of computer optimization, but is barred if the level of c is unknown. Evolutionary Operation
-4short digression is now made to discuss a method of process improvement known as evolutionary operation ( 2 , 5 ) . \Ye then show how the principles used in this technique may be adapted to provide a method of automatic optimization which does not require detailed knowledge of the process function or measurement of variables such as c. In the method of evolutionary operation a routine of plant Operation is set up where the levels of the process variables are changed by small but discrete amounts in a prearranged pattern based on the statistical theory of the design of experiments (7: 8). When evolutionary operation is used on continuous processes, the time interval between changes is so chosen that, for a given set of conditions, a sufficient period is allowed for the transients to die down and for the steady-state levd of the response to be achieved. T o avoid interference with regular production, only small changes are introduced, such as would produce effects of the same order of magnitude as the random errors inherent in the observations of the response. By averaging the observations from a number of experiences of the repeated
pattern, the local nature of the steady-state process function can be determined. From this information the form that the changes in the levels of the process variables should take to achieve a better response is deduced. The method is normally employed to check process operating conditions systematically (in groups of two or three variables at a time) over a period embracing the whole life of the process. I n this use the method is concerned with variables in which the process function is essentially static. The method has occasionally been used in a different situation. iYhen there are uncontrollable variables, such as catalyst activity, whose levels si.i:adily change with time, evolutionary operation may be used to follow the moving optimum. I n this use it provides a form of manual adaptive optimization. Considerable experience exists in the successful application of evolutionary operation (70, 77). This procedure, therefore, provides a good starting point from which to consider methods of adaptive optimization of continuous processes. Essentially, two new factors arise: (1) the continuous nature in time of the deliberate changes and disturbances and (2) the necessity to consider the dynamic behavior of the process. Discrete Changes of Process Variables. No Dynamic Effects. I n Figure 2, suppose the process function for the response, 7 , with respect to the process variable x does not change with time and evolutionary operation is being conducted for the variable, x . T h e plant is run alternately at two levels, .xu and x b , and the averaged differences in observed responses at these two levels are used to indicate the direction in which to move. The true response, 7, will, in general, never be observed because of random fluctuations or errors inherently present in the ‘observations. The observed response! y, of q always includes error e in accordance \vith the relation JJ(.W)= q ( x ) e. Let yo be an observed response at level xu and yb an observed response at level xh. The difference bet\veen two averages of runs made at conditions x a and xb is
+
.yb - P(,
= 76 - 7 a $. AZ
(1)
where it is assumed that A 2 is independent of the true difference in response q b - va. T o ensure that the estimate of y b - qa is as reliable as possible under the experimental conditions, the averaged error AZ must be as small as possible for a given amount of repetition. Suppose a total of N runs are made, N/’2 runs at level xu and z\r/2 runs at level X b . There are various ways in which such a
sequence could be performed-for example, a t one extreme, all the xu runs in sequence followed by all the x b runs in sequence. At the other extreme an xu run might be followed by an xb run alternately. Both are specific examples of the more general situation in which n successive runs at level x, are followed by n successive runs at level x b followed by n successive runs at level x,, and so on. If A T = 2m, then r measures the number of alternations of runs employed at level xu and at level x b . Then in the first arrangement mentioned, r = 1, n = N 1 2 ; in the second arrangement 7 = Ar/2, n = 1. Various arrangements between these two extremes exist. iVith N fixed, consider the problem of how the value of n should be chosen so that the variance of AB is a minimum. This depends on the nature and characteristics of the error, e. The error committed in the zth run would often not be independent of the error committed in previous runs. Assume that the time sequence of errors follows a stationary process. One method of characterizing the relationship between successive errors is to consider the correlation coefficients between errors one run apart, two runs apart, three runs apart, and so on These are called serial or lag correlation coefficients and are denoted by p l . Thus pt = E [ ( y , - vt)(yt-t - v , + i ) l / o ~ = E[eLX t t + t l / u ;
(2)
where y a is the observed response in the zth run and is the value that would be obtained if there were no error. Similarly y a f Z is the observed value 1 runs later and q l T l is the corresponding true value. The “true” values referred to are the hypothetical mean or expected valueq if infinite repetition were possible at exactly the same conditions. The symbol E denotes expected or mean value. c: is the variance of the observations and on our assumption is constant for all z and all 1. I t is defined by u: =
E [ y z - 1 J 2 = E[e:l
= E[y,+i
- vr+tl2 =
(3)
E[e,’+t ]
If the errors are completely random-that is, the value of ei is in no way connected with the value of ei+ I-then pl = 0 for 1 = 1, 2, 3, The variance of ( g b - Qa) is then completely independent of the choice of n and r. Here are two extreme situations:
.. ..
Errors of runs close together are strongly positively correlated (Figure 3,C). I n this case a plot of errors typically shows slow trends” as in Figure 3,A. Thus changes should be made
:i)q PIDSCSI funclion
2 k-
5:
I
1 .
1
.e
PIOCCISv ~ r l a b l c I Flow r a t e
i
Figure 1. Response as a function of flow rate and catalyst activity
Figure 2. Effect discrete changes process variable VOL. 1
NO. 1
of of
FEBRUARY 1 9 6 2
3
Y
I
In a similar manner, if N = 20 and
A
i.
1
= -0.4,
then
- ga) = 0.352 0:.
% 6. Carry over
A . Trend
corrrlalion
In this case design ii yields d variance which is about five times less than that obtained using design i. Thus, even for a moderate amount of correlation, the arrangement or design of the scheme is of considerable importance. In many practical situations the observations are such that a slow trend is often intermixed with a carry-over situation. I n such an event the value of n will be a compromise. A detailed analysis of the strategy involved, when the more complete correlation pattern is taken into account, is described elsewhere (9). Only in very rare circumstances will there be no correlation between the errors in successive observations.
Correlation
I*
I
If n = 1, V(W
pl
C , Correlogram
D. Correlogram
Figure 3. Characterization of noise in the discrete case
Adaptive Optimization
from one level to the other as often as possible, so that the differences y b - ya are associated with small error changes. Errors of runs close together are strongly negatively correlated (Figure 3 9 ) . This arises when the product of a process is physically carried over from one batch to the next because of incomplete emptying of vessels, pumps, and lines. This plot of errors typically shows rapid fluctuations about the true value of the response. If all error is due to carry-over effects, one might expect high and low results to occur alternately (Figure 3,B). In this situation a series of batches should be run at one level before changing to the other. T o illustrate this mathematically, suppose that the dominant correlation is that between successive observations-that is, relative to p l the higher lag correlations, p2, pa, .,can be neglected. Let us refer to n batches performed in sequence at condition x, followed by n batches in sequence at condition xb as an alternation. Suppose r repetitions are made of this alternation so that N = 2 m and furthermore let the response of the ith member of j t h alternation under condition x, be denoted by yacl. Similarlyyat, is the response of the ith member of the jth alternation under condition x b . Hence
..
(4)
Then from the well known expression for the variance of a linear function, we obtain
Continuous Change of Process Variable. No Dynamic Effects. In practice evolutionary operation is run on continuous chemical processes by treating the process as if it were producing separate batches of product. After each change in the levels of the variables has been made, a period of time is allowed for the process to settle down before measuring performance. I n this way the necessity for considering explicitly the dynamic characteristics of the process is avoided. However, when a permanent system of continuous adaptive optimization is planned, a more refined technique and careful analysis involving dynamic behavior are needed. To make this analysis the consequences of introducing continuous cosinusoidal changes in the level of a single process variable are considered. At first the frequency of the cosinusoidal input signal is supposed to be very low, so that frequency response effects may be ignored. The output response then contains periodic terms of the same frequency as that of the input (fundamental) together with those of the second and higher harmonics. Specifically, let the true unperturbed steady-state response, 7,be a function only of the level, x , of a particular process variable when all other conditions are maintained constant. Suppose x is centered at some particular value, xo, and in the neighborhood of X O , the process function can be represented by an equation of the second degree of the form
Let a cosinusoidal variation of amplitude 6 x and frequency w be imposed on X O , so that x = xo
I t is clear that if p l is positive, the variance will be minimized by making n as small as possible-that is, n = 1. When pi is negative, the variance is minimized by making n as large as possible-that is, n = N / 2 . For example, if N = 20 and p l = +0.4, then i.
If n = 1, V(gb - g,) = 0.048 u:.
ii.
At the other extreme ifn =
N -, 2
V(&. - ga)
=
0.336
gi.
Hence, using design i the variance is seven times smaller than that obtained using design ii. Therefore seven times as many observations are required using design ii to attain the precision obtained using design i. 4
l&EC FUNDAMENTALS
+ 6x cos
wt
(7)
Since frequency w is sufficiently low for the frequency response characteristics to be neglected, the perturbed steady-state response is
T h e cosinusoidal forcing function thus produces a n output containing a fundamental whose amplitude is proportional to Pl--i,e., the slope of the process function at x = X O . This fundamental will be in phase with the forcing function when the slope is positive at ro and 180’ out of phase when the slope is negative at xo. T h e curvature of the process function as represented by will give rise to a harmonic at twice the frequency of the forcing function. Thus, as shown in Figure 1, if the current catalyst activity was c2 and the flow rate was maintained at x 2 , there would be no signal at the fundamental frequency induced in the output when a cosinusoidal *ariation was imposed in the flow rate about x2. However, if the catalyst activity now changed to c1, a signal at the fundamental frequency would be induced in the output which would be 180’ out of phase with the input. Similarly if the catalyst activity was c3, there ivould be a fundamental in the output wihich is in phase with the input. Therefore, if the frequencies are low, the estimation of P i can provide a basis for a process of adaptive optimization. The level of the process variable, .Y, can always be adjusted so that there is no fundamental in the output corresponding to x = x i when c = c1, to x = x 2 when c = c2, and to x = x 3 , when c = c3. T h e direction of change in x is determined by the sign of Pi. Effect of Noise, No Dynamic Effects. The output signal is, in practice, more 01- less obscured by noise generated in the system. The effect of this noise is now considered. Let it be assumed that at time t , the observed output, y , can be represented by the additive setup y d x ) = vt(x)
+
Therefore, in order to estimate b l with as great a degree of precision as possible, it is necessary to find those conditions which maximize this ratio. Now for a given value of 7, et cos ut dt
1’
is a measure of the power spectral density
(I), P(u). of the noise, e, at frequency w . Consequently, on this simple theory for which the frequency response characteristics of the system are ignored, the optimal frequency to be used is that frequency a t which the power spectral density is the least. Continuous Change of Process Variable with Dynamic Effects. So far the discussion has been limited to cases where the forcing cosine wave was of sufficiently low frequency that dynamic effects could be ignored. At higher frequencies the frequency response characteristics of the system have to be considered. Specifically, as the frequency was increased, the output at the fundamental frequency would undergo a phase shift and would be subject to attenuation-that is, a reduction in amplitude. Other things being equal, the system of optimization discussed earlier can still be used. Three methods come to mind. METHODI. If at frequency w the frequency response characteristics are described by the attenuation, y, and phase shift. p, then that part of the output signal at the fundamental frequency is yt
= vp18x cos (ut
+ p) +
If a signal is now introduced, and if over the region xo - 6 x to xo 6 x an equation of second degree in x is adequate to rep, observed output will be resent the process function ~ ( x ) the yc =
(PO
+3
+
-t8x81 cos wt
c ~ ~ i i )
1
+ et
~il8x2cos2wt
p cos at
-
+
ct
(1 5)
Now two correlator signals are introduced zt =
cos w t
and
+
1
cos
= vP18x
v P 1 8 x sin p sin w t
(9)
et
et
w t = sin w t
and the quantities
(10)
Following standard statistical theory, the least squares estimate of is provided by introducing a cosine wave whose value at time t is z t == cos ut which will be referred to as a correlator signal. Then the quantity
and
are calculated.
Now cos
E(m1) =
vP1
E(m2) =
-v/31
p
and 1 =
61
+ 4bl)
6.~~81 cos 1 2wt
+ e t ) cos
wt
dt
(11)
is calculated, where e ( hl) =
(a
an integral number of periods. variance of b 1 are supplied by
6x7
) -’
Le,
cos wt dt and 7 is
Then the expected value and
E(bi) = 81
= u*
(19)
and r n ~and m2 are uncorrelated. For the estimate of VPI the ml can be used. Such an estimate quantity M I = d m y appears to have an advantage, since it can be calculated without knowing p. METHOD11. If cp is known, the least squares estimate, MII, of up1 may be calculated where
+
A411
(12)
(18)
sin p
var ( m l ) = var ( m t )
=
ml cos
p
+ m2 sin
p
(20)
which has a variance, 2. METHOD111. Alternatively, if p is known, we may use the single correlator signal z’t =
cos ( a t
O n the average, 61 will be equally distinguishable from the noise contribution for a’llvalues for which the ratio
when an estimate,
is the same.
This estimate too has a variance,
M111,
of
VOL. 1
up1
+
p)
(21)
can be obtained directly from
$.
NO. 1
FEBRUARY 1962
5
I n practice it appears that only the last method is available, since in actual dynamic systems of interest the simple theory outlined above does not immediately apply. This is illustrated later in the case of simple consecutive reactions. The simple theory so far discussed may be adapted, however, for use in such systems, provided Method I11 is used to estimate the process function gradient. A Simple Reaction System. To introduce dynamics specifically, a simple type of reaction system is now considered in which a reactant N1 decomposes to N2, which in turn decomposes to N3. The reaction is supposed to be conducted in a n isothermal, ideally stirred tank reactor into which N 1 is continuously fed (Figure 4). For this system the yield of Nz passes through a well defined maximum as the flow rate is increased. If the specific rate constants are k1 and k2, the system may be represented by Ni
ki
nomic criterion such as profitability calculated from yield and the process conditions would usually be considered. Mathematical Analysis. Differential Equations. As shown in Figure 4, let the instantaneous outlet concentrations of Ni, ND:and N3 be 71, 72, and 73 moles per liter, respectively. Let only S1 at a concentration of 1 mole per liter be fed into the reactor at a flow rate of q liters per minute. The stoichi72 73 = 1. If ometry of the system is represented by 91 the volume of the reactor is V liters, the simultaneous differential equations describing the reaction system are
+ +
ki
+ N2 +
Na
v
where P = x min.-’
For this study other reactants which may be necessary for N 1 and NSto decompose are supposed present in large excesstheir effect will enter as a multiplicative factor only with k l and k2. I n the first analysis, it is assumed that the instantaneous rates of decomposition of N 1 and Nz, respectively, follow firstorder kinetics-that is, the rates are proportional to the first powers of the instantaneous concentrations of K1 and Nz, respectively. Later the treatment is applied to more complex sequential reactions of higher order kinetics. The importance of this type of system is that it might be expected to approximate certain situations found in the chemical process industry where the response of interest actually passes through a maximum value. Typically, a “product of interest” is formed which then is decomposed to some “undesirable product.” I n practice not yield itself but some eco-
=
space velocity. For convenience it
is assumed that V = 1 liter. LAPLACE TRASSFORM. With the Laplace transform defined as
and remembering that
where q I , 0 refers to the initial conditions, we obtain the linear equations 71
112
?3
p+x+ki
= p111,o
-kl
.b+x+kz
- k2
p
+
+
x
= P112,?
=
x
(24)
POW
UNPERTURBED STEADY-STATE Sor UTION. The steady-state solutions at t = m for the specific flow rate x = xo are obtained by substituting@ = 0 in Equations 24. The solutions are
+ kl)
41,s
= xo/(xo
02.8
= klXo/(Xo
7 ~ 3 = , ~
+ kl)(xo + kz) +
k l k 2 / ( ~ 0 kl)(xo
+ k2)
=
kltll,s/(xo f k z )
(25)
where subscripts refers to steady state. SOLUTION WHEN 7 2 , 8 Is AT ITS MAXIMUM. O n setting b72,,/bxo equal to zero, we obtain the maximum value of q 2 . s occurring at x: = d k l k z
(26)
when
e Y , 0 1 concentrotion
n,, inlet N
:I
mol liter“
I
71t
Figure 4. reactor 6
It+
N,
L Na
Isothermal
1 s ‘71.
inlet‘l’
ideally
l&EC FUNDAMENTALS
=
dL/(dh+ 4
77:,9
=
kl/(dh
v ; , ~=
It
?I
?I:..
i
stirred
)
+ dW d k l k z / ( d & + dk;)’
GENERAL EXPRESSION FOR GRADIENT OF From Equation 25
Reoction
TO x o . Concrntrat ion Stoichiometry
continuous
2
flow
q 2 , s WITH
(27)
RESPECT
GJis the Laplace transform of the displacement in T~ from the steady-state value. Also 6x is the Laplace transform of the perturbation in x , since
LAPLACE TRANSFORM FOR SYSTEM INITIALLY I N UNPERTURBED STEADYSTATE. Suppose that “initially” the unperturbed state situation already exists at the flow rate x = X O , so that the subsequent behavior of the system for perturbations in the flow rate is described in Equation 24 with n,,o set equal to 7, ,,-that is,
-71
52
)+x+kz
fi
- kz
Ps1.a
=
Ps23s
(29)
en.
L ( X 0
+ 6x1 - L j x o )
6x
= .bt/a3s
+ X
67;
G = - = A-I(u - n,)
Then Equation 29 can be written in an obvious matrix notation as ,4ii =
=
Consequently, if G is the vector of transfer functions for 71, 72, and qs-that is, a vector whose elements are the ratios of the transform of the output displacement to the transform of the input displacement-we have, explicitly,
+x
=
+ 6x1 - xn]
+ 6x) - xn=
= (xn
53
$+x+ki
-ki
L j 6 X ) = L((x0
+ xu
(35)
6.7
I t can be readily shown that A-I is given by
(30) 1
+ +k i ) ki (P + + k i ) ( $ + + k2) (P
A-1 =
x
x
where the matrix A
=
p
[
+ x +- k ki i
$+x+k2
-kz
the vector ii
= (;I,
p f x
1
(P
x
-
1
;iz, ; i d )
dii g = - = An-’(u dx
Limit A-l(u - n,) = limit (Ao
If the flow rate is maintained at x = XO, then the unperturbed steady-state conditions exist and the solution of Equations 29 and 30 is iiz.8
= 72.8;
ii3.a
- n,)
(37)
since
the vector u ’ = (1, 0, 0)
7i,5;
x
Kow if 6x + 0, then Z / 6 x d i / d x and the corresponding instantaneous transfer functions by the vector g = (gl,gz,93) when
the vector n’, = (si,8,mS, ?s,~)
ii1.8 =
+ + kz)
+ GxI)-’(u - n,) = Ao-I(u - n,) (38)
6S+O
Thus
(31)
13.8
where the situation can be expressed by the relationship Ann, = pn,
+ xou
(32) (1
here An is the matrix A with xo replacing x. TRANSFER FUNCTIONFOR DISPLACEMEh r FROM STEADY STATE AS A RESULT O F CHANGING FROM x = xo TO x = x o +ax. After conditions have been established at the steady state with x = X O . the Laplace transform for q, changes from yj,, to v j , , 4- Gj. Then Equation 29 becomcs M
p
+ijld ++ST +- k i xa
6x
T2.S
kl
{J
+ ;6;
+ 6K fi + + 93,.
+ + 6 x +- kakz xn
xo
=
psi,,
= psz,. ax = p s 3 , s
which can be written in matrix notation as (Ao
+ 6xI)(n, + 8 7 )
=
pn,
+ xou +
6xu
where I is the identity matrix. On subtracting Eqluation 32 from Equation 34 we cxplicitly 2=p - n,) 6X
(34)
Obtain
- tl1,s)khz - qz.akz(P + x o + ki) -
+ + k2XP + xo -t k i ) + k J P + xo + kZ)(P + xo
~a.s(P
(P + xn
=
If we p u t p
=
0, the correct value of the gradient is obtained-
sj(p = 0 )
+ + 6.x 20
xu)
= gL0) =
btdh
(40)
(33)
Response to Input Cosine Wave. To obtain the response at the fundamental frequency to a cosinusoidal input of frequency w and amplitude ax, p is replaced in gl(p), gz(i;), and g&) by iw to obtain g l ( i w ) , g z ( i w ) , and g3(iw), Limiting attention to the intermediate product, Nz, whose steady-state process function shows a maximum with respect to x: we have gdiw) =
ki(1
(iw
+ ki) + xn + k i ) ( i w + + kz) 7 ~ 3 )
tlz..iw
t/z.s(~n
xn
(41 1
- , d o ) - iw7z..a1az
+ iwa1)(l + iwaz) where al = l/(xo + kl), a2 = l/(xo + k z ) , and gz(0) is the slope (1
where the vector (u - ns)’ = (1 - v,,,, and the matrix A = A. 6x1. Since = L { q l ( x o i- 6x1 - L ( vj(x0)
+
Fj
= L(vi(x0
+
6x1
-
-~
-
2 , ~ ,q3,J
of the process function.
]
qj(xo)} = L ( 6 V i j
Thus g2(iw) = Re’q = R cos VOL.
1
‘p
+ iR sin
NO. 1
(42)
‘p
FEBRUARY
1962
7
where R is the “gain” or amplitude ratio and p is the phase change of the fundamental in the output and both R and p a r e functions of xo and W . Now for the input x(t) =
xg
+ 6x cos wt
(43)
the output 7 2 at the fundamental frequency is nz(t) = =
72,s ~
2
+ 6xR
COS ( a t
+ 6 x ( A cos ,
~
wt
(44)
- B sin w t )
Thus coefficients A and B are the real and imaginary parts of and the gZ(iw). Now when vr,*is at its maximum value corresponding value of xo = x:, the transfer function is
,)I:(
-
Except a t zero frequency ( w = 0), g l ( i w ) will not be zero. I t will in general contain both real and imaginary parts. Consequently, even at the maximum, at x:, where the slope of the process function is zero, there will be a transfer of signal A* cos ut - B* sin w t of amplitude R* and phase change cp**
Analysis in Terms of Phase-Shifted Components. Now instead of analyzing the output in terms of wt we analyze in terms of the phase-shifted function (ut e). Let us introduce a phase shift, 8, such that
+
and (47)
I n terms of the phase-shifted signal, the output a t the fundamental frequency is 7 * ( t )=
+ 6 x [ A ’ cos
(wt
+ e ) - B’sin
(wt
+ e)]
(48)
+
Wy
W t / ~ , ~ Q i ( l ~ (W l 2Q~)-1’2(1
+
(51)
W2Q;)-112
Thus the output at the fundamental frequency may be written in the form m(t) =
+ V)
where A = R cos p and B = R sin p so that R = (AZf B2)ll2 and p = tan-’(B/A).
wherep* = tan-‘
and
~
2
+ 6x[vg*(O) ,
~
+ +
cos ( ( ~ t 0)
WY
sin ( w t
+ 011
(52)
I n practice this means that by using a correlator signal fit
=
COS ( w t
+ e)
(53,
we can estimate g ~ ( 0 )= d?2,s/dxg = P I , the slope of the process function, from the coefficient of the phase-shifted cosine term in Equation 52. This component is the useful component and the phase-shifted sine term is the useless component. A practical difficulty is that the appropriate phase shift, 8 , is itself a function of X O , so that strictly speaking the phase shift to be introduced will change at different levels of X O . I n practice, however, it wou!d appear satisfactory to use a constant phase shift, e*, the value appropriate at 2,” corresponding to maximum ~ 2 , ~ I .n the next section a number of numerical values indicate the practicality of this approach. At x:, the coefficient of the phase-shifted cosine component is zero, so that
- B’* sin (ut + O * ) = R* cos ( w i from which it follows that e* = q* - 90’.
+ v*)
I n any actual industrial process, the exact mechanism of the process and the values of the constants may not be known. Consequently it would generally be necessary to calibrate the system by experiment and, in particular, to determine a compromise phase shift. I t is believed that the phase shift appropriate for the “average” optimum conditions would usually be found satisfactory. Numerical Investigation. Figure 5 shows the unperturbed steady-state process function
for k l = 0.1 min.-’ and k2 = 0.05 min.-l. maximum in
The value of the
where A’ and B’ are the real and imaginary parts of g 2 ( i w ) e i e . NO\V
= vgl(0)
- iwy
where Y
=
(1
+ w2a:)-’*(1 + w2a;)-li2
=
attenuation oftheslope
(50)
occurring at x i = dk,O = 0.070711 min.-l. For reference x,* and the two values x ’ ~ = 0.03300 min.-l and Po = 0.15000 min.-l have been marked or indicated. The range x ‘ ~to x ’ ’ ~ includes conditions such that ~ 2 is, within ~ -10% of its maximum value. Figure 6 shows the value of the amplitude of the phaseshifted sine component (useless) as w is changed. I n general the quantity, wy, has a maximum at the value of w = [(YO
0.01
0.02
004
Space
Figure 5. 8
I&EC
0.1
0.06
velocity
x0
0.2
mi”-’
Steady-state process function
FUNDAMENTALS
0.4
+ kl)(XO + k d P 2
obtained by differentiating Equation 51 with respect to w This frequency may be called the resonance frequency for the response, 72, with respect to the process variable, x . Figure 7 shows the attenuation, Y , of the slope of the process function occurring in the phase-shifted cosine component (useful) as w is increased. As the resonance frequency is approached, a marked attenuation occurs, the amplitude of this useful signal having been reduced almost exactly by a factor of 0.5. By the time that a frequency of about ten times the resonance frequency has been reached, this amplitude has been reduced by a factor of -0.1.
The decision as to the best frequency a t which to operate may be reached by considering the signal to noise ratio. Specifically, for a fixed input amplitude, ax, we should max) or, equivalently, the difference imize the ratio ~ z ( w /P(co) 2 In ~ ( u) In P ( w ) where P ( w ) is the power spectrum. Figure 8 shows the optimal phase shift, 8, over the range do to xB0 calculated from Equation 47 for representative values of w. Little is lost by employing a fixed phase shift, e*, appropriate to maximum value of ~ 2 over , ~ the range doto x”,, once w has been fixed. Although, in practice, an experimental approach would be adopted to determine the best phase shift, it is perhaps helpful to have a mechanistic picture of each kind of phenomenon which may be responsible for the shift of the maximum. Thus, for example, let us suppose that the object of adaptive optimization of a reactor such as the one discussed in this section is to maintain the maximum concentration of species Nz in the ki
Figure
6.
Amplitude of useless component
k2
catalytic reaction S1 -+ Nz N3, where catalyst decay is occurring with respect to time. Typically, the change in the properties of the catalyst can be represented by a change in time in the effective rate constants, kl and kz, which might follow a course like that shown in Figure 9. I n such a system, a compromise phase !shift corresponding to some “average” situation in the period iof interest would be used. -f
Extension of Theoretical Model
Multicomponent First-Order Sequential Reactions. This type of analysis can readily be extended for any particular ki kr
ks
ki
component in the multicomponent scheme N1+ Nz
-
+
NB-+
ks
, etc. For the moment the assumption that the specific rate constants refer to first-order kinetics is retained. If our object is to optimize 73, which shows a maximum with respect to x o in the multicomponent system, the differential mass balance equation Ibr Q Q N4 + T\;S
Frequency
Figure 7.
o , min-’
Attenuation of useful component
must be considered simultaneously with those for 71 and 72. Analysis in Terms of Phase-Shifted Components. T h e output at the fundamental frequency in terms of the phaseshifted time function, w t 8, is easily shown to be v3(t)=
VJ
+
+
Gx(Y[gJ(o)
+
q38
~
*
~
l COS ~ 2
+ e )+
(ut ~ j
wy
sin
~
(ut
+ e)}
(55)
+ 5).
where u, = 1/(q The cosine component contains not only g3(0), the slope of the process function, but also a term in w2. Unlike the situation for Q?, it does not appear possible to isolate completely a component which is directly proportional to the gradient. Possibly if this extra term is troublesome it could be corrected by suitable calibration. Numerical 1nvestiga.tion. Figure 10 shows the unperturbed steady-state process function
for k l = 0.1 min. -I, k2 = 0.05 min.-l: and k3 = 0.03 min. -I. As before, we have indicated x i = 0.025 min.-’ and the two values x f 0 = 0.014 min.-’ and x ” ~ = 0.045 min.-’, so that ~ -10% this range includes conditions for which ~ 3 is , within of its maximum value. Figure 11 shows how the amplitude of the useless phaseshifted sine component changes with w. Here again this amplitude has a maximum for each value of x o at a particular
Figure 8.
Optimal phase shift
Fig. 9. Typical change constants with time VOL.
1
NO. 1
of
effective
rate
FEBRUARY 1962
9
SPOC* velocity
,
io mi"*!
Fieguency
Figure
10. Steady-state process function
Figure 1 1 ,
value of w which is termed resonance frequency for 7 3 with respect to the process variable, x . Figure 12 indicates the attenuation, V , of the useful phaseshifted cosine component. As the resonance frequency is approached, v changes rapidly and at the resonance frequency attenuation of the useful signal is reduced almost exactly to 0.5. As w increases further, the attenuation asymptotically approaches zero. The best frequency at which to operate for a given value of 6x is that at which the signal of interest to ) lnP(w), is a maximum. noise ratio, measured as 2 In ~ ( w Figure 13 shows how the amplitude of the cosine component-namely, gS(0) w2?3,,ala2a3, plotted as function of xo-changes as the frequency is increased. \\'hen w 0, this is simply the slope of the process function. As w increases, an appropriate correction must be applied to obtain g3(0). The optimal value of the phase shift, e, for various frequencies in the range to of operating conditions is shown in Figure 14. The indications are that the error will not be great if a fixed phase shift, e*, corresponding to maximum is used for each frequency over the range X ' O to
w , min-'
Amplitude of useless component
where al > 0, so that ai is the order of the reaction. I n the appendix, if we can assume that 6x, the amplitude of the input cosine wave, is sufficiently small so that the perturbations from the steady state are not too large, then the final equations describing the system will be of the same form as those for the system with first-order kinetics but with the kl replaced by Kj = k,aj~,,~4-1
(57)
The phase-shifted transfer function for the j = J t h component is
-
+
+ ....
X"0.
Analysis When Orders of Reactions in Multicomponent System Are Not Necessarily Unity. Let us denote by N, where j = 1, 2, 3, 4, ., n, a particular component of
..
ki
k2
ki
concentration n, in the sequence of reaction NI + N B+ N3 + kn-1
N4e.. . N, occurring in an ideally stirred isothermal reactor. Suppose the rate of reaction of the j t h component is klqlai,
Frequeney
Figure 12. 10
w , mi"-'
Attenuation of useful component
l&EC FUNDAMENTALS
(1
+
WzQ:)1/2
+
where a j = 1/(n K,,J. The observed output response, V J , for the input x ( t ) = x o 4-6x cos w l in terms of the phase-shifted components is ~ ( t =) vJ,3
+ 6 x ( A ' ~cos ( w t + e ) - B; sin
Figure 13. tenuated)
(wt
+ e))
(59)
Amplitude of cosine component (unat-
absolute temperature in degrees Kelvin; both F, and H, are assumed to be independent of temperature. 'The process variable in this system is denoted by Xwhere I
X = 1/T
I
(62)
I n the neighborhood of Xo,for small changes in X-that is, when X - Xo is small compared to Xo-the value of the rate constant at Xis given by
I---:
From Equation 61 I
1
I
I
space veIocity a.,
Figure 14.
I
0 03
0 02
0001
l
l
0 04
mm-'
Simultaneous Differential Equations. The set of simultaneous differential equations representing the instantaneous. situation in the reactor is
Optimal phase shift
where A ' J and B'J are the real and imaginary parts of EJ(iw)e"eJ, so that A'J =
gJ(0)
+ terms involving powers of w when J > 2 J
(1 j=1
+ w2Q,?)'/'
(60)
where d , =~ kl,xo(l Once again for the third and higher components the useful cosine component contains terms associated with w which need to be corrected for. Optimizing with Respect to Temperature. So far the optimization of the yield of an intermediate species, Nj, with respect to the flow rate has been considered. A similar development can be used where the system can be optimized with respect to temperature with all the other process variables, including flow rate: maintained constant. Joint optimization with respect to several process variables is not considered here. Arrangements are macle to introduce a sinusoidal disturbance tll
+ H,Xo)
and x o = qo/V = constant space velocity. Laplace Transform for System Initially in Unperturbed Steady State. We define the Laplace transform as in the previous cases. Let us suppose that steady-state conditions have already been reached at the specific value of X = Xo. Hence initial conditions now refer to steady-state conditions at Xo-that is, ~,,o = qj,$. Equation 65 can therefore be transformed into the set of linear equations 1 2
1 3
= Ptl1.8
+
xo
= Ptl2.a
only in the temperature of the system. Furthermore, at this stage heat transfer and other such phenomena are not included and a t the frequencies used, it is assumed these are unimportant. Mathematical Analysis. For the simple first-order conki
which can be rewritten in matrix notation as Aii = pn,
+ xou
(67)
where
k2
secutive reaction scheme: N1+ Nz + Na occurring in an ideally stirred reactor, a fixed flow rate of Nl is maintained at unit concentration into the reactor. As the temperature of the ideally stirred reactor is varied with respect to time, concentrations 71, q2, and 773 will change because the specific rate constants, kl and kz, vary with temperature and the steady value of q 2 will show a well defined maximum when plotted against the temperature. Now in general
and u ' = (1,0,0). If, however, X is maintained a t X = Xo,then steady-state conditions exist. The situation is described by AOfi, = p n,
(68)
where A, is the matrix A with Xo replacing X and the solution is VI.,
where F, is the so called nonexponential factor, H, is the ratio of the activation energy to the gas constant, and T is the
+ xou
= 7i.a;
-
72.8
=
7l2,a;
? 3 1 ~= V a s a
Unperturbed Steady-State Solution. At the specific value of X = XO, the steady-state responses are VOL.
1
NO.
1
FEBRUARY 1962
11
+ di,o + C1,oXo) + k1,xJ) + ci.oXo)/(Xo + dz,o + C Z , O X O ) vi,.&i,xo/(xo + k z , x o ) v2,Jdz.o + CZ,OXO)/XO vz,&z,xO/Xo (69)
71,s = X O / ( X O
vz,a
= tli.e(di,o
vas.
=
=
XO/(XO
=
=
T h e slope of the steady-state process function, qz,s, with respect to Xo is dtlz,l
d
+ ki,xo) + ki,xo)(xo + kz.xo) --ki,xovi,dfixo + kz.xoqz,Jfn(xo + k i , x o ) + k i , x J ( x o + kz.xC)
~~,SCI.OXO
x =
(XO
=
~Z.SCZ.O(XO
(XO
=
[“’I
where the matrix Co =
-ci,o
pn.
c;,O -c2.0
+ xou
(71)
:]
+ + +
+
cl,ovl.s(iw X O ) - c z . 0 ~ 2 ,( ~ i w +XO ki.xo) (iw xo h . x o ) ( i w xo kz.xo)
=
- gd0)
+ +
6;)
+
(70)
+
+ SXCO)(
