1974, Vol. 29, pp. 1193-1204.
in Great Britain
NONLINEAR STOCHASTIC SIMULATION OF STIRRED TANK REACTORS N. J. RAOt Department of Electrical Engineering D. RAMKRISHNA* Department of Chemical Engineering and J. D. BORWANKER Department of Mathematics, Indian Institute of Technology, Kanpur 208016, India (Received 3 February 1973; accepted 18 September 1973) Abstract-An algorithm for solving nonlinear stochastic differential equations due to the authors has been used for simulating the behaviour of stirred tank reactors subject to random parametric fluctuations. Two examples have been chosen from the literature derived from the works of Pell and Aris, and Aris and Amundson who approached these problems through linearization techniques. The limitations of linearization, qualitatively pointed out by these authors, have been evaluated in quantitative terms in this work. Several interesting predictions reveal that while linearization provides an effective techniaue in some situations. in several other instances, powerful limitations exist necessitating direct analysis of the nonlinear equations. INTRODUCTION Successful operation of industrial processes depends on consistent qualitative and quantitative characteristics of the operating parameters in-
volved in these processes. Frequently, such operating parameters are susceptible to random statistical fluctuations which may or may not significantly affect the process. An a priori assessment of the requirement of control of any specific parameter hinges on one’s ability to evaluate the effect of random fluctuations in that parameter on the process performance. Alternatively, the deleterious effects of fluctuations of a parameter may be overcome by controlling an interacting variable. In either case, it is clearly essential that prior to control, one has information about the relation of process performance to fluctuations in process parameters. Such information necessitates the solution of equations which result from appropriately modeling the process with the incorporation of statistical fluctuations in the pertinent parameters. Solution of such stochastic (random) equations has +Present address: School of Automation, Indian Institute of Science, Bangalore 12, India. STo whom correspondence may be addressed.
been an exceedingly difficult task save for certain special situations in which the noise is idealized and the model equations are linear. Until recently, nonlinear equations have resisted attempts at numerical solution, and had forced earlier workers to adopt various linearization techniques with their obvious limitations. Rao et al. [I] have presented an algorithmic solution of nonlinear stochastic differential equations which effectively provides for simulation of non-linear processes subject to random parametric fluctuations. The objective of this paper is to investigate the behaviour of stirred tank reactors under fluctuations of selected reactor parameters. Further, it will probe into the merits of linearization which has formed the basis of past work. Specifically, we shall be concerned with two examples of reactor simulation. The first of these owes its origin to the work of Pell and Aris. They analyzed the effect of random temperature fluctuations on the yield of a useful reaction intermediate from a sequential set in a perfectly mixed continuous stirred tank reactor. The second example is obtained from Aris and Amundson who considered the effect of fluctuations in flow rate and feed temperature on conversion in a CSTR entertaining a second order reac-
J. RAO et al.
tion. The assumptions, on which the cited works are based, are relaxed in the present work, so that the full implications of the nonlinearity in the Arrhenius reaction rate expression are brought forth. As a prelude, we shall present briefly some aspects of stochastic modeling and the evolution of an algorithm developed by Rao et a/.[ I] for the solution of stochastic differential equations.
conclusive discussions[S] on the relative merits of the two models represented by Eqs. (1) and (5) in regard to their applicability to physical problems. One may be tempted to lean on Eq. (5) for physical applications because its mathematical basis is sound, although Eq. (1) has been assigned a meaning through a somewhat devious reasoning by Wong and Zakai. The net result is that Eq. (1) has an ‘Ito form’ given by
problem on hand is the solution of the differential equation The
x(t) = x(td +
f tacds), I kl
where t(t) is the familiar characterized by
+ TN= S(T)
b(x, s)dW(s) (5)
in the form
dx = a(x(t),
t) dt + b(x, t) dW(t)
where dW(t) is a Brownian process increment the interval dt, with the properties E[d W(t)] = 0 E[dW(t)2]
Further, the integrals in Eq. (5) are stochastic integrals for a background of which we refer to Wong[41. The second integral in Eq. (5) is due to Ito and is called the Ito integral . Since, informally, d W(t) = t(t) dt
h(x(s), s) dW(s)
or dx(t) = [a(~, t) +fb(x,
where 6 is the Dirac delta function. An alternative stochastic model to Eq. (l), which represents the white noise model, is the Ito integral model x(t) = x(t,)+
EK(r)l = 0 EMr)5(t
x(to) = x0
it would appear that Eqs. (1) and (5) are identical. This, however, is erroneous and has led to some in-
+ b(x, t) d W(t).
b(x, t)] dt (11)
It is interesting to note that if in Eq. (5) one interprets the Ito integral term to be an integral in the sense of Stratanovich then the solution processes for Eqs. (1) and (5) are identical. The Stratonovich integral differs from the Ito integral only in the constitution of the approximating sum. Such features of stochastic calculus have no deterministic counterparts. However, in this paper, we will deal with the Ito interpretation of the integrals concerned. While there appears to be no rational a priori basis at the present stage for a relative assessment of the white noise model and the Ito model, the computations presented for the examples herein appear to favor the Ito model. However, one interesting aspect that deserves mention is that from a practical standpoint, neither the white noise process nor the Brownian process are truly realized. As a result, the noise term t(t) in Eq. (1) must be replaced by a more realistic stochastic process, possibly a ‘band limited noise’ n(t) obtainable from dv dt = - Kq + K’.$(t)
where K and K’ are constants. In such a case the ‘ambiguity’ between Eq. (1) and Eq. (5) disappears and both would yield the same solution process. Notice further that the differential equation
Nonlinear stochastic simulation
possesses the same solution process x(t), as its Ito form since b is independent of x. We shall be concerned with the solution of Eq. (5). The method should be useful for solving Eq. (1) through solution of its Ito form (10). The solution process of an It0 integral equation of the type (5) under fairly general conditions can be shown to be a Markov process (with continuous sample paths)  with a probability density f satisfying the Fokker-Planck equation
$f(x, tl x3,
t)f(x, tlxn s)l
subject to the initial condition f(x, SIX,,s) = 6(x -x,1.
A complete characterization of x(t), the solution of Eq. (5), may be obtained if Eq. (14) can be solved subject to (15). Unfortunately, it is not often possible to solve the Fokker-Planck equation and alternate methods must be found. Rao et al. have derived an algorithm for solving equation (5) by generating discretized sample paths [ 11. The derivation of this algorithm is based on discretizing the time interval [to, t] into a large number N of equal subintervals of length h and rewriting Eq. (5) in the form
a(x(s), s) ds
b(x(s), s) dW(s).
a function of the solution process of equation then it is shown that du =$dti$dx+;$b2d,
which is clearly different from its deterministic counterpart by the last term on the right side. It is readily established by including the term ~(a’u/ax’)dx’ in the expansion and pretending that d Wz = dt so that d W = (dt)“2 from which it follows that d Wdt = 0(dt)“2. Although, this equality of differentials of varying order may sound strange at first, it may be understood from a physical point of view, as the net result of random forward (positive) and backward (negative) motion characteristic of a Brownian particle. It is precisely this aspect that plays an important role in the selection of terms for inclusion in the algorithm for a given order of error arising from the neglected terms. From (17) it is possible to write for the function u(x, t) u(x(t), t) = u(x0, to1
+ t: [4(x(s), s)
u*(x(s), s)a(x(s), s)
which holds for every sample path of the stochastic process x(t). Subscripts have been used in (18) to denote partial derivatives. From (18) a simpler version of the algorithm of Rao et al. [l] can be readily derived which brings forth the essential groundwork for the development of the more complex and accurate algorithm. Thus one may write from (18) u (x(t),
t) = u (xo,
to> + ux(x0, to)b(xo, h)
In (16) the integrands are replaced by their Taylor series expansions about (x., t,) where x. is an approximate value of X(L) and x0 = x (to). The Taylor series expansions should obviously be truncated suitably, indicating the order of the neglected terms. Here too certain aspects of stochastic calculus are encountered that may at first seem ‘idiosyncratic’ to one who is used to ordinary calculus. Reference to Eq. (8) shows that a second order differential is expected to be of the order of the first order differential dt, which lies at the root of some of this deviated behavior. Thus if u(x, t) is
+ o(lt - toI1
in which the stochastic integral j’:, d W(s) may be evaluated [ 1l] as a normal random variable at time t, say Z,, with expectation zero and variance (t to). Applying (19) to the functions a and b in (16) one gets x,+,=x.+a.h+b.Z,.+b..b,~:.+‘dW(t)~:
+ o(h) where b,. represents
the value of b, at x., t,.
that generates another stochastic integral in the last
N. J. RAO et al.
term on the right which is found [ 1l] to be given by
j-;*‘dWt) j-1d W(s)
= f(.Z:, - h)
where 2,. = J::*I d W(t) has mean zero and variance h. Thus Eq. (20) transforms to x.+1
x. + a,h + b,Z,, + b,,b,f(Z:,-
h) + o(h) (22)
Algorithm (22) is easily implemented on a computer since the normal random variable Z,,, is readily generated. Thus discretized sample paths are obtained from which the moments of the solution x(t) are calculated. Algorithm (22) represents an interesting deviation from the corresponding one for the deterministic case. If the error has to be minimized to o (h2) a similar development results in a proliferation of terms consisting of several stochastic integrals representing random variables that are either normal or nearly normal[l, 111. We refer to 111for a complete statement of this algorithm which is omitted here since it requires a substantial investment of space besides leading to an explosive growth of the nomenclature. The situations dealt with here require the extension of the above algorithm to the vector case with scalar noise. This extension is not particularly difficult even if the final algorithm is considerably augmented in length. EXAMPLE
tion of the resulting problem Pell and Aris expressed the rate constants k, and kz in the form ki(T(t)) = &(a’)+ k;(t),
where &(u2) was obtained by assuming T’ to be normal as given by Pell and Aris ; K is a perturbation from the mean, as a result of temperature fluctuation, with the properties E[k#)]
+ T)] = 20&9)6(7).
Using (2.5) Eqs. (23) and (24) may now be written as
+ (Plcl - PzCB)5(t)
which represent the linearized white noise model for the system. In Eqs. (26) and (27) we have set pi = m. Aryaratnam and Graefe [ lo] derived the Fokker-Planck equation for the density of the solution processes of such systems from which equations could be derived for their moments. Equations (26) and (27) submit to a numerical solution through their representations in the Ito forms by the algorithm presented herein. It is also possible to consider the following linear Ito model for the system.
Pell and Aris  analyzed the sequential reaction A&Bk’-C in a continuous stirred tank reactor fed by a stream containing A at a concentration C,. The continuity equations for C, and Cs take the form dCa _=-
- [++ k,(T)]G.
The temperature in the reactor was assumed to fluctuate about a mean value with variance a2 the effect of which was to be ascertained on the concentrations C, and Cs within the reactor. Since a direct substitution of the nature of temperature fluctuations into Eqs. (23) and (24) then precluded a solu-
- [; +
,I(PICA - pzG)
Cl?(t) = c&(h) +
Note that this model is different from the white noise model. Equations (28) and (29) could also be solved by the algorithm presented. Since the linear models were obtained by linearization it would be of interest to ascertain the consequences of the nonlinear model in which Eq. (25) is ignored and the temperature fluctuation T’ is considered as a band-limited noise satisfying the
Nonlinear stochastic simulation differential
!g =- KT’
whether this is the white noise or the Ito model no longer exists. Of course, the real model for the situation considered here, shoald also account for the energy balance seeking the ‘source’ of the temperature fluctuation. This indeed is a tractable problem with the present algorithm, which is however omitted for this example since the second example encompasses this wider perspective.
whose Ito form is T’(t) = -
KT’(s) ds +
= G(t”) +
Computations were made using the present algorithm for various values of u’, the variance of the temperature fluctuation, on the linear white noise model (which is a repetition of the computations of Pell and Aris since identical numerical values of the parameters have been used in this work) the linear Ito model and the nonlinear model. These results are displayed in Fig. 1, and in Table 1. The algorithmic solution of the linear white noise model is identical to the results of Pell and Aris. The steady state concentrations of A and B drop with increasing temperature fluctuations. This trend is understandable since one would expect a preponderance of positive fluctuations in the rate constants about their respective values at the mean temperature, by the very nature of the Arrhenius expression. The linear white noise model tends to show substantial variation with cr2. This is all the more interesting with the nonlinear model exhibiting less of such variation. The linear Ito model however is closer to the results of the nonlinear model than the white noise model, which earns some support to the Ito model at least in the present example. For larger variances the linearization becomes irrelevant and the nonlinear model predicts a smaller drop in the mean concentrations of A and B. In fact the linear
If the value of K is sufficiently large then T’ may be considered approximately white since it has a reasonably flat spectrum over a large frequency band. Equation (30) with a large value of K becomes necessary, since the solution of Eqs. (23) and (24) with a direct substitution of a white noise process for T, does not come under the purview of the algorithm. To complete the nonlinear model, the concentration equations may be written as c,(t)
El Ca ds (32) R(T + T’) )I
(33) where the rate constants have been replaced by their Arrhenius expressions. Equations (31), (32) and (33) completely define the nonlinear model. Note that the ambiguity of
Fig. 1. The effect of temperature
noise on the steady
N. J. RAO
Table 1. Variation of expected values of steady state concentrations functions of U* mes
White noise linear model
0.333 0.326 0.319 0.312
0.333 0,326 0.318 0.310
5 6 7 8 9 10 15 20 30
0.298 0.291 0.285 0.278 0.272 0.267 0.240 * *
0.298 0.298 0.282 0.276 0.270 0.262 0.234 0.208 0.165
Nonlinear Ito model 0.333 0.320
0.305 0.285 0.28 0.262 0.248 0.222
White noise linear model
Linear It0 model
Nonlinear It0 model
0.333 0.323 0.313 0.302
0.333 0.325 0.314 0.304
0.278 0.265 0.25 1 0.236 0.220 0.202 0.077 * *
0.285 0.273 0.263 0.254 0.245 0.234 0.191 0.154 0.099
0.282 0.267 0.252 0.224 0.198
*The value of msS becomes negative, suggesting that the linear whitenoise model is no longer applicable in this range of oz. white noise model shows negative concentrations for large a*. From the foregoing discussion, it is clear that linearization could lead to significant errors even when the random perturbations are not particularly large. Thus for as low a temperature fluctuation of about 4” there is a considerable difference between the predictions of the linear and the non-linear models. EXAMPLE
which may be non-dimensionalized
using the dimensionless
This example is derived from the work of Aris and Amundson who considered a linearized analysis of a continuous stirred tank reactor in which was carried out the second order reaction 2A +products. Their analysis comprised an evaluation of the effects of fluctuations in the feed flow rate, coolant flow rate, the feed temperature and the coolant temperature on the reactor performance. In this paper, we will concern ourselves with only fluctuations in the flow rate, and shall consider both the linearized and the nonlinear equations describing the reactor operation. The mass and energy conservation equations take the form
z = CIC, ;
q = Th/(- AH)C,;
P = VR,IqC, ;
P = 414;
q< = Tch/(- AH)Cf ; Linearized
70 = T,h /(- AH)C, ;
U = VU*/(-
Aris and Amundson equations
2 =- (1 + (u)x -
py + (1 - &)A
(G - C) - 8(kKZexp
(T, - T) +$(-AH)/&
exp (- E/RT)
The fluctuation in the flow rate, A which was taken to be a white noise process by Aris and Amundson, is assumed here to be given by the Ito equation dh(T) = - KA dr +uK
26.3 &2exp (26.3-y)
- et: exp
26.3 - y
67.5 r)S+ Y
time interval of 0.005 for the choice of the parameters listed in Table 3. The choice of 100 for K corresponds to around 16 cycles/residence time while u = O-1 is equivalent to a 10 per cent fluctuation in the flow rate. The computations extended over a period of about five residence times beyond which the accumulated error prohibits fruitful continuation. The steady states about which linearizations were performed are shown in the customary plots of heat generation and removal vs temperature in Figs. 2-5 from which it should be clear that the entertained situations involve multiple steady states. All linearizations have been about the high temperature steady state. In Figs. 6a, 6b, 8a and 8b are represented the expected transient responses for the linear and non-
model is defined by the equations
+ (& +x)* exp [
where u and K are to be suitably found as before. Equations (38), (39) and (40) define the linearized model, which can be solved by the algorithm herein.
(a) 320 -
where we have substituted the values for the Arrhenius reaction rate expression, the heat transfer coefficient, heat transfer area, reactor volume and feed concentration used in the calculations of Aris and Amundson. All other numerical values have also been borrowed from this source. These have been listed in Table 2. Table 3 lists the values of the steady state reactor temperature and concentration for various choices of the coolant temperature and holding time. Equations (41) and (42) are of course to be considered together with Eq. (40).
T C = 450’R 6 225secs q
RESULTS OF COMPUTATION
linear and the nonlinear models were solved by using the algorithm with (+ = 0.1, K = 100 and a
Fig. 2. Heat generation and removal as functions of temperature.
(lb m:;e/ft3)(&) 1
(A) 2.7 x 10”
N. J. RAO etal.
linear models. In Figs. 7a, 7b, 9a and 9b the standard deviations about the expected values are plotted vs time.
TC =4CQ’R 8 I 250 sets
Fig. 3. Heat generation and removal as functions of temperature.
Fig. 4 Heat generation and removal as functions of temperature.
TC = 500”R 0 = 25Osecs
Fig. 5. Heat generation and removal as functions of temperature.
can be captured by examining, say, Figs. 6a, 6b, 7a and 7b. For the operating conditions represented in Fig. 3, consider the predictions of the linearized and the nonlinear models shown in Figs. 6a and 6b. Whereas the linear model shows fluctuations about the high temperature steady state, the nonlinear model predicts a considerable reduction in the mean reactor temperature. Also Figs. 7a and 7b show a high time-dependent standard deviation about the mean behavior for the nonlinear model and a low standard deviation for the linear model. These observations are also true for operating conditions in Fig. 4 as seen from 8a, Sb, and 9a and 9b. Such behavior is certainly in order from a perusal of Fig. 3 which reveals sensitivity of the reactor steady state to small variations in the slope of the heat-dissipation curve. On the other hand, for the operating conditions in Fig. 5, which represents a higher coolant temperature and an identical holding time, the predictions of the linear and nonlinear models are identical. These observations also hold for the operating conditions in Fig. 2, as reflected by Figs. 8a, 8b, 9a and 9b. Thus the behavior of the reactor, for a given fluctuation in flow rate, is seen to depend much on the sensitivity of the steady state. The nonlinear model predicts substantially varying behavior under selected operating conditions. The general
spirit of the predictions
The limitations of linearized analysis of the effects of random disturbances of reactor parameters on reactor performance have been qualitatively recognized by Ark and Amundson and by Pell and Aris. These limitations have been demonstrated quantitatively by solving the nonlinear reactor equations, using an algorithm developed by Rao et al.[ I] for numerical solution of nonlinear stochastic differential equations. Even for ‘small’ perturbations, the linearized analysis is shown to falter if the reactor steady state is ‘sensitive’. The algorithm used herein should also be useful in the solution of stochastic optimal control problems. An interesting aspect of stochastic analysis is the incorporation of noise terms that are more realistic than the idealizations used normally. Thus a fruitful area of investigation is the identification of ‘noise models’ for observed stochastic processes
CS =O.l246lb Non-linear
TC = 500”R 8 =25osecs
Fig. 6a. The change in concentration
as a function of time in the presence of noise.
Fig. 6b. The change in temperature as a function of time in the presence of noise.
0.20 - 1 Non-linear 0.16 _ 2 Linear 3 Non- linear
TC = 4OO”R 0=250secs T S = 809.6’R TC= 500”R 0= 250secs T s = 848.23aR
Fig. 7a. The change in standard deviation in concentration
29 NO. 5-J
as a function of time.
N. J. RAO et al.
Fig. 7b. The change in standard deviation in temperature as a function of time.
TC:4500C 0 = 325 sets TSr790.2’R CS=O.l515lb
TS = 45O”R 8 : 225secs TS = 843.2”R
Fig. 8a. The change in concentration
as a function of time in the presence of noise.
Fig. Sb. The change in temperature as a function of time in the presence of noise.
TC:450°R 8=325secs TS:7902”R CS=O.1515Ib.n
TSz450.R 8=3255erTSz843.2”F cs=o.o775
Fig. 9a. The change in standard deviation in concentration
as a function of time.
60 k 40 -
Fig. 9b. The change in standard deviation in temperature as a function of time.
heat transfer area. A$: species undergoing reaction a, h functioned defined in Eq. (1) C concentration D auto and cross correlation functions defined below Eq. (24) E activation energy f density function for the solution process of Eq. (5) h volumetric heat capacity AH heat of reaction K constant in Eq. (12) constant in Eq. (12) K’ k reaction rate constant defined immediately following P quantities A
Eqs. (28) and (29); also a dimensionless quantity below Eq. (37) appearing in Eqs. (36) and (37) P dimensionless quantity defined below Eq. (37) q volumetric flow rate; with an overbar it represents the mean flow rate R gas constant R, reaction rate of example II s dummy variable of integration t time; with a subscript zero it represents the initial time T absolute temperature T’ temperature fluctuation u heat transfer coefficient U dimensionless quantity below Eq. (37) U* rate of heat removal by cooling system V volume of reactor W Wiener process
N. J. RAO et al.
1204 x y
state variable; also dimensionless concentration fluctuation dimensionless temperature fluctuation
A = O(h)
means lim 4
Greek symbols a, p, y dimensionless
quantities defined below Eq. (39) 6 Dirac delta function A dimensionless flow rate fluctuation (r standard deviation 7 dimensionless time in Eqs. (36) and (37) defined below Eq. (37) 0 average holding time in continuous reactor .$ dimensionless concentration q dimensionless temperature
Subscripts c f n s
refers to coolant used in reactor cooling refers to feed to continuous reactor integer steady state
means lim$z=O h-0
111 _ _ Rao N. J., Borwanker J. D. and Ramkrishna A. SIAM J. Control (in press). r21 - - Chem. Funda_ _ Pell T. M. and Aris R.. Ind. Ennng mentals 1%9 8 339.  Aris R. and Amundson N. R., Chem. Engng Sci. 1958 9 250.  Wong E., Stochastic Processes in Information and Dynamical Systems, McGraw-Hill, N.Y. 1971.  Stratonovich R. L., Conditional Markov Processes and Their Application to the Theory of Optimal Control, American Elsevier, New York 1969. [clWang E. and Zakai M., Inc. J. Engng Sci. 1965 3 213. [71 Strat&ovich R. L., SIAM J. Control 1966 4 362. 181 Varadan S. R. S., Stochastic Processes. Courant Institute of Mathematical Sciences, New York University, N.Y. 1968. [91 Ito K., Lectures on Stochastic Processes. Tata Institute of Fundamental Research, Bombay, India 1960. I301 Aryaratham S. T. and Graefe P. W. U., ht. J. Control 1965 1 239; 1965 2 295. [Ill Rao N. J., Ph.D. Thesis, Indian Institute of Technology, Kanpur, India 1972.