E. STANLEY LEE
Invariant A Versatile Computational Concept This comparatively unknown mathematical approach eases the solution of many chemical engineering problems by reformulating the descriptive differential equations hen the complete initial conditions are given,
Wmodern digital computers are efficient tools for solving differential equations. But, unfortunately, many problems in chemical engineering are boundary value problems in which the conditions are not completely given at one point. Problems of this type are most subtle and difficult and are not suited for digital computers. Recently, various systematic approaches have been devised to reduce the boundary-value problem into a form more suitable for digital computers. These approaches can be classified into two general classes. The first class will be referred to as the usual or classical approach. A very promising technique in this class is the quasilinearization technique. The concepts used in the classical approach are basically the commonly used ones. I n general, the computational procedure is iterative. T h e second class of approach uses the invariant imbedding principle, which is completely different from the usual or classical approach to the boundary value problem. Instead of considering a single problem with a fixed duration or length of the independent variable, as used in the usual or classical approach, the invariant imbedding approach considers a family of problems with durations ranging from zero to the duration of the original problem. By imbedding this family of problems, expressions for the missing conditions of a two-point boundary-value problem can be obtained. Invariant Imbedding Concept
The invariant imbedding principle has its origin in the theory of semigroups. Ambarzumian (I) seems to have been the first one to employ this principle in any significant fashion. Chandrasekhar (6) generalized
the results of Ambarzumian and used this principle elegantly in treating the theory in radiative transfer. This principle has been studied extensively by Bellman, Kalaba, Wing, and coworkers (20). Invariant imbedding is only a concept; it is not a technique or method. This concept can be applied to a variety of different problems. Because of its completely different approach, it frequently gives some different insights to the same physical problem previously treated by the usual or classical method. Furthermore, these new formulations often have distinct advantages over the original formulations both computationally and theoretically. The functional equation of dynamic programming is a good example. Although dynamic programming technique has its own difficulties, it possesses certain distinct advantages over the original classical calculus of variations approach to optimization problems. T h e functional equation of dynamic programming is essentially the invariant imbedding equation with the addition of maximum or minimum operation. This paper serves to introduce the invariant imbedding concept and to discuss some of the more promising applications of the concept to chemical engineering problems. Emphasis will be placed on the use of the basic concepts of invariant imbedding as a computational tool. We shall not discuss the more elegant formulations and the theoretical or analytical applications of the invariant imbedding concept (72, 20). The approach will be introduced by considering a simple boundary-value problem. As a numerical example, the missing initial condition of a tubular reactor with axial mixing will be obtained recursively. We shall see that this recursive expression, resulting from the VOL. 6 0
NO. 9 S E P T E M B E R 1 9 6 8
55
application of the invariant imbedding concept, is tirneconsuming to solve numerically. In general, a different approach will be used. However, the numerical example is important in that it illustrates the invariant imbedding concept clearly. Furthermore, this direct solution of the recursive relation has certain applications. For example, it can be used to solve differential equations which are unstable when the usual numerical integration techniques are used. The recursive relations, which are the discrete expressions for the missing boundary conditions, can be reduced to a set of quasilinear partial differential equations of the initial value type. Further simplification of the partial differential equations to ordinary differential equations of the initial value type can be made for a large number of problems. Most of the applications discussed in this paper are based on this simplification. I n many engineering applications such as the parameter estimation problem and the numerical solution of semi-infinite boundary-value problems, the solution of a family or a sequence of two-point boundary-value problems instead of a single problem is required. We shall see that by using the invariant imbedding approach, the solution of a sequence of problems can be reduced to the solution of a single problem in which the sequences of missing boundary conditions are the unknown dependent variables. Obviously, it is impossible to discuss all of these applications in detail. Furthermore, some of these applications are still in the development stage. Our purpose is to show the type of problems to be treated by the concept and to emphasize the concept’s versatility.
different values from zero to t,, say a = 0, A, 2 A, . . , , then there will be a family of problems. Each member of this family has a different starting value of a and is represented by Equations 1 and 3. Let us consider obtaining the missing initial conditions y (a) for this family of problems. The idea is that neighboring processes are related to each other. The missing initial condition for the original problem y(0) may possibly be obtained by examining the relationship between the neighboring processes. Note that the missing initial condition y ( a ) for this family of processes is not only a function of the starting point of the process a, but also a function of the starting state or the given initial condition e. Define : r(c, u ) = the missing initial condition for the system represented by Equations 1 and 3 where the process begins at t = a withx(a) = c Then, obviously
y(4
=
r(c, a )
(4)
Note that ~ ( u andy(a) ) represent the starting state of the process. We shall consider r as the dependent variable, and c and a as the independent variables. An expression for r in terms of c and a will be obtained. Considering the A, the missneighboring process w-ith starting value a ing initial condition of this neighboring process can be related toy (a) by the use of Taylor’s series
+
y(a
+ A)
+ r’(a)A + O(A)
= y(a)
(5)
where O(A) represents higher order terms or terms involving powers of A higher than the first. At the starting value a, Equation 1 becomes
Illustration of Invariant Imbedding Approach = f [ x ( a ) , r(a), a1 = f [ c , d c , a ) , a1
(64
r’(a) = g [ x ( a ) ,y ( a ) , a1 = g[c, ~ ( c , a ) , a1
(6b)
.’(a)
T o illustrate the invariant imbedding approach, consider the nonlinear two-point boundary-value problem
Substituting Equations 6b and 4 into Equation 5, we have
+ g [ c , +, a), a l a 3. O(A> (7) O n the other hand, the following expression can be obtained for this missing initial condition y ( a + A) from y(a
with boundary conditions
+ A)
= r(c, a)
Equation 4 y(a
4- A)
= r[x(a
4- A>, a
+
with 0 5 t 5 t f . T o avoid the various computational difficulties in solving the above boundary-value problem, we shall convert it into an initial value problem. In other words, the missing initial condition y(0) will be obtained by using the invariant imbedding concept. To do this, consider the more general boundary conditions X(.)
= c
with a 2 t 5 t f ; a is the starting value of the independent variable 8 . However, it should be kept in mind that a also controls the duration of the process. If a assumes 56
INDUSTRIAL AND ENGINEERING CHEMISTRY
+ AI
(8)
Again, the expression x ( a A) can be related to its neighboring process .(a) = c by Taylor’s series x(a
+ A)
= ~ ( a -I)
x ’ ( n ) A -I- O(A) = c
+ fk,r(c,
a),
ala
+O(4
(9)
Thus, Equation 8 becomes y(a
+ A)
=
r{c +f[c,
r(6,
a), ala
+ O(A), + A ) a
(10)
Equating Equations 7 and 10, we obtain the desired relation r(c,a) +g[c,r(c,a),aIA = r { ~ + f [ ~ , r ( ~ , a ) , a l A , a f A }
omitting the terms involving powers of A higher than the first. The difference equation, Equation 11, can be used directly to obtain the missing initial condition r (c, a). Alternately, a partial differential equation can be obtained from Equation 11. Expanding the right-hand side of Equation 11 by Taylor's series, we have
r { c + f [ c , r(c, a>, ala, a
+A}
tion x , represents the concentration of A before entering the reactor. Equation 16 can be written as
dx - =Y dt
9 = Np,y + Np,Rx2 dt
=
(12) I n the limit as A tends to zero, the following first-order quasilinear partial differential equation is obtained from Equations 11 and 12 :
(19)
Instead of using the original boundary conditions, the following simpler boundary condition will be used x(0) = 1
Y(1) = 0
(20)
with 0 It I 1. I t will be seen later that if the range of x ( 0 ) used is large enough, the solution of Equation 11 for
If the process had zero duration, the missing initial condition would be equal to the given final conditions. Thus (14) t,) = 0 with c equal to all the possible values of x ( t ) . The missing initial conditions, r(c, a ) , for the family of processes with the starting values of the independent variable, a, from zero to t,, can be obtained by solving the systems in Equations 13 and 14. T(C,
Numerical Example
Since Equation 13 is a partial differential equation, various techniques can be used to solve this equation. However, to illustrate the invariant imbedding approach, the original difference equation, Equation 11, will be solved. The solution of this difference equation illustrates the concept clearly. Consider the simple chemical reaction (72)
A+A+B (15) which is taking place in a tubular reactor with axial mixing. The following equation can be obtained by material balance on reactant A : 1 d2x - dx -- Rx' 0 Npe dt2 dt
the systems of Equations 19 and 20 will also yield the original boundary condition as given by Equations 17 and 18. Instead of the above one boundary-value problem, consider the family of boundary-value problems with the following family of boundary conditions : x(a) = c
Y(1) = 0
(21)
with a t 5 1. Now, Equation 11 can be used for the systems in Equations 19 and 20. Identifying with the nomenclature of the previous section, we obtain f(X,
Y, t)
=
Y
g(x, Y, t ) = NPey
+ NP&'
and
t, = 1 Thus, Equation 11 becomes
r(c, a) f N d c , a>A-I- Np0Rc2A= r [ c
+ r(c, .)A, a + A ] (22)
Solving for r(c, a), one gets
with boundary conditions (79) (23) If A is small, the following approximation can be used
r [ c 4- r(c, a)A, a f A ] = r [ c f r(c, a where NPe is the Peclet group, R is the reaction rate group, x is the concentration of A, and t is the dimensionless reactor length which varies between 0 and 1. The concentra-
E. Stanley Lee is Associate Professor of Industrial Engineering at Kansas State University, Manhattan, Kan. This work was partly supported by the Air Force Ofice of Scient$ic Research, Ojice of Aerosjace Research, U.S.A.F., under contract F44620-684-0020, and by Federal Water Pollution Control Administration Grant WP-07 747-07. AUTHOR
+ A)A,
a
f A] (24)
and Equation 23 becomes
r(c,
a) =
1
+
1 NPeA
f r(c,
a
+ A M , + AI a
N ~ , R C ~ A(25) ] with the final condition
r(c, t,) = r(c, 1) = y(1) = 0, VOL 6 0
NO.
9
0
5
c
c
SEPTEMBER 1968
(26) 57
DURATION OF PROtfSS
0
.1 .2
0.3
INITIAL CONDITIONS
R = 2, Np. = 6, A = 0.1,
8
0.1
(28)
At the final point, a = t , = 1, the value of r is known for all values of c. I n the present case, this value is a constant value and is equal to zero. Next, let us obtain r(c, 1 A), or the values of r at a = 1 - A. Substituting a = 1 - A into Equation 25, we obtain
-
r(c, 1 Equation 25 can be solved in a badrward recursive fashion starting with the final condition, Equation 26, at
r(c, 1
Since we clearly cannot evaluate r at aU values of c, some discrete values of c can be chosen, say
8,
28,
...
(27)
and r i8 evaluated only a t these discrete values of c. Thus, there are grid points in both dimensions of the problem. The initial state of the system, c, is divided into c/8 grid points and the duration of the process, t,, is divided into t,/A grid points. As a concrete illustration, let us consider the following numerical values: 58
INDUSTRIAL A N D ENGINEERING CHEMISTRY
1 (rE + N 4
Since r(c, 1) = r [c
a+A=t,.
c = 0,
- A) = 1
- A)
+
r(C,
1 1 4 11 - N P ~ A } (29)
+ r(c; l)A, 1] = 0, thus =
1
1
+ NPJ I-NP&+AJ
(30)
.. . .
The values of r a t a = 1 A can be obtained from EquaThe results are listed in tion 30 for c = 0, 8, 28, Table I. With r at a = 1 A known, the values of r at a = 1 2 A can be obtained. Substituting a = 1 2 A into Equation 25, we have
-
-
-
r(c,l - 2 A ) =
+ NPJ (r[c + r(c, 1 - A)A, 1 - A] - Np@A) 1
1
(31)
Since the values of r at a = 1 - A are listed in Table I, the right-hand side of Equation 31 is completely known. The problem for c = 0.0 is trivial. For c = 0.1, from Table I we obtain ~(0.1,1
- 4)
=
-0.0075
with x, = 1. Now we wish to search the first and last columns of Table I11 for values of c and r(c, 0) so that the right-hand side of Equation 34a equals 1. For x ( 0 ) = c = 0.8, we find that x8 =
0.8
or r[0.1
+ ~(0.1,1 - A)A, 1 - A ] = ~(0.09925,1 - A)
.
T(C,
1
- 34)
1
{r[c
1 4-h T d
1
+ r(c, 1 - 2A)A,
- 2A] - Np,Rc2A)
(33)
Equation 33 can be solved in the same way as Equation 31--i.e., by using Table I1 with linear interpolation. This procedure can be continued until the values of r(c, 0) are obtained. Note that in calculating r(c,l -kA), only the values of r[c,l - (k - l ) 4 ] are needed. There is a table similar to Tables I and I1 for each r(c,l - kA), k = 1 , 2 , . . . , 10. These tables cover all the grid values of c and a. The last table for r(c, 0) for the above problem is shown in the first two columns of Table 111. The influences of the sizes of the intervals, A and 6, on the accuracy of the obtained missing initial conditions are shown in Table 111. The results in the last column are obtained with different values of A and 6. The missing initial condition for the original problem y (0) is given by the last row of Table 111, where c = 1.0 and a = 0. The calculated values of r(c, u ) with A = 6 = 0.01 are shown in Figure 1. Only a few seconds computation time are needed to calculate and print out all the tables for A = 6 = 0.1. However, over half a minute is needed to perform the same calculations when A = 6 = 0.01. An IBM 7094 computer was used in this work. The solution of Equation 25 not only gives the missing initial condition of the original system represented by Equations 19 and 20, but also gives the missing initial conditions for all the problems with boundary conditions represented by Equation 21 with 0 5 a It , and 0 5 c I c (Figure 1). Since all the initial conditions for all the interested values of u and c are obtained, the missing inital condition for the problem with the given boundary conditions of Equations 17 and 18 also can be obtained. Equations 17 and 18 can be written as xg = x ( 0 )
(0) - Y-
= 0.95761
remembering thaty(0) = r(c, 0). find that
(32) The value of r in Equation 32 can be obtained from Table I by linear interpolation. Once the value of r in Equation 32 is known, the value of r(0.1, 1 - 2A) can be obtained from Equation 31. This same procedure can be used to obtain r(c, 1 - 24) for c = 0.2, 0.3, . ., 1.O, and Table I1 is formed. Once the values of r a.t a = 1 - 2A are known and tabulated, the values of r(c, 1 - 3A) can be obtained. From Equation 25, we get
- -0.94565 6
X,
= 0.9
- -1.1687 6
For x ( 0 ) = 0.9, we
=
1.09478
By linear interpolation, we obtain the following corresponding values for x , = 1 : ~ ( 0= ) 0.83090, ~ ( 0 = ) -1.01458 The actual values obtained in a previous paper are (73): ~ ( 0= ) 0.83129, ~ ( 0 = ) -1.0122
The above computational scheme is exactly the same as the dynamic programming technique except that no search procedure is needed. T o obtain a fairly accurate value for the missing initial condition, the values of 4 and 6 must be fairly small. This computational scheme is obviously time-consuming. Furthermore, for problems with a large number of differential equations, the fastaccess memory requirement would severely limit the application of this scheme. This is exactly the same dimensionality problem encountered in the dynamic programming technique. Although Equation 13 can be solved easily by the method of characteristics, the numerical solution of the general case which consists of a system of first-order partial differential equations obviously is not simple. We shall see that a system of partial differential equations generally results from a boundary-value problem with a general system of ordinary differential equations. However, there are various ways to reduce the partial differential equations into a system of ordinary differential equations of the initial value type.
(344
NPe
YO)
=
0
(34b) VOL. 6 0 NO.
9
SEPTEMBER
1968
59
Systems of Equations Consider the following general system of nonlinear differential equations
functional equation of dynamic programming is essentially the invariant imbedding equation with the addition of maximum or minimum operations. To see this, let us consider the following problem in the calculus of variations: Find that function z ( t ) such that the function x ( t ) given by the differential equation dx
-= dt i
...,n
1,2,
=
x(0) = c
Xd(0)
=
Y&r)
=Y
XI0
(36) J(z)=
with 0 It 5 t,. T o obtain expressions for the n missing initial conditions, let us consider the problem with the more general boundary conditions : &(a)
=
Cd
Y&,)
=
Yi
(37)
1, 2, . . ., n
with a 5 t 5 t,. Following the approach used in the previous sections, we define r I ( c l , cg, . , . , c,; u ) = the missing initial conditions for the system represented by Equations 35 and 37 where the process begins at t = a with xi(.) = c$, 0bviously yi(a) =
.
rl(c1, c2,
.,;,c
,
a),
i
=
1, 2, . . .,n
(38)
The invariant imbedding equations representing these missing initial conditions are (72) :
gi[c, r(c, a ) ,
i
=
.I
(39)
1,2, , . . , a
where c represents the given inital vector with components c1, cg, . , . , c,; and r represents the n-dimensional missing initial vector. T o obtain the final conditions for Equation 39, observe again that if the process had zero duration, or a = t,, the missing initial conditions would be equal to the given final conditions. Thus T ~ ( C I , CZ,
. . .,
c;,
t,)
=
ytf,
Z
= 1, 2,
,
.
,,
TI
(40)
Dynamic Programming
T o show the connection between dynamic programming and invariant imbedding, the functional equation of dynamic programming will be formulated. Emphasis will be placed upon the use of the invariant imbedding concept. Some of the details of the derivation will be omitted. These details can be found in published monographs (2, 12, 76). The dynamic programming technique probably is better known than invariant imbedding. However, the 60
(42)
maximize the integral
d
. . . ,n
i = 1,2,
=
2)
and the initial condition
with boundary conditions
i
f(X,
INDUSTRIAL A N D E N G I N E E R I N G CHEMISTRY
1"
h(x, z)dt
(43)
The function x ( t ) represents the state of the system and the function z ( t ) is the control or decision variable. If this problem is treated by the classical calculus of variations, the result will be a two-point boundary-value problem. T o use the invariant imbedding concept, observe that when the maximum of Equation 43 has been obtained, the integral is a function only of the initial condition c and the duration of the process t,. Thus, we wish to imbed the original problem with particular values of c and t , within a family of processes in which c and t , are parameters. Define g(c, t,) = the maximum value of J where the starting state of the process is c and the total dura-
tion is t f Thus
The maximization is executed by choosing the proper value of z over the interval (0,t,). The function g will be referred to as the optimum return and J, which in general is not the optimum or maximum value, will be called the return or nominal return. The control variable z ( t ) also is known as a policy. The optimum value of z ( t ) is the optimal policy. Let us perturb the duration of the process and consider the original process with duration from t = 0 to t = t,, and a neighboring process with duration from t = A to t = t,. However, instead of relating these two processes as has been done in the previous sections, a different approach, using the property of g(c, t,) and the additive property of the integral, will be employed. The original process can be assumed to be composed of two different processes. The first process has a duration of t = 0 to t = A and the second or neighboring process has a duration oft = A to t = t,. We may write
In the limit as A + 0, Equation 53 becomes The second term represents the maximum return from the second process. Obviously, the starting state for this second process is
which is obtained from Equation 41. From the definition ofg and Equation 44, one can see that the maximum return from the second process is
Substituting Equation 47 into Equation 45, we have PA
The terms under the integrals may be approximated by PA
J
h(x, z)dt = h[c, z(O)]A
0
(49) with terms involving A2 and higher orders of A omitted. Equation 48 now becomes g(c, t,) = max ( h b , 4 0 ) IA ~ ( 0A),
+ g { + fb,4 0 )14 t f - A 1) c
(50) I t is interesting to compare the above equation with Equation 11. The present optimum return function g corresponds to the missing inital condition r of Equation 11. If we use a procedure completely similar to that used in obtaining Equation 13, the continuous version of the functional equation of dynamic programming can be obtained. Thus, using Taylor series, we have g {c
+
f[C,
4 0 ) 14 t , - A
1 = g(c, t,) +
Equation 50 becomes
Sinceg(c, t,) is independent of the choice of z, put it outside of the maximum operation sign. Thus,
0 = max (“[c, z(0, A)
z(O)]A
+ f [ c , z(O)]A Q d c
t,)
-
The above derivations are well known (2, 76). Our purpose is to interpret the derivations in terms of the invariant imbedding concept and to compare these derivations with the derivation of the invariant imbedding equation, Equation 13. To obtain the initial condition, observe that if the process had zero duration, the profit would be zero
Equation 54 is the desired relationship. This equation should be compared with the functional equation of invariant imbedding, Equation 13. Note that the initial condition, Equation 14, is also for a process with zero duration. Equation 54 is essentially a mathematical expression of the well known property of an optimal process, namely the principle of optimality. The advantages and disadvantages of dynamic programming are well known and are discussed in published literature (2, 76). Notice how we have changed a basically two-point boundary-value problem into a n initial value problem by the invariant imbedding approach. Estimation Problem
A basic model-building problem in chemical engineering is : Given the forms of the differential equations and the noisy measurements of the dependent variables of the process as functions of the independent variable, estimate the parameters or coefficients of the differential equations so that these equations can describe the process accurately. Note that we not only have the noisy measurements of the process, but also know that certain differential equations can represent the process. I n other words, we assume a certain knowledge concerning the process. This is much more practical than the usual black box approach. For a physical process, we almost invariably have a certain amount of knowledge concerning this process. Very little work has been done to solve this important nonlinear estimation problem until recently (3, 9-72, 74, 75, 77, 78). This is essentially an extension of the well known linear prediction problem treated by Kalman and Bucy ( 9 ) . If the differential equations are fairly simple and can be solved analytically, the above estimation problem can be solved fairly easily by the usual estimation procedures. I n engineering applications, very few equations can be solved analytically. Since the parameters to be estimated cannot be measured directly and only the dependent variables can be measured, this estimation problem is very difficult to solve. The quasilinearization technique has been used to solve this problem (72, 74). I n this section, we shall discuss how the invariant imbedding concept can be used to estimate the parameters or coefficients in the differential equations. Since the invariant imbedding approach is different from the usual VOL. 6 0 NO.
9
SEPTEMBER 1968
61
approach, we shall see that several advantages are obtained. Before we discuss the approach, let us consider a simple example to illustrate the problem. Consider the estimation of the rate constant from noisy measurements of the concentration of reactant A A+A+-B
(56)
If we assume second-order reaction, the reaction rate equation is
(57) where x represents the concentration of A, k is the reaction rate constant, and t represents time for a batch reactor or volume of reactor divided by the feed rate for a flow reactor. The dependent variable in this example is the concentration, x; the parameter or coefficient to be estimated is the rate constant, k . Inasmuch as Equation 57 can be solved analytically, there is no particular difficulty in estimating k from noisy measurements on x. However, many chemical reactions are much more complex. For example, the differential equations resulting from the analysis of catalytic reactions based on the Hougen-Watson-Langmuir-Hinshelwood models (8) are much more complex and generally cannot be solved analytically. Since the differential-reaction rate data frequently cannot be obtained, there is a need for a method to estimate the rate constants directly from concentration measurements in an integral reactor. More discussion on these aspects can be found in Reference ?2. The estimation of parameters or coefficients in differential equations has fairly general practical values. Not only must the dynamic model of a chemical process be represented by differential equations, but differential equations frequently are needed to represent even the steady-state process. One example would be the packed bed reactor, where not only the reaction rate, but also the axial diffusion rate represented by the Feclet group must be estimated (72). To illustrate the approach, let us consider a simple estimation problem. Assume that from our knowledge concerning the process we know that it can be represented by the nonlinear differential equation dx =A dt
X ,
t>
where f is a known function. From the measured data on x as a function oft, we wish to estimate the value of x. There are two kinds of noise or disturbance involved in the measured data. The first kind is the unknown disturbance on the input. Owing to the presence of this disturbance, the process must be represented by
(59) where u ( t ) is the unknown disturbance on the input and g ( x , t ) is a known function. The second kind of dis62
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
turbance is the measurement error. value of x ( t ) is
zl(t)
= x(t)
Thus, the measured
+ (measurement errors)
(60)
I n practical situations, we frequently find that x ( t ) cannot be measured directly, and only a certain function of x can be measured. Call this measurable function h(x, t ) ; we have z ( t ) = h(x, t )
+ (measurement errors)
(61) with h as a known function of x and t. By using the classical least-squares criterion, this estimation problem can be stated as: Based on the measuret t,, estimate the ment z ( t ) at various points oft, 0 unknown condition, x ( t f ) , for differential Equation 59 such that the following integral is minimized :
<
