Nomenclature
t 6%
= normalized dead time
a aT b
= =
bT
= = = = =
transportation lag or dead time normalized minor time constant minor process time constant, 0 < b < 1 process output normalized, transformed process output constants of integration process error, e = R - C
Pe
“‘dtz
= = = = = = =
= = = =
K, L
= =
aT) = = =
R* T
= = = =
= = =
initial crror after set point change error ai: first switch, e ( t 2 ) error dlerivative at first switch, e ( t 2 ) functions defined by Equation 3 forcing function process transfer function Hamiltonian function upper limit on M, a constant lower limit on M, a constant process steady-state gain process load variable manipulated variable normaliized transformed manipulated variable manipulated variable a delay time in the past desired value of C; controller set point initial set point predominant process time constant time, measured at controller time zero of set point change time for first switch tz a?” time for final switch, t 5 aT minimum time to reach origin from e(0) = el, e(0) = 0
+
-
x X
A&)
h
= time for error to reach and remain within
=t5% el of origin for set point changes (within i570L for step load changes) = state variable = state vector = conjugate system variable defined by Equation 4 = conjugate system vector
literature Cited (1) Balakirev, V. S., Automation Remote Control 24, 948 (1963).
(2) Biery, S. C., Boylan, D. R., IND.ENG.CHEM.FUNDAMENTALS 2, 44 (1963). (3) Bushaw, D. W.? “Differential Equations with a Discontinuous Forcing Term,” in “Contributions to the Theory of Nonlinear Oscillations,” Vol. 4, Princeton University Press, Princeton, N. J., 1958. (4) Cotter, J. E., Takahashi, Y., Chem. Eng. Progr. Symp. Ser. 46, 59, 119 (1963). (5) Coughanowr, D. R., Koppel, L. B., “Process Systems Analysis and Control,” McGraw-Hill, New York, 1965. (6 Desoer, C. A., Information Control 2, 333 (1959). (71 Eckman, D. P., “Automatic Process Control,” Wiley, New York, 1958. (8) Feldbaum, A. A., Automatika Telemekhanika 14, 712 (1953). (9) Grethlein, H. E., Lapidus, L., A.Z.Ch.E. J . 9, 230 (1963). 10) Hougen, J. O., Chem. Eng. Progr. Monograph Ser. 4, 60 (1964). 11) Kurzweil, F., ZEEE Trans. Auto. Control AC-8, 27 (1963). (12) Lapse, C. G . , I S A J . 3, 134 (1956). (13) .Latour, P. R., “Application of Pontryagin’s Maximum Principle to Time Optimum Control of a Second Order Overdamped System with Transportation Lag,” thesis for M.S. degree, Purdue University, W. Lafayette, Ind., June 1964. (14) Pontryagin, L. S., BoltyanskiI, V. G., Gamkrelidze, R. V., Mishchenko, E. F., “The Mathematical Theory of Optimal Processes,” Wiley, New York, 1962. (15) Williams, T. J., et al., “Approximation Models for the Dynamic Response of Large Distillation Columns,” Fourth Joint Automatic Control Conference, Minneapolis, 1963. RECEIVED for review July 20, 1964 ACCEPTED May 3, 1965
I
MONTE CARLO SIMULATION
OF REACTION KINETICS G E R A L D D. C O H E N l Operations Research Department, Allied Chemical Corp., New York, N . Y. Complex set!; of equations describing the concentration-time profiles of a chemical reaction can be solved by random number (Monte Carlo) simulation techniques. The basis for a valid solution is a stochastic model of a chemical reaction in which the reaction path is not fixed but allowed to vary from trial to trial. An example is presented to illustrate the convergence of Monte Carlo simulation to results predicted by the laws of mass action. A relationship between Monte Carlo sampling and numerical integration is demonstrated which should prove useful in the numerical solution of complex reactions (involving a series of equations) difficult to solve by analytic procedures.
differential equations describing concentration-time changes in a chemical reaction can often be formulated, but not solved by direct analytic means. Various methods have proved useful in this situation. T h e two most common are analog computer simulation and numerical integration on a digital computer. A third method, Monte Carlo simulation on a digital computer, can also be used, but has only recently received attention for this application ( 5 ) . I t has the important ability to include the realities of an industrial situation, such as changing temperatures and pressures, the addition or removal of catalysts and inhibitors, and general process charac-
T
1
HE
Present address, Mat:hematica, Princeton, N. J.
teristics. The simulation of reaction kinetics under these conditions can provide information that may reduce the amount of laboratory or plant experimentation. I t could also be used for experiments which cannot be physically performed, and a fruitful area of application of this idea is in the field of polymer chemistry. I t is possible to consider reactions containing several hundred species (representing different molecular weights and types, etc.), and arrive a t results either difficult to obtain or undeterminable by direct measurement, such as molecular weight distribution. The basic idea in a Monte Carlo simulation is to devise a set of sampling rules which, when followed sequentially, provide a realization of the reaction mechanisms. The VOL. 4
NO.
4 NOVEMBER 1 9 6 5
471
sampling is performed by selecting reaction units (corresponding to molecular concentration) one a t a time, at random, from a predetermined sample size, and deciding whether and how each “molecule” is to react. The reaction decision is dependent on the value of the reaction rate coefficient. High coefficients permit more frequent reactions and hence a faster rate of reaction. As an example of how the sampling works, consider the simplest case of a first-order reaction a t constant volume
There are two mutually exclusive ways in which X molecules of type A can remain at the end of the n l t h trial.
+
+ +
Either (X 1) existed before the trial and a successful reaction occurred, reducing the total to X, an event with probability (X 1 ) / T of selecting an A type molecule, times the probability of reaction, p , Or X already existed at the nth trial and no reaction occurred in the next trial. This has the probability of [(l - b) X / T ] .
P(X,n where A, = initial concentration of component A at time zero, moles per liter Bo = initial concentration of component B at time zero, moles per liter k = reaction rate coefficient, minute-’ at constant temperature A deterministic description of the concentration-time profile for this reaction is provided by the law of mass action.
+ 1) = p
P(X
+
~
T
+ 1,n) +
P(X,,O)= 1
(7)
n = 0
Equation 6 is a doubly discrete linear death process ( 3 ) . Its solution is facilitated with a Z transform defined as
Applying Equation 8 to 6, The solution of Equation 2 is
A(t)
=
(2
- l)Q(X,Z) - ZP(X,O) =
(3)
A,e-kt
The sampling procedure for this process requires that 1 “liter” of space be subdivided into T pockets. Into these pockets we put X o “molecules” of type A material, so that the fraction of filled pockets at time zero is proportional to the initial concentration, A,. (4) By induction : The process starts by first selecting a pocket at random. If it contains an A-type molecule, a decision is made as to whether this unit reacts. Successful reactions are recorded by reducing X , which represents the concentration of A, by 1 unit, and increasing Y by 1 unit (representing the concentration of B). The current trial number is increased by 1 after every selection, regardless of the outcome, and the process 1s repeated until completion. The decision whether a unit reacts is based upon the reaction probability computed for a small unit of time, At. Reaction probability p =
KAt -?
P
p
5
1
(5)
On a digital computer a pocket is examined for an A-type unit by generating a random number, r l , which is uniformly distributed between 0 and 1, and comparing it to the fraction X / T . If rl 5 X / T , the pocket contains an A unit. Similarly, for another random number, 72, if r2 5 p , the decision to react is successful.
and in general Q(X,
- r,Z)
=
(E) T
X,! X (X,- r ) !
~
Stochastic Description
If different sets of random numbers are used in the sampling just described, the course of the reaction will be different. We cannot speak of having exactly X molecules of type A at time t, but rather a probability of X molecules at time t. Let P(X,n) = probability that X molecules of type A remain after the nth sample is drawn p = probability that a given selected molecule of type A will react T = total number of pockets in system X, = initial number of type A molecules 472
I&EC FUNDAMENTALS
The inverted solution of Equation 13 is
[l
-1 (X, T
- i)]’
(14)
In order to relate real time to the sampling trial number, n, it is necessary from Equations 4 and 5 that
The incremental change in real time for each trial must be inversely proportional to the sample size. This is reasonable because with larger sample sizes more trials are required to make an equivalent amount of material react; hence the smaller effect of each trial. Using the time relation in Equation 1 4 yields in the limit
Tzm [1
Txt-8 I.1
Ai
(x,- i)]
--f
e-k(xo-i)t
(16)
accuracy may be acceptable. results. The conditions were
A , = 2 moles per liter T = 1000 pockets X, = 1000 molecules k = 0.10 minute?
The computing time was approximately 60 seconds on a digital computer (IBM 7074). Alternate Sampling Procedure 1
The amount by which each reactant changes during a trial of the basic sampling procedure is a variable, 6, with value 0 or 1. When a pocket is sampled and contains an A-type unit the conditional probability of the change variable is
The fit is good even for moderate size T , so that
P(6 P(6
E ( 6 / A )= 1 X
= 01A) = 1
-p
(22)
p
+ 0 X (1 -
p) = p
(23)
Therefore, instead of deciding which value the change variable takes, by drawing and comparing to a random number, we can credit the trial immediately with its expected change. This avoids the variability incurred by the use of random numbers. The entire procedure is
stands for
X,! (X,- X)!X! This is the same result arrived at by Bartholomay (7), who reasoned from the physical rather than the sampling process. The expected number of A-type molecules a t time t is by definition E(XIt') = 21 = XP(X1t)
Select a pocket at random (with probability X / T ) . If it contains an A-type unit, reduce the amount of these units by p, and increase B a like amount. The difference equation is
P(X,n -k 1 ) =*X
X
x,e
= 1/A)= p
The expected value of the change is
which is an expansion of binomial terms and reduces to
So that
Table I shows some typical
T
P(X
+ p,n) +
P(X,,O) = 1
2t
-
At
= A,ePkt
-k1
If we let X, = p M and X = p r , the solution can be written as (19)
This is seen to be identical with Equation 3. Thus, the interpretation of the deterministic equation based upon the law of mass action is that it describes the average behavior of a first-order process. The sampling procedure is therefore consistent with usually observed results. I t has been postulated that Equation 1 8 is a general description of first-order reactions (2) and can be used to explain variations from results predicted by deterministic analysis. Experimentally these variations may be measurable in very dilute solutions, and may also be observable in highly complex polymerization reactions when the formation of ultralong chains has extremely low probability. The standard deviation of process (Equation 1 8 ) around its average is
The average and standard deviations are
Table 1.
Trial No. 0
The variation of the sampling process (Equation 21) depends upon the absolute, not the relative, size of the sample, and for large sample sizes can be made negligible. This, however, requires a larger number of trials and computer time. The sampling procedure described for the example is correct, but there are various modifications which would reduce its variance. Two of these are described in the following section. However, evsen for samples of moderate size its
i , ooo
2,000 3; 000
Simulation of First-Order Reaction Time, Concentration of A Min. Simulated Au. Predicted 0 2.000 2.000 1 .o 1.806 1.810
2.0 3 .O
4,000 5,000
4.0 5.0
6,000 7,000
6.0 7.0
8 000 9,000 10,000
8.0 9.0
10.0
VOL.
1.620 1.474 1.356 1.198 1.078
1.000 0.896 0.792 0.710
1,637 1.482 1.341 1.213 1.098 0.993 0.899 0.813
0.736
4 NO. 4 N O V E M B E R 1 9 6 5
473
Since, p 5 1, the standard deviation of Equation 28 is less than Equation 21, and depending upon the scaling factor and T may be considerably less-Le., if p = 0.01, Equation 28 is only 1/10 the magnitude of Equation 21. Alternate Sampling Procedure 2
Before a pocket is selected, the unconditional probability of a change in the amounts of the reactants is
Higher Order Reactions
I n second-order reactions the basic sampling procedure is that described, but two pockets are selected a t random, and if they contain, for example, an A type and B type, a decision is made on whether they react. The average of this stochastic sampling process is not consistent with the deterministic process predicted by the law of mass action ( 3 ) . The difference, however, can be made insignificant. I n the reaction
X P(6 = 1) = - p T P(6 = 0 ) = 1 -
k
A+B+C
(29)
X - p T
(32)
the deterministic rate equation is (33)
and
E(6) =
X T
If the average of each side is computed,
- p
Proceeding as before, we can avoid the necessity of actually selecting a pocket and just credit each trial with its expected change, so that
dC-- k E [ ( A , - C ) ( B , -
kAOB0 and
P(X,,O) = 1
E(&)
=
T' ):.)
X , (I
=
-
1
$n
=
k [ ( A , - C)(B,-
=
+
I&EC FUNDAMENTALS
S~+Sj+&'~+j+RH
0
The result is a completely deterministic solution, with no variance, and is equivalent to numerical integration. I n this case X,,I = X,, - p X,,/T, which is convertible to = A 2 - h dAt/dt. These are the first two terms of the Taylor h ) , and this is the basis for the Euler series expansion of A ( t method of numerical integration. Other sampling plans can be devised which simulate the results of different methods of numerical integration. In situations where there are a number of reactants, all systems of numerical integration require that the incremental change for each reactant be computed, at least once, during each trial. The entire population is then credited with the change increments, and new increments are calculated from the new state of the system. The number of terms in each change calculation may be considerable, and involve not only all other reactants but previous changes as well. I n large scale reactions such as those present in high polymers, these calculations may be prohibitively time-consuming. The essence of the sampling procedures is to select a subset of the total terms necessary to progress from one time interval to the next in such a manner that although each repetition of the process gives different results the average remains unbiased. I n large scale systems (of second order) we have found that relatively few repetitions are necessary to give adequate results if judicious consideration is given to devising the sampling plan.
C)] + k V a r . C
which differs from the deterministic equation for the reaction by the factor Var. C. I n sampling experiments with the basic procedure using T = 10,000 the effect of the variance term was not observable. A stratified sampling plan can be easily devised for secondorder reactions. One of the two species can be predetermined and a sample drawn only to determine the other. The use of this method is illustrated in the reaction k
CX/t
474
(34)
C2is added and subtracted from Equation 34, noting that
dt
-
+ kE[C2 - (A, + BoIC1
ECz = Var. C - C2,
5 P(,(l
C)] =
dt
i = 1,
. . . .N
(36)
Instead of selecting two pockets, set i = 1 and sample one pocket at random, store the result, then set i = 2, sample again, store the result, and continue for all species. When i = N , update the state of the system with the stored results, increase the system time, and repeat for the next time increment. An efficient procedure can do this in 2 N multiplications and approximately 5,V additions. The simplest numerical integration would require about N ( N 1)/2 multiplications and N(N 1)/4 additions. Even with a number of repetitions of the sampling solution, so as to average the final result, the entire computing time is considerably less than numerical integration. On the basis of the multiplication time the sampling process could be repeated N ( N 1)/2/2'V A X/4 times. If, for example, N = 500, N / 4 = 125, but only 9 or 10 repetitions may be sufficient to yield a reasonable average. The net result is a reduction in computation time by a factor of 12. This factor can be crucial in simulating reactions which include incomplete mixing patterns or diffusion mechanisms because the number of possible reactions is compounded by a spatial identification.
+
+
+
Conclusions
It is not difficult to formulate sampling rules which, when executed, satisfy various mechanisms of chemical reactions. Realistic process conditions can be simulated and a large
number of reaction components considered. The major difficulty with the method, which has slowed its acceptance, is an incomplete understanding of the statistical techniques needed. I t is apparent from Equation 21, for example, that the sampling error is a function of the square root of the sample size. Large increases in samplc size have a diminishing return and other methods are needed if the variance is to be kept within acceptable bounds. Three variance reduction techniques, which increase the efficiency of the sampling process, have been described. These were based upon statistical analogies and in effect substituted a statistically equivalent sampling process for the direct method. A combination of methods is necessary to make large scale process simulations computationally feasible. There are interesting aspects of simulation as a research
tool ( 4 ) ,but its greatest immediate potential appears to be in the area of process optimization. Here the chemistry is largely known but the profusion of variables requires elaborate experimentation, and this is precisely where “computer experimentation” is useful. literature Cited (1) Bartholomay, A. F., Bull. Math. Biophys. 20, 175-89 (1958). (2) Zbid., 21, 363 (1959). (3) Bharucha-Reid, A. T., “Elements of the Theory of Markov Processes and Their Applications,” McGraw-Hill, New York,
1960. (4) Bunker, D. L., Sci. American (July 1964). (5) Schaad, J., A . I . Ch. E . J. 85, 3588 (1963). RECEIVED for review October 30, 1964 ACCEPTEDMay 7, 1965
COM MUN I CAT1ON
EFFECT OF HOLDUP ON BATCH DISTILLATION 0 PTI M IZATIO N The problem of determining the optimum distillate draw-off rate for batch distillation with holdup i s considered. The presence of holdup causes the optimum control to be on-off, and causes computational difficulty when exact rnethods are used. An approximate method was used and several cases were evaluated.
I N THIS work we are concerned with the following problem. I n batch distillation with holdup considered, what distillate rate policy should one use to obtain the maximum amount of overhead product of a specified purity in a specified time? The above problem statement can be written symbolically as follows.
lT lT lT D dt
Max: Subject to:
yD dt
=
nT
J
D[l
+ A@*
-J)]
(2)
dt
(3)
I t is important to realize that the above optimization is also constrained by the difkrential equations which represent the batch distillation. On-Off Control
Formally the prob1e:m is solved by choosing the control, D ( t ) , so that the Hamiltonian function, H ( t ) , is maximized a t every instant. Neglecting holdup, Converse and Gross (3) found
Hi
=
D {1- A
where y = y(X,,U). becomes
I.Y~‘ - YI)
+ ;[Xs-
Ly*
- KX,] - pl
(Xi- KX,)
+(X,- XI) f S P2
f
N
- (K
+ 1)‘5 + 4 , l
(7)
O n comparing these two expressions it is important to recognize that the inclusion of holdup has caused H to change from a nonlinear to a linear function of the control, D ( t ) . Hence when F # 0, D must take on one of its extreme values, 0 or V , depending whether F is < or > 0, respectively. Socalled on-off (bang-bang) control is optimal. If F = 0, D is not defined by the property that H is maximized-i.e.,the singular case. This situation is of interest only if it is of finite duration; otherwise it is merely a switching boundary. Hence we are interested in conditions for which
F=P=O
(8)
Now it can be shown (4) that H is constant; hence if P = 0, then G = 0. k and 6 are both functions of D but not identical functions (4). Hence, we conclude that D cannot satisfy both requirements simultaneously and therefore the singularLe., F = 0-solution does not exist. Computation
- Dp2
~ l p l
(4)
With holdup considered, the Hamiltonian El2 = FD $I G
S
-
P j + 2 [ W - l
0
D
V
E - p2
1
j=1
I t is convenient to combine the constraint with the objective function by employing a Lagrange multiplier, A, to obtain Max:
G
E
(1)
D dt
y*
where F
(5)
As is well known ( I ) , the differential equations that must be solved so that H can be evaluated, form a two-point boundary value problem. The initial conditions of the state variables-Le., concentrations on the individual trays and amount of material in the still-are known. Unfortunately, it is the final conditions of the adjoint variables that are known. VOL. 4
NO. 4
NOVEMBER 1 9 6 5
475
