- Email: [email protected]
Chapter 7 Quasilinearization and Invariant Imbedding
Chapter 7
QUASILINEARIZATION AND INVARIANT IMBEDDING
1. Introduction
The quasilinearization and invariant imbedding procedures can be combined in at least two different ways [ l , 21. A predictor-corrector formula can be formulated in which invariant imbedding predicts the missing boundary condition and the quasilinearization procedure corrects this predicted value. As can be seen from Table 6.3, the predicted missing initial condition given by the invariant imbedding equations is not very accurate unless a very small step size d is used. On the other hand, the problem may not converge to the correct solution when the quasilinearization procedure is applied if the guessed initial value is too far removed from the correct value. Thus, invariant imbedding can be used to predict a better initial approximation for the quasilinearization procedure, while quasilinearization can be used to correct the predicted initial condition to obtain a more accurate solution. A second way to combine these two procedures is to use invariant imbedding to avoid the numerical solution of algebraic equations in the quasiliniearization procedure. We have fairly accurate schemes for obtaining the solution of a set of ordinary differential equations. However, as we have seen in the previous chapters, serious difficulties exist in solving a set of linear algebraic equations. T h e invariant imbedding approach can be used to avoid these difficulties. Some computational algorithms by the combined use of dynamic programming and quasilinearization are also discussed. It is shown that by the combined use of these two techniques the dimensionality for a large number of problems can be reduced to one. T h e aim of this chapter is to present a variety of computational algorithms by the use of quasilinearization, invariant imbedding, and 217
218
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
dynamic programming. No effort is made to cover all the combinations. Only some of the more promising ones are discussed. Most of the numerical examples discussed in this chapter have been solved in previous chapters by quasilinearization. Thus, a comparison can be made between the present procedures and those discussed in earlier chapters. 2. The Predictor-Corrector Formula T o illustrate this procedure, consider the tubular flow reactor with axial mixing. This problem has been solved by using the difference equation of invariant imbedding in Section 3, Chapter 6 . Let us see how the predictor-corrector formula can be used to improve the results obtained there. The system is represented by the following equations (see Section 3, Chapter 6): x'=y
(la)
+ NpeRX'
( 1b)
x(u)
=c
(24
Y(1)
=0
(2b)
y' = Npey
By considering a family of processes with different starting points, a, the imbedding equation is obtained. In difference form the equation is
with final condition r(c, t t ) = r(c, 1) = 0
0
(4)
where r(c, a ) represents the missing initial conditions of the different processes with different starting points, a, and with different values of the starting state variable, c. The system represented by Eqs. (1) and (2) also can be solved by the quasilinearization technique. These equations can be linearized by the generalized Newton-Raphson formula 4+1
= Y%+l
2. THE PREDICTOR-CORRECTOR FORMULA
219
with boundary conditions %+l(4 = c Yn+1(1) = 0
(6b)
However, the quasilinearization approach only can solve these different processes with different values of a and c separately. Thus, if the following grid points for a and c are considered
,..., 0
a=1,1-Ad,1-2A c = 0 , 6 , 26 )...,c
+
(7) (8)
+
there are ( l / d 1) (c/S 1) different processes to be solved independently by the quasilinearization technique. However, as shown in the results in Section 3, Chapter 6, the processes with a = 1 for all values of c and the processes with c = 0 for all values of a are trivial. Consequently, there are ( l i d ) (c/S) different processes only. With the numerical values listed in Eqs. ( 3 . 2 ~ and ) (3.10) of Chapter 6 and A t = 0.01 (9) the above equations are solved by the predictor-corrector approach. The Runge-Kutta numerical integration scheme has been used and A t is the integration step or integration interval. T h e missing initial conditions predicted by invariant imbedding for processes with starting point a = 1 - d are listed in Table 6.1. These missing initial conditions can be corrected by the quasilinearization technique. From Table 6.1, the missing initial condition for the process with c = 0.1 is obtained: T(C,
1 - 0) = y(1
-
A)
=
-0.0075
(104
The given initial condition for the process is x(l - 0 ) = c
=
0.1
(lob)
< <
with 1 - A t 1. T h e problem represented by Eqs. (1) and (10) is an initial-value problem that can be integrated easily. Using the results from this integration as the initial approximations for x and y , more accurate solution for the problem represented by Eqs. (l), (2b), and (lob) can be obtained by using the recurrence relations, Eq. (5), with the boundary conditions xn+Jl - 0 ) = c = 0.1 (114 Yn+1(1) = 0
(11b)
220
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
Equations ( 5 ) and (1 1) constitute a linear boundary-value problem. If the state vector ~ % + ~is( tdefined ) as the two-dimensional vector with components ~ ~ + ~and ( t~ )% + ~ (the t ) ,general solution for these equations is represented by xn+,(t) = XP(n+l)(t) Xh(n+l)(t)an+l (12) ~ ) (the t ) homogeneous solution The particular solution vector ~ ~ ( ~ + and matrix Xh(n+l)(t)can be obtained by integrating (5) and its homogeneous form, respectively, with the following initial values:
+
Xh(n+1)(1- 4
(14)
=1
with c = 0.1. Since one initial condition is missing, only one set of homogeneous solutions is needed. The iterative procedure discussed in Section 16 of Chapter 3 can be used to obtain the more accurate solutions. More accurate solution for the next process with c = 0.2 can be obtained by the same procedure using T(C,
1 - 0) = ~ ( 1 A)
=
-0.030
(154
as the missing initial condition in obtaining the initial approximations. Obviously, the given initial condition is x(1
A)
= c = 0.2
(154 with 1 - A t 1. This procedure can be continued for all the processes with a = 1 - A and c = 0.3, 0.4,..., 1.0. The correct missing initial conditions from these calculations are shown in Table 7.1. -
< <
TABLE 7.1 MISSING INITIALCONDITIONS v(c, a) FOR a = 1 - A AFTER CORRECTION BY QUASILINEARIZATION C
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .o
Y(C,
1 - A)
-0.0089706 -0.035674 -0.079805 -0.14107 -0.21919 -0.31388 -0.42489 -0.55195 -0.69483 -0.85328
2. THE
PREDICTOR-CORRECTOR FORMULA
22 1
Only one iteration of the quasilinearization technique is needed to obtain these missing initial conditions starting with the predicted values listed in Table 6.1. Further iterations do not change the values shown in Table 7.1. With r at a = 1 - d known, the values of r at a = 1 - 24 can be predicted approximately by Eq. (3.13) of Chapter 6 and Table 7.1. The computational procedure is the same as that described in Chapter 6. These predicted values are listed in the second column of Table 7.2. TABLE 7.2 MISSING INITIAL CONDITIONS r(c, a) FOR a = 1 - 24
C
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .o
Predicted r(c, 1 - 24)
Corrected r(c, 1 - 24)
-0.01 3056 -0.05 1701 -0.1 1518 -0.20277 -0.31379 -0.44760 -0.60358 -0.781 14 -0.97972 -1.1988
-0.013735 -0.054027 -0.11962 -0.20938 -0.32228 -0.45741 -0.61392 -0.79104 -0.98806 - 1.2043
More accurate values can be obtained from these predicted values by the same procedure as that used in obtaining Table 7.1 with the exceptions that a = 1 - 24 and that the duration of these processes equals 24. Thus, from the second column of Table 7.2, the missing initial conditions for c = 0.1 is obtained T(C,
1 - 24) = ~ ( 1 24) = -0.013056
(164
and the given initial condition for this process is x(1 - 24) = c = 0.1
< <
(16b)
with 1 - 24 t 1. Initial approximations can be obtained from Eqs. (1) and (16). Then the quasilinearization procedure can be applied to Eq. (5) with the boundary conditions xn+l(l - 24)
= c = 0.1
rn+l(l) = 0
222
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
This procedure, again, can be continued for all the processes with a = 1 - 24 and c = 0.2, 0.3,..., 1.0. These more accurate values are shown in the third column of Table 7.2. These values are obtained from the values listed in the second column of Table 7.2 with only one iteration. Further iterations do not change these values. T h e correct values for r at a = 1 - 3 4 also can be obtained by the same procedure. This procedure can be continued for all the grid points of a. T h e last table for r(c, 0) is shown in Table 7.3. Similar tables exist for a = 1 - 34, 1 - 4d,... . TABLE 7.3 MISSING INITIALCONDITIONS ~ ( c a) , FOR a
C
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .o
=
0
Predicted 7(c, 0)
Corrected
-0.019028
-0.018822 -0.071737 -0.15482 -0.26529 -0.40104 -0.56040 -0.74199 -0.94464 -1.1674 - 1.4093
-0.072429 -0.15617 -0.26743 -0.40406 -0.56436 -0.74695 -0.95067 -1.1766 - 1.4225
7(c,
0)
I n the above calculations, not only the missing initial conditions for all the processes, but also the complete solutions for Eqs. (1) and (2) for all the values of a and c with 4 = 6 = 0.1 have been obtained. A total of 100 processes have been solved. Approximately over 10 seconds are needed to perform these computations on the IBM 7094 computer. T h e solution for the process with a = 0 and c = 1 is listed in Table 7.4. T h e improved values listed in Table 7.3 should be compared with the values predicted by invariant imbedding, listed in the second column of Table 6.3. T h e initial conditions for the original problem with boundary conditions (3.16) of Chapter 6, again, can be obtained from the third column of Table 7.3 by linear interpolation. With x, = 1, these conditions are ~ ( 0=) 0.83104 y ( 0 ) = -1.01378 T h e actual values from Table 3.1 are ~ ( 0=) 0.83129
y ( 0 ) = -1.0122
3.
223
DISCUSSION
TABLE 7.4 SOLUTION FOR
THE
PROCESS WITH a
=
0 AND
c =
1
tk
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.o
1.o
- 1.4093
0.87490 0.77537 0.69462 0.62807 0.57258 0.52606 0.48733 0.45623 0.43412 0.42508
- 1.1095
-0.89227 -0.73026 -0.60595 -0.5073 1 -0.42500 -0.34988 -0.27006 -0.16582 0.30517 x
Values almost as good as those shown above have been obtained in Chapter 6 with A = 6 = 0.01. However, only the missing initial conditions, not the complete solutions have been obtained there in spite of the fact that more computing time has been used. 3, Discussion
T h e quasilinearization procedure does not have to be applied to every process. Suppose only the process with a = 0 and c = 1 is to be considered; then the missing initial condition can be predicted by invariant imbedding as has been discussed in Section 3 of Chapter 6. From the second column of Table 6.3 at c = 1.0, we obtain the predicted missing initial condition for this particular process: r(1,O)
= y(0) =
-1.4274
Using this value, initial approximations for x and y can be obtained and thus the quasilinearization procedure can be used to obtain a more accurate y(0). At the same time, the complete solution for this process with a = 0 and c = 1 also is obtained. As has been discussed in Section 14 of Chapter 6,-the use of the difference equations of invariant imbedding to obtain the missing initial conditions is limited to problems with small number of variables only. T h e computational procedure to be discussed in the next section has a much wider applicability. This procedure is based on linear differential equations which can be obtained from nonlinear equations by the generalized Newton-Raphson formula.
224
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
4. Linear Boundary-Value Problems
T h e invariant imbedding approach can be used to replace the numerical solution of algebraic equations in the quasilinearization procedure. To illustrate this, consider the nonlinear second-order differential equation d2x + q1(t) -ax& dt2
+42(W
=
O
with boundary conditions x(0) = 1
with 0
< t < tf . Equation (1) can be written as dx
=Y
dY
dt = -4dt)Y
-
Applying the generalized Newton-Raphson formula, we obtain x:+1
= Y,+l
Yh+l = -41WYfi+1
- 2&)XnXn+1
+ 42(t)X:
with boundary conditions
%+do) = 1 Y n + d f f )= 0
Equations (4) and ( 5 ) represent a linear two-point boundary-value problem. Two different approaches have been used to solve this problem in Chapter 3. T h e first one is based on the superposition principle; and the second one is the finite-difference method. Both of these approaches have limitations. Because the rapid-access memory in current computers is limited, the finite-difference method is limited to problems with small number of variables. T h e problem of ill-conditioning in solving the algebraic equations can cause considerable difficulties when the first approach is used. Invariant imbedding can be used to form the third approach. This approach avoids the numerical solution of algebraic
4. LINEAR
225
BOUNDARY-VALUE PROBLEMS
equations. However, both this and the first approaches may have stability problems caused by the marching integration technique. T o use the invariant imbedding approach, consider again the problem with the more general boundary conditions (64
%+,(a) = c
YTZ,l(tf) = 0 (6b) An expression for the missing initial condition for the system represented by (4)and (6) can be obtained by considering a family of processes with different starting points, a. Let Yn+1(4 = rTZ+,(c, 4
(7)
+
where Y,+~ represents the missing initial condition for the (n 1)st quasilinearization iteration. We can obtain the invariant imbedding equation from Eq. (12.5) of Chapter 6
=
-q1(4rn+1(c,
4 - 2qz(a)xn(a)c + 4 2 ( a ) g ( 4
(8)
remembering that the variables with subscript n are known variables. T h e final condition for Eq. (8) is rn+,(c, 5) = 0
(9)
Since Eq. (4)is linear, the function r,+, can be expressed as rn+,(c, 4 = c 5 n + 1 ( 4
+
rln+1(4
(10)
Substituting (10) into (8) and equating terms involving like powers of c, we obtain
Equation (11) could have been obtained directly from Eq. (13.5) of Chapter 6. From Eqs. (9) and (lo), we see that 5n+1(tr) = 0
.ITZ+dtr) = 0
(12)
Equations (11) and (12) constitute a nonlinear initial-value problem. T h e missing initial condition for the original systems (4) and (6) can be obtained by integrating these equations.
226
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
Equations (1) and (2), which are now represented by the recurrence relations (4),( 5 ) , (1l), and (12), can be solved in the following iterative manner. With an assumed initial approximation for xnSo(t)= xnz0(u), Eq. (1 1) can be integrated backward starting at the final condition (12). Thus, we have the functions {,+,,,(a) and ~ ~ + , , ~ (for a ) all the processes with different starting points, a. We can then obtain the missing initial condition for the original process with a = 0 from the equation
remembering that the value of c for the original process is x,+,(O) = 1. Now, Eqs. (4),(5a), and (13) with n 1 = 1 become an initial-value problem and the results of the first quasilinearization iteration can be obtained by integrating these equations with the assumed initial approximation xo(t). With xl(t) and yl(t) known, Eq. (1 1) can be integrated with n = 1 to obtain the functions {n+l=2(a) and ~ , + ~ , ~ ( aand ) hence the missing initial condition, ~ , + , = ~ ( 0Thus, ) . the results for the second iteration can be obtained. This procedure can be continued until no more improvements on the values of x and y can be obtained. Notice that there are two independent variables in the above equations. T h e independent variable in systems (4)and (5) is t, while in the invariant imbedding equations (1 1) and (12) it is a. Notice also that the missing initial condition is obtained by integrating (1 1) backward starting with (12), which obeys the given final condition. Then, on the basis of this newly obtained missing initial condition, the solutions of (4)and ( 5 ) are obtained by integrating (4)forward. Thus, the calculated value of the given final condition, y,+,( tr), depends on two integration processes which are performed in series. Consequently, the present procedure is very sensitive to the accuracy of the integration processes. Obviously, this procedure is not suited to problems which are unstable when the marching integration technique is used. However, since no system of algebraic equations is solved, it seems to be more suited to problems with ill-conditioned algebraic systems.
+
5. Numerical Results
T h e problem formulated above is solved with the following numerical values: t , = 1.0 4 d t ) = -6 (1) q z ( t ) = -12 At da = 0.01
5.
227
NUMERICAL RESULTS
where d t and da are the integration step sizes for Eqs. (4.4) and (4.1 l), respectively. With the initial approximation equal to the given initial condition, X,Jtk)
=
1.0
Fz
=
0, 1, 2 ,..., N
(2)
TABLE 7.5
CONVERGENCE RATESWITH INVARIANT IMBEDDING
Iteration 1
0 0.2 0.4 0.6 0.8 1 .O
1 .O 0.78880 0.66685 0.59664 0.55740 0.54222
-1.3723 -0.79251 -0.45732 -0.261 74 -0.13737 -0.1io10 x 10-4
-2.7445 -2.7442 -2.7409 -2.7083 -2.3933 0.0
1.3723 1.3721 1.3705 1.3541 1.1967 0.0
Iteration 2
1 .O
1 .O 0.77547 0.62838 0.52591 0.44961 0.38042
- 1.4091 -0.89141 -0.60511 -0.43401 -0.34280 -0.38498
0 0.2 0.4 0.6 0.8 1.o
1 .O 0.77537 0.62807 0.52608 0.45633 0.42546
- 1.4093
0 0.2 0.4 0.6 0.8 1 .O
1 .O 0.77537 0.62807 0.52608 0.45634 0.42559
- 1.4093
0 0.2 0.4 0.6 0.8
1 .O 0.77537 0.62807 0.52608 0.45634 0.42560
- 1.4093
0 0.2 0.4 0.6 0.8
-2.5155 -2.1346 - 1.9007 - 1.7281 - 1.4025 0.0
1.1064 0.76436 0.59321 0.49566 0.38497 0.0
Iteration 3 -0.89226 -0.60591 -0.42483 -0.26931 -0.17480
-2.5003 -2.0697 -1.7571 - 1.4901 - 1.0994 X
lo-'
0.0
1.0909 0.71251 0.49766 0.35904 0.23208 0.0
Iteration 4 -0.89226
-0.60591 -0.42481 -0.26921 0.38191 x lo-'
-2.5002 -2.0695 -1.7583 - 1.5015 - 1.1418 0.0
1.0908 0.71235 0.49838 0.36489 0.25083 0.0
Iteration 5
1 .O
-0.89226 -0.60590 -0.42481 -0.26919 0.38776 x lo-'
-2.5002 -2.0695 - 1.7583 -1.5016 -1.1419 0.0
1.0908 0.71235 0.49839 0.36490 0.25087 0.0
228
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
with t,,, - t, = A t , the results listed in Table 7.5 are obtained. The Runge-Kutta integration scheme is used for both Eqs. (4.4) and (4.11). Since y,(t) does not appear in Eq. (4.4), the initial approximation yn=O(tk) is not needed. Each iteration takes approximately less than half a second to compute on the IBM 7094 computer. Table 7.5 shows that the functions 5 and 7 at a = 0.8 have the slowest convergence rates. Furthermore, it shows that the accuracy of y ( t ) at t near tf is not very good. T h e value of y(tf) should be zero, but a value of 0.38776 x is obtained. Thus, it seems that the accuracy of the solution for x and y is not as good as the results obtained by the superposition principle. With the numerical values listed in Eq. ( l ) , this problem is the same as that solved in Section 3 of Chapter 6 or Section 2 of Chapter 7 for c = 1 and a = 0. Thus, Table 7.4 can be compared with Table 7.5. T h e accuracy can be improved by a reduction in the step size, A t , used. T h e influence of the initial approximation on the convergence rates is shown in Table 7.6. Only the missing initial and missing final conditions are shown. It can be seen that the convergence rates are almost independent of the initial approximations within the range tested. No convergence has been obtained with xo(t,) = 5, k = 0, 1, 2, ..., N . Except for the differences in the given initial conditions, the present problem is the same as that solved in Section 5 of Chapter 3. Thus, approximate comparisons also can be made with the results obtained there. That convergence has been obtained with xo(t,) = 10 in Chapter 3 seems to suggest that the present approach has a smaller interval of convergence. I n order to avoid interpolations, the two step sizes, A t and d a , should be equal. Equation (4.11) must be integrated backward. It can be TABLE 7.6 AS
CONVERGENCE RATESOF x(1) AND y ( 0 ) FUNCTIONS OF INITIAL APPROXIMATION, INVARIANT IMBEDDING
W0(tk) =
Iteration
x(1)
0.0 y(0)
XO(tk) =
0.01
XO(tk)
=
0.1
XO(tb) =
2.0
41)
Y(0)
41)
Y(0)
41)
Y(0)
0.96751 0.46717 0.38156 0.42560 0.42559
-0.03945 -1.3748 -1.4091 - 1.4093 - 1.4093
0.74046 0.14186 0.41168 0.42561 0.42559
-0.35713 -1.3919 -1.4093 - 1.4093 - 1.493
1.0 0.54222 0.38042 0.42546 0.42559
0.0 -1.3723 -1.4091 - 1.4093 -1.4093
~ _ _ _ _
~
1 2 3 4 5
1.0 0.54222 0.38042 0.42546 0.42559
0.0 -1.3723 -1.4091 -1.4093 -1.4093
6.
OPTIMUM TEMPERATURE PROFILES IN TUBULAR REACTORS
229
shown easily that the same Runge-Kutta formula can be used for backward integration,. provided that the integration step size d a is replaced by --da.
6. Optimum Temperature Profiles in Tubular Reactors
In order to test the combined quasilinearization and invariant imbedding procedure and also to compare it with other approaches, the problem solved in Section 3 of Chapter 5 will be solved by this approach. T h e equations for the state variables and Lagrange multipliers are (see Section 2 of Chapter 5 )
dhl = k,h,
dt
-
k,X,
with boundary conditions X,(O)
= x;
T h e control variable T ( t ) is obtained from the following algebraic equations: XlklEl(h1
-
A,)
+ x,k,E,X,
=0
(3)
and
Equations (1)-(4) represent the system to be solved. T h e above equations have been linearized in Section 2 of Chapter 5, by the generalized
230
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
Newton-Raphson formula. These linearized equations have been listed in (2.17) of Chapter 5 and can be represented symbolically as follows:
d xl.n+l 4 1 i . n ( t ) ~ i . n + i+ qi2.n(t)Xzmn+i + ~ i i . n ( t ) L n + i dt
dxz,n+l
-dt
+ -
qZl,n(t)Xl.n+l
+
+
+ hz,n(t)
+
vlz.n(t)&,n+l
4 .n+l - ~ z l . r r ( t ) ~ l . n + l -dt
+
+
(54
+
q Z ~ . n ( t ) ~ ~ . n + l~ l , n ( t ) h l , n + l
Pzz,n(t)hz,n+1
4 .n+l - Wll.n(t)Xl.n+l -dt
+
+h l M
P1z,n(t)hz.n+1
“zz.n(t)~z,n+l
wlZ.n(t)X*,nil
+
(5b)
+
Vll.n(t)l,.n+l
(5 4
~l,n(t)
Wzz.n(t)Xz,n+i
+ %(t)
+
Wz1,n(t)Ln+1
(54
with boundary conditions xl,n+l(o)= x;
(64
4
(6b)
XZ.,+l(O)
=
hl,n+l(G)
=0
(64
1 (64 The control variable T(t) has been eliminated from the differential equations by using Eqs. (3) and (4) before linearization. The coefficients q, p , w , v , h, and u are functions of the variables x ~ , x~ ~, , Al,n, ~ , A,,n, which are known functions and are obtained from the previous nth iteration. The exact forms of these coefficients can be obtained by comparing Eq. (5) with Eq. (2.17) of Chapter 5. The system represented by Eqs. (5)and(6)has been solved in Chapter 5 by the superposition principle. This system also can be solved by obtaining the missing initial conditions by invariant imbedding. Thus, the original boundary-value problem becomes an initial-value probiem that can be solved by marching integration techniques. T o use the invariant imbedding approach, the more general boundary conditions hz,n+l(tf)
=
Xl.ra+l(4
= c1
(74
Xz,n+1(4
= CZ
(7b)
6.
OPTIMUM TEMPERATURE PROFILES I N TUBULAR REACTORS
23 1
< <
with a t tf again can be considered. Equations ( 5 ) and (7) are essentially the same equations as those listed in (13.1) and (13.2) of Chapter 6 with M = 2. T h e h in Eq. (5) corresponds to y in Eq. (13.1) of Chapter 6. If 4*n+l(a) = T1.n+&1 > c2 > 4 (84 Az*n+1(4 = Tz.n+1(C1
3
c2
9
4
(8b)
the invariant imbedding equations are
Equations (9) and (10) are obtained from Eqs. (13.5) and (13.6) of Chapter 6. T h e functions 5 and 77 are defined by Tl*n+l(Cl
9
CZ >
TZ*n+l(Cl
t
CZ
4 = C1511.n+1(4 4 = C1521.n+1(4
+ +
C251z.n+1(4 C2522.n+1(4
+ +
%,@+1(4
(114
rlz,n+1(4
(Ilb)
With the nth iteration known, the missing initial conditions for the (n 1)st iteration, Al,n+l(0)and A2,n+l(0), can be obtained by integrating (9) backward starting at the final conditions, Eq. (10); and by using Eqs. (8) and (11) at a = 0. Once these two missing initial conditions are obtained, the results for the (n 1)st iteration can be obtained by integrating ( 5 ) numerically. This iterative process can be continued.
+
+
232
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
7. Numerical Results
Using the equations and procedures discussed in the preceding section, some numerical experiments have been performed. This procedure will not converge if the numerical values listed in Eqs. (3.1)-(3.3) of Chapter 5 are used with (1)
da = d t
The values of the functions <(a) and ~ ( a increase ) rapidly during the numerical integration of the first iteration and extremely large values for these functions have been obtained at a = 0. However, convergence can be obtained if the initial approximations are obtained by using the following assumed temperature profile: T,(t,) = 335°K
k
= 0,
1, 2 ,...,N
(2)
These initial approximations have been obtained by first integrating Eqs. (6.la) and (6.lb) with (6.2a) and (6.2b) as the initial conditions and with T(tk)given by (2). T h e values of the Lagrange multipliers are then obtained by integrating Eqs. ( 6 . 1 ~ and ) (6.ld) backward with (6.2~)and (6.2d) as the final conditions. With the initial approximations obtained in this manner and with the numerical values listed in Eqs. (3.1) and (3.2) of Chapter 5 , the results listed in Table 7.7 are obtained. A figure very similar to Fig. 5.2 would be obtained if the convergence rate of the temperature profile were plotted. T h e optimum values for the functions and 77 are shown in Figure 7.1. Table 7.7a should be compared with
<
a
FIG.7.1. Optimum values of 1 and 7 , tt
=
8.
7.
233
NUMERICAL RESULTS
Table 5.2. Notice that the convergence rates for these two approaches are almost the same. However, since no convergence was obtained when the given boundary conditions were used as the initial approximations, the present procedure seems to have a smaller interval of convergence than the procedure used in Section 3 of Chapter 5. TABLE 7.7a CONVERGENCE RATESWITH t, Iteration
0 1 2 3 4
XI@,)
0.21 120 0.17760 0.17058 0.17044 0.17043
X&f)
0.67248 0.67737 0.67935 0.67943 0.67943
=
8
AND
To&) = 335°K
4(O)
UO)
T(O)
T(tf)
0.55788 0.60293 0.60946 0.61000 0.61000
0.87629 0.84683 0.82873 0.82823 0.82823
335.00 342.55 340.90 340.82 340.82
335.00 337.06 335.94 335.89 335.89
TABLE 7.7b
CONVERGENCE RATES
1 2 3 4 5
0.11490 0.12861 0.13607 0.13590 0.13591
WITH
-0.15062 -0.16847 -0.16547 -0.16514 -0.16515
tf
=
8
AND
-0.15062 -0.16847 -0.16547 -0.16514 -0.16515
TO&)= 335°K
0.35652 0.22436 0.20126 0.20072 0.20072
0.60680 0.61374 0.60903 0.60898 0.60898
0.77336 0.82154 0.82939 0.82945 0.82945
Some numerical experiments also have been performed with the more severe reaction conditions listed in (3.5) of Chapter 5. No convergence was obtained when the initial approximations were sbtained from the assumed temperature profile listed in (3.7) of Chapter 5. I n fact, even with the better initial temperature profile obtained from Eqs. (3.8) and (3.9) of Chapter 5, convergence still cannot be obtained. To obtain convergence, a smaller integration step size is needed, and thus a better accuracy in integration is obtained. Instead of Eq. (3.5) of Chapter 5, the following numerical values are used: t,
=
10 minutes
x! = 0.95 mole/liter
x i = 0.05 mole/liter
A t = 0.01
(3)
234
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
T h e initially assumed temperature profile is obtained by using Eqs. (3.8) and (3.9) of Chapter 5 ; and the initial approximations are obtained from this initially assumed temperature profile by integration. Use of the reaction rate constant listed in (3.1) of Chapter 5 gives the results shown in Table 7.8. T h e optimum values of the functions 5 and 7 are shown in Figure 7.2. Table 7.8a should be compared with Table 5.5. T h e conver-
0
2
4
6
10
8
a
FIG.7.2. Optimum values of 1 and 7,tt
=
10.
gence rates for these two procedures, again, are quite similar. However, Table 7.8 was obtained with a much smaller integration step size than that used in obtaining Table 5.5. T h e obtained value of A,($), which is Since this value should be zero, not shown in Table 7.8, is 0.23 x an accuracy of 0.23 x only is obtained at this final point. A much TABLE 7.8a
CONVERGENCE RATESWITH t,
0 1 2 3
4 5
0.16542 0.17709 0.17253 0.17265 0.17264 0.17264
0.67805 0.68069 0.68025 0.68012 0.68012 0.68012
0.67630 0.67704 0.67846 0.67859 0.67858 0.67859
=
10
0.71152 0.70862 0.70938 0.70948 0.70948 0.70948
345.00 360.44 359.96 359.94 359.94 359.94
335.00 335.72 336.07 336.11 336.1 1 336.11
8.
235
DISCUSSION
TABLE 7.8b
CONVERGENCE RATESWITH
1 2 3 4 5
0.85653 0.13502 0.11586 0.11549 0.11551
t, = 10
x -0.82157 x -0.82159 X 0.45777 x lo-' -0.20195 x lo-' -0.20195 X lo-' 0.36437 x lo-' -0.20494 x lo-' -0.20494 x lo-' 0.36339 x lo-' -0.20458 x lo-' -0.20458 X 10-l 0.36325 x lo-' -0.20460 x lo-' -0.20460 X lo-' 0.36326
0.66894 0.67819 0.67851 0.67851 0.67851
0.68652 0.71035 0.71078 0.71076 0.71076
higher accuracy has been obtained by using the superposition principle even with a much larger integration step size (see Table 5.4). No convergence can be obtained if the following values are used to replace the values listed in Eq. (3.8) of Chapter 5: To(0)= 360°K
T o ( t f )= 335°K
(4)
The other numerical values used remain the same. Notice that Figs. 7.1 and 7.2 are very similar. In fact, if the right-hand portion of Fig 7.2 for a = 2 to a = 10 is compared with Fig. 7.1, we find that they are almost the same. This is not surprising. T h e initial conditions for the problem with tr = 8 are xIo = 0.53 and x20 = 0.43; while from the problem with tf = 10, we obtain x1 ( t = 2) = 0.55255 and x2 ( t = 2) = 0.41405 (see Table 5.4). Thus, the problem with tf = 8 is almost the same as the problem with tf = 10 except that the former problem starts 2 minutes later. T h e initial temperature profile listed in Eqs. (3.8) and (3.9) of Chapter 5 has been estimated from the optimum temperature profile of the problem with tr = 8. T h e Runge-Kutta integration scheme has been used for all the numerical integrations discussed in this section.
8. Discussion
T h e combined quasilinearization. and invariant imbedding procedure is very sensitive to the numerical errors of the integration procedure. Thus, a much smaller integration step size is needed when this procedure is used. When the difficulties in solving the problem lie in the numerical integration process, this procedure is not suitable. On the other hand, it should be used when the difficulties in solving the problem are connected with the numerical solution of the algebraic equations.
236
7. QUASILINEARIZATION
A N D INVARIANT IMBEDDING
9. Dynamic Programming and Quasilinearization-I
For optimization problems, dynamic programming also can be used to form various combinations of numerical procedures [l, 3, 61. Quasilinearization and dynamic programming can be combined in various ways. In the first place, dynamic programming can be used to replace invariant imbedding in the predictor-corrector formula discussed in Section 2. A search procedure must be added when dynamic programming is used. Second, the original problem can be divided into a series of subprobbms. Quasilinearization can be used to obtain the optima of these subproblems first. Then dynamic programming can be used to obtain an approximate solution of the original problem by combining these subproblems in an optimum manner. To illustrate the procedure, consider the minimization of the integral
J
=
J t ’ f ( x , x’) dt 0
with the terminal conditions x(0) = cg
(24
4 t f ) = Cf
(2b)
Instead of solving the above problem directly, let us consider the problems with the shorter duration
(3) with t, = 0 and tN = tf . In other words, we have divided the original problem into (N - 1) subproblems. Now, let us see what the terminal conditions of these subproblems are. The terminal condition at t, is given by (2a). However, the terminal condition at t, is unknown. A series of terminal conditions will be allowed at this point. Suppose the value of x at t, can be reasonably assumed to be within the range xa
< <
xb
(4)
Then we can assume m terminal conditions at t, with the ith terminal condition as q t * ) = x,
+- i
x b - xa
m-1
i
= 0,1,2,..., ( m - 1)
9.
DYNAMIC PROGRAMMING AND QUASILINEARIZATION-I
237
Thus, in the interval [tl = 0, t,], m optimization problems must be solved by the quasilinearization technique. Each problem is represented by Eq. (3) with k = 1. T h e initial conditions for these problems are given by (2a) and the final conditions are given by ( 5 ) with i = o , 1 , 2 ,..., ( m - 1 ) . T h e next subproblem with the interval [t, , t,] can be treated in the same manner except that now both the initial and final terminal condi) be tions are unknown and all the possible values for x ( t z ) and ~ ( t , must considered. This procedure can be continued until the last subinterval in which the initial terminal is unknown. But, the final terminal is known and is given by Eq. (2b). Once the optima for all the subproblems are obtained, the optimum of the original problem with the interval [0, tr] can be obtained easily by dynamic programming. Consider the first subinterval, From co , we can go to any of the m points at t, . However, only one of these points gives the minimum value of (1). From these points at t , , we can go to any of the points at t, . This process can be continued. T h e problem is to find the best way to go from co to ct . This is a special case of the general routing problem which is especially suited for dynamic programming and has been considered in detail in various publications [4, 51. T h e functional equation of dynamic programming can be established easily. Let hidtk)
=
I
the minimum of
Jk with terminal conditions x t i l ( t k )and ~ ‘ j ) ( t ~ + ~ )
and gi(‘?k)
=
+
1
+ +
the minimum of the sum h,(t,) hil(tk+l) ... h i j ( t ~ - l ) ) where the process begins at tk with the starting terminal diI(tk)j
Note that the subscripts i and j in the definition of g are optimum values at times tk+l , t,,, ,..., t N - , . At time t, , the value o f j is optimum, but the value of i is not. On the other hand, at time t N P l ,the subscript i is optimum, but j must correspond to the given terminal condition ct . T h e principle of optimality yields the functional equation gdtk) = mp[hij(tk)
fg j ( t k + l ) ]
=
2,***,( N - 2,
(6)
T h e optimum of the last subinterval can be obtained by considering a process with duration t = t,-l to t = t , , gi(tN-1)
= hidtN-1) = hicf(tiv-l)
(7)
238
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
Thus, the optimum of the original process can be obtained by minimizing Eq. (6) recursively starting at the last subinterval represented by Eq. (7). This is essentially an optimization problem with N - 1 stages. The optimum obtained is, in general, not the true optimum, but only an approximate one. This approximate result can be used as the initial approximation and the quasilinearization procedure, again, can be used to obtain better results. 10. Discussion
As has been shown in Section 12 of Chapter 2, there are two ways to ensure convergence of the quasilinearization procedure. The first is to obtain an initial approximation which is sufficiently close to the desired solution; and the second is to make the interval of interest of the independent variable [ t o ,tr] sufficiently small. The simplest way to make this interval small is to consider a process with shorter duration first, and then use the results of this process to solve a process with longer duration. As has been discussed in Section 7, this approach has been used to obtain the initial temperature represented by Eqs. (3.8) and (3.9) of Chapter 5 from the results of the process of shorter duration with tf = 8. The procedure discussed in the preceding section constitutes a more sophisticated way to reduce the interval of interest of the variable. It is true that this procedure is fairly time consuming when x is multidimensional. However, if only an approximate initial approximation is needed, fairly small values of N and m can be chosen and hence the time required can be reduced considerably. 11. Linear Differential Equations
We have seen in previous chapters that the quasilinearization procedure is a useful tool for solving nonlinear boundary value problems. In the next section it will be shown that this procedure also can be used to reduce the dimensionality of an optimization problem when dynamic programming is used. However, before we discuss this approach, some well-known results concerning the solution of linear differential equations will be summarized. Consider the linear differential equations with variable coefficients
12.
DYNAMIC PROGRAMMING AND QUASILINEARIZATION-I1
239
I n vector-matrix notation, Eq. (1) becomes
with initial conditions x(0) = xo
(3)
where Q ( t ) is an M x M matrix, and p(t) and x are M-dimensional vectors. T h e solution of (2) is [7, 81
+ 1 X(t)X-'(s)p(s) ds t
x(t) = X(t)xo
0
(4)
where X ( t ) is an M x M matrix and is the solution of the matrix equation dX dt
- = Q(t)X
with initial conditions X(0)
where I is a unit matrix. At t
=
=
I
tr , Eq. (4) becomes
where c
=
X(tf)XO
K(s)
=
X(tt)X-'(s)
(7)
and c is an M-dimensional constant vector and K(s) is an M x M matrix. Equation (6) can be rewritten as
12. Dynamic Programming and Quasilinearization-I1
Now we wish to show that by the combined use of dynamic programming and quasilinearization a reduction in the dimensionality of an optimization problem can be obtained. Consider the problem of maximizing the function +I(b),
XZ(tf)!..., x m ( 5 ) )
(1)
240
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
over the control variables z(t), which are related to the state variables x by means of the nonlinear differential equations dxi -=ffi(x,z) dt
i = 1 , 2,..., M
(2)
1, 2, ..., M
(3)
with initial conditions X i ( 0 ) = x:
i
=
<
where m < M , 0 \< t tf , and x and z are M-dimensional vectors. I n addition, the problem must satisfy the constraints zi.min
< z i ( t ) < zi,max
i
= 1, 2,.-, M
(4)
For simplicity, the number of control variables has been assumed to be equal to the number of state variables. Obviously, the procedure to be discussed can also be applied if this is not the case. If the functional equations of dynamic programming are used, this problem involves the computation and storage of functions of M variables. However, Bellman and Kalaba [9, 101 have shown that if Eq. (2) were linear, the above problem could be treated by dynamic programming involving sequences of functions of m variables. Since m is equal to one or two for a number of significant problems, this is a very important reduction in the dimensionality of the problem. This is especially true in view of the fact that Eq. (2) can be linearized by the quasilinearization procedure. Thus, an iterative scheme can be formed by the combined use of dynamic programming and quasilinearization. Applying the generalized Newton-Raphson formula, Eq. (2) can be linearized
I n vector-matrix notation, Eq. (6) becomes d%+l dt - f (xn 9 zn)
+ J(xn
zn)(xn+i - xn)
+ Jz(xn
>
zn)(zn+1
- zn) (7)
where Jz(x,, z), is the same Jacobi matrix as that in Eq. (17.10) of Chapter 2, except that the differentiation is with respect to z not x.
12.
24 1
DYNAMIC PROGRAMMING AND QUASILINEARIZATION-I1
I n obtaining Eq. (7), we have considered both x and z as unknown functions. Comparing Eqs. (11.2) and (7), we obtain Q(t>
=
J(xn
(8)
2,)
+
2,) - J ( x n G)X, JAxn > ~ n ) ( ~ n + l2,) (9) Using Eqs. (8) and (9), the solution of (7) can be represented by Eq. (1 1.8). Introducing the Lagrange multiplier A, the problem becomes the maximization of
~ ( t= ) f (xn
J
7
= + ( ~ l ( t f )x*~ t f ) , * - *Xvm ( t r > ) -
J”
tf 0
f ( z )dt
(10)
over all z ( t ) ,satisfying Eqs. (2)-(4). Using Eq. (11.8), Eq. (10) becomes
where z,+,(t) are the unknown control variables after n iterations. If we consider X as a known parameter, the maximum value of Eq. (11) depends only on c1 , c2 ,..., c, , and tj . Furthermore, by examining Eqs. (11.5), (11.7), and (8), it can be seen that c is independent of the control variables zn+l(t). Thus, we wish to imbed the original problem with particular values of cl, c2 ,..., c, and duration tf within a family of processes in which cl, c2 ,..., c, and t j are parameters. Notice that if the explicit solution, Eq. (1 1.8), were not used, the maximum value of Eq. (10) would depend on c l , c2 ,..., cM and tf . We have reduced the number of variables from M to m. If m is equal to one or two, a feasible computational procedure has been obtained. Following the approach used in Chapter 6, we define g(c,
i
the maximum value of J where the process begins
> ~2 >.*.>
cm
2
a) = at t
=
u with starting state c1
,..., c,
.
1
Since the process is nonstationary, we have fixed the final time tf . A family of processes with different starting points a will be considered. T h e new maximization problem is
242
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
T h e maximization is excuted by choosing the proper values of z over the interval [a, tf]. Applying the principle of optimality, we obtain the desired recurrence relation
T h e terms under the integral sign may be approximated by
M
cm
+ c hmj(a)Pj(a) j=l
a
+A)]
(16)
T o obtain the final condition for Eq. (16), observe that if the process had zero duration or a = tf , then the maximum value of Eq. (11) would be equal to zero. Thus g(c, > cz ~
. * cm * ~
7
tj) = 0
(17)
Notice that we have assumed that the duration of the process tf is divided into small intervals of A width. Let tf = AN, then a = 0, d, 24, ..., N d . Thus, Eq. (16) can be solved in a backward recursive fashion starting with the known final condition, Eq. (17), at a = tf . T h e computational procedure can now be summarized as follows: (a) Estimate a reasonable control policy ~ ~ = ~satisfying ( t ) , Eqs. (4) and ( 5 ) .
13.
FURTHER REDUCTION I N DIMENSIONALITY
243
(b) Calculate xn=,(t)from Eqs. (2) and (3), using the newly obtained values of z,=,(t). (c) Obtain z,,,(t) by maximizing Eq. (1l), using the newly obtained values of z,=,(t) and ~ , = ~ ( tT) h. e values of z,,,(t) must satisfy Eq. (4)(d) Return to (b) with n = 1. Equation (11) can be maximized by using the recurrence relation, Eq. (16), remembering that the solution of the linearized equation, Eq. (7), has been used in obtaining this recurrence relation.
13. Further Reduction in Dimensionality
Owing to the limited rapid-access memory of current computers, the above algorithms cannot be used if m is larger than three. However, for a large number of problems, a further reduction of the dimensionality can be obtained. Consider the problem of maximizing the function
Let us introduce a new state variable, xM+,(t), defined by
Differentiating Eq. (2) with respect to t , we have
T h e initial condition is X M + m = H(x(0))
(4)
If we consider Eqs. (12.2) and (3) as the system of differential equations, the objective function, Eq. (12. l), becomes a function of one variable
4 = XM+l(tf)
(5)
If the algorithms obtained in the previous section are used, a problem with a dimensionality of one is obtained. This is a significant reduction in terms of computational requirements.
244
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
14. Discussion
T h e above procedure can be generalized easily to problems with more general objective function. For example, the problem of maximizing the integral
J
=
/”f(x, 0 z) dt
(1)
can be treated by the above procedure if we introduce a new state variable, ~ ~ + ~defined ( t ) , by
with initial condition %4+1(0) = 0
(3)
T h e problem now becomes the maximization of ~ ~ + ~A( variety t ~ ) .of other forms of objective functions also can be treated by the algorithms obtained in the previous section. More discussion can be found in the references listed at the end of the chapter [ll-131. Let us consider briefly how Eq. (12.16) may be solved. I n order to obtain the values of c and K(a), the homogeneous differential equation, Eq. (1 1.5), must be solved first. Since Q ( t ) does not contain the unknowns x , + ~and z,+~ , Eq. (1 1.5) can be solved easily. Thus, in actual computations c and K(a) in Eq. (12.16) can be considered known. T h e functions f and p contain the unknown control variables Z ~ + ~ ( U T ) . h e problem is to find Z,+~(U) so that the expression inside the square bracket on the right-hand side of Eq. (12.16) is maximized. Since dynamic programming is especially suited for discrete processes, the present approach in the reduction of dimensionality appears to be a promising tool for solving stagewise processes. REFERENCES 1. Bellman, R., Kagiwada, H. H., and Kalaba, R., Numerical studies of a two-point nonlinear boundary-value problem using dynamic programming, invariant imbedding, and quasilinearization. RM-4069-PR. RAND Corp., Santa Monica, California, March, 1964. 2. Bellman, R., and Kalaba, R., Dynamic programming, invariant imbedding and quasilinearization: Comparisons and interconnections. RM-4038-PR. RAND Corp., Santa Monica, California, March, 1964; see also Bellman and Kalaba in “Computing Methods in Optimization Problems.” Balakrishnan, A. V., and Neustadt, L. W., eds., Academic Press, New York, 1964.
REFERENCES
245
3. Bellman, R., and Kalaba, R., “Quasilinearization and Nonlinear Boundary-Value Problems.” American Elsevier, New York, 1965. 4. Bellman, R., and Dreyfus, S., “Applied Dynamic Programming.” Princeton Univ. Press, Princeton, New Jersey, 1962. 5. Kalaba, R., On some communication network problems. “Combinatorial Analysis.” Am. Math. SOC.,Providence, Rhode Island, 1960. 6. Kalaba, R., Graph theory and automatic control, in “Applied Combinatorial Mathematics” (E. F. Beckenbach, ed.). Wiley, New York, 1964. 7. Bellman, R., “Introduction to Matrix Analysis.” McGraw-Hill, New York, 1960. 8. Bellman, R., “Stability Theory of Differential Equations.” McGraw-Hill, New York,
1953. 9. Bellman, R., Some new techniques in the dynamic programming solution of variational problems, Quart. Appl. Math. 16, 295 (1958). 10. Bellman, R., and Kalaba, R., Reduction of dimensionality, dynamic programming, and control processes, J . Basic Eng. 83, 82 (1961). 11. Rozonoer, L. I., The maximum principle of L. S. Pontryagin in optimal system theory, Automatika Telemekhhanika 20, 1320, 1441, 1561 (1959), [English transl. Automation Remote Control 20, 1288, 1405, 1517 (1960)l. 12. Katz, S., Best operating point for staged systems, Ind. Eng. Chem. Fundamentals 1, 226 (1962). 13. Fan, L. T., “The Continuous Maximum Principle.” Wiley, New York, 1966.
QUASILINEARIZATION AND INVARIANT IMBEDDING
1. Introduction
The quasilinearization and invariant imbedding procedures can be combined in at least two different ways [ l , 21. A predictor-corrector formula can be formulated in which invariant imbedding predicts the missing boundary condition and the quasilinearization procedure corrects this predicted value. As can be seen from Table 6.3, the predicted missing initial condition given by the invariant imbedding equations is not very accurate unless a very small step size d is used. On the other hand, the problem may not converge to the correct solution when the quasilinearization procedure is applied if the guessed initial value is too far removed from the correct value. Thus, invariant imbedding can be used to predict a better initial approximation for the quasilinearization procedure, while quasilinearization can be used to correct the predicted initial condition to obtain a more accurate solution. A second way to combine these two procedures is to use invariant imbedding to avoid the numerical solution of algebraic equations in the quasiliniearization procedure. We have fairly accurate schemes for obtaining the solution of a set of ordinary differential equations. However, as we have seen in the previous chapters, serious difficulties exist in solving a set of linear algebraic equations. T h e invariant imbedding approach can be used to avoid these difficulties. Some computational algorithms by the combined use of dynamic programming and quasilinearization are also discussed. It is shown that by the combined use of these two techniques the dimensionality for a large number of problems can be reduced to one. T h e aim of this chapter is to present a variety of computational algorithms by the use of quasilinearization, invariant imbedding, and 217
218
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
dynamic programming. No effort is made to cover all the combinations. Only some of the more promising ones are discussed. Most of the numerical examples discussed in this chapter have been solved in previous chapters by quasilinearization. Thus, a comparison can be made between the present procedures and those discussed in earlier chapters. 2. The Predictor-Corrector Formula T o illustrate this procedure, consider the tubular flow reactor with axial mixing. This problem has been solved by using the difference equation of invariant imbedding in Section 3, Chapter 6 . Let us see how the predictor-corrector formula can be used to improve the results obtained there. The system is represented by the following equations (see Section 3, Chapter 6): x'=y
(la)
+ NpeRX'
( 1b)
x(u)
=c
(24
Y(1)
=0
(2b)
y' = Npey
By considering a family of processes with different starting points, a, the imbedding equation is obtained. In difference form the equation is
with final condition r(c, t t ) = r(c, 1) = 0
0
(4)
where r(c, a ) represents the missing initial conditions of the different processes with different starting points, a, and with different values of the starting state variable, c. The system represented by Eqs. (1) and (2) also can be solved by the quasilinearization technique. These equations can be linearized by the generalized Newton-Raphson formula 4+1
= Y%+l
2. THE PREDICTOR-CORRECTOR FORMULA
219
with boundary conditions %+l(4 = c Yn+1(1) = 0
(6b)
However, the quasilinearization approach only can solve these different processes with different values of a and c separately. Thus, if the following grid points for a and c are considered
,..., 0
a=1,1-Ad,1-2A c = 0 , 6 , 26 )...,c
+
(7) (8)
+
there are ( l / d 1) (c/S 1) different processes to be solved independently by the quasilinearization technique. However, as shown in the results in Section 3, Chapter 6, the processes with a = 1 for all values of c and the processes with c = 0 for all values of a are trivial. Consequently, there are ( l i d ) (c/S) different processes only. With the numerical values listed in Eqs. ( 3 . 2 ~ and ) (3.10) of Chapter 6 and A t = 0.01 (9) the above equations are solved by the predictor-corrector approach. The Runge-Kutta numerical integration scheme has been used and A t is the integration step or integration interval. T h e missing initial conditions predicted by invariant imbedding for processes with starting point a = 1 - d are listed in Table 6.1. These missing initial conditions can be corrected by the quasilinearization technique. From Table 6.1, the missing initial condition for the process with c = 0.1 is obtained: T(C,
1 - 0) = y(1
-
A)
=
-0.0075
(104
The given initial condition for the process is x(l - 0 ) = c
=
0.1
(lob)
< <
with 1 - A t 1. T h e problem represented by Eqs. (1) and (10) is an initial-value problem that can be integrated easily. Using the results from this integration as the initial approximations for x and y , more accurate solution for the problem represented by Eqs. (l), (2b), and (lob) can be obtained by using the recurrence relations, Eq. (5), with the boundary conditions xn+Jl - 0 ) = c = 0.1 (114 Yn+1(1) = 0
(11b)
220
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
Equations ( 5 ) and (1 1) constitute a linear boundary-value problem. If the state vector ~ % + ~is( tdefined ) as the two-dimensional vector with components ~ ~ + ~and ( t~ )% + ~ (the t ) ,general solution for these equations is represented by xn+,(t) = XP(n+l)(t) Xh(n+l)(t)an+l (12) ~ ) (the t ) homogeneous solution The particular solution vector ~ ~ ( ~ + and matrix Xh(n+l)(t)can be obtained by integrating (5) and its homogeneous form, respectively, with the following initial values:
+
Xh(n+1)(1- 4
(14)
=1
with c = 0.1. Since one initial condition is missing, only one set of homogeneous solutions is needed. The iterative procedure discussed in Section 16 of Chapter 3 can be used to obtain the more accurate solutions. More accurate solution for the next process with c = 0.2 can be obtained by the same procedure using T(C,
1 - 0) = ~ ( 1 A)
=
-0.030
(154
as the missing initial condition in obtaining the initial approximations. Obviously, the given initial condition is x(1
A)
= c = 0.2
(154 with 1 - A t 1. This procedure can be continued for all the processes with a = 1 - A and c = 0.3, 0.4,..., 1.0. The correct missing initial conditions from these calculations are shown in Table 7.1. -
< <
TABLE 7.1 MISSING INITIALCONDITIONS v(c, a) FOR a = 1 - A AFTER CORRECTION BY QUASILINEARIZATION C
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .o
Y(C,
1 - A)
-0.0089706 -0.035674 -0.079805 -0.14107 -0.21919 -0.31388 -0.42489 -0.55195 -0.69483 -0.85328
2. THE
PREDICTOR-CORRECTOR FORMULA
22 1
Only one iteration of the quasilinearization technique is needed to obtain these missing initial conditions starting with the predicted values listed in Table 6.1. Further iterations do not change the values shown in Table 7.1. With r at a = 1 - d known, the values of r at a = 1 - 24 can be predicted approximately by Eq. (3.13) of Chapter 6 and Table 7.1. The computational procedure is the same as that described in Chapter 6. These predicted values are listed in the second column of Table 7.2. TABLE 7.2 MISSING INITIAL CONDITIONS r(c, a) FOR a = 1 - 24
C
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .o
Predicted r(c, 1 - 24)
Corrected r(c, 1 - 24)
-0.01 3056 -0.05 1701 -0.1 1518 -0.20277 -0.31379 -0.44760 -0.60358 -0.781 14 -0.97972 -1.1988
-0.013735 -0.054027 -0.11962 -0.20938 -0.32228 -0.45741 -0.61392 -0.79104 -0.98806 - 1.2043
More accurate values can be obtained from these predicted values by the same procedure as that used in obtaining Table 7.1 with the exceptions that a = 1 - 24 and that the duration of these processes equals 24. Thus, from the second column of Table 7.2, the missing initial conditions for c = 0.1 is obtained T(C,
1 - 24) = ~ ( 1 24) = -0.013056
(164
and the given initial condition for this process is x(1 - 24) = c = 0.1
< <
(16b)
with 1 - 24 t 1. Initial approximations can be obtained from Eqs. (1) and (16). Then the quasilinearization procedure can be applied to Eq. (5) with the boundary conditions xn+l(l - 24)
= c = 0.1
rn+l(l) = 0
222
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
This procedure, again, can be continued for all the processes with a = 1 - 24 and c = 0.2, 0.3,..., 1.0. These more accurate values are shown in the third column of Table 7.2. These values are obtained from the values listed in the second column of Table 7.2 with only one iteration. Further iterations do not change these values. T h e correct values for r at a = 1 - 3 4 also can be obtained by the same procedure. This procedure can be continued for all the grid points of a. T h e last table for r(c, 0) is shown in Table 7.3. Similar tables exist for a = 1 - 34, 1 - 4d,... . TABLE 7.3 MISSING INITIALCONDITIONS ~ ( c a) , FOR a
C
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .o
=
0
Predicted 7(c, 0)
Corrected
-0.019028
-0.018822 -0.071737 -0.15482 -0.26529 -0.40104 -0.56040 -0.74199 -0.94464 -1.1674 - 1.4093
-0.072429 -0.15617 -0.26743 -0.40406 -0.56436 -0.74695 -0.95067 -1.1766 - 1.4225
7(c,
0)
I n the above calculations, not only the missing initial conditions for all the processes, but also the complete solutions for Eqs. (1) and (2) for all the values of a and c with 4 = 6 = 0.1 have been obtained. A total of 100 processes have been solved. Approximately over 10 seconds are needed to perform these computations on the IBM 7094 computer. T h e solution for the process with a = 0 and c = 1 is listed in Table 7.4. T h e improved values listed in Table 7.3 should be compared with the values predicted by invariant imbedding, listed in the second column of Table 6.3. T h e initial conditions for the original problem with boundary conditions (3.16) of Chapter 6, again, can be obtained from the third column of Table 7.3 by linear interpolation. With x, = 1, these conditions are ~ ( 0=) 0.83104 y ( 0 ) = -1.01378 T h e actual values from Table 3.1 are ~ ( 0=) 0.83129
y ( 0 ) = -1.0122
3.
223
DISCUSSION
TABLE 7.4 SOLUTION FOR
THE
PROCESS WITH a
=
0 AND
c =
1
tk
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.o
1.o
- 1.4093
0.87490 0.77537 0.69462 0.62807 0.57258 0.52606 0.48733 0.45623 0.43412 0.42508
- 1.1095
-0.89227 -0.73026 -0.60595 -0.5073 1 -0.42500 -0.34988 -0.27006 -0.16582 0.30517 x
Values almost as good as those shown above have been obtained in Chapter 6 with A = 6 = 0.01. However, only the missing initial conditions, not the complete solutions have been obtained there in spite of the fact that more computing time has been used. 3, Discussion
T h e quasilinearization procedure does not have to be applied to every process. Suppose only the process with a = 0 and c = 1 is to be considered; then the missing initial condition can be predicted by invariant imbedding as has been discussed in Section 3 of Chapter 6. From the second column of Table 6.3 at c = 1.0, we obtain the predicted missing initial condition for this particular process: r(1,O)
= y(0) =
-1.4274
Using this value, initial approximations for x and y can be obtained and thus the quasilinearization procedure can be used to obtain a more accurate y(0). At the same time, the complete solution for this process with a = 0 and c = 1 also is obtained. As has been discussed in Section 14 of Chapter 6,-the use of the difference equations of invariant imbedding to obtain the missing initial conditions is limited to problems with small number of variables only. T h e computational procedure to be discussed in the next section has a much wider applicability. This procedure is based on linear differential equations which can be obtained from nonlinear equations by the generalized Newton-Raphson formula.
224
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
4. Linear Boundary-Value Problems
T h e invariant imbedding approach can be used to replace the numerical solution of algebraic equations in the quasilinearization procedure. To illustrate this, consider the nonlinear second-order differential equation d2x + q1(t) -ax& dt2
+42(W
=
O
with boundary conditions x(0) = 1
with 0
< t < tf . Equation (1) can be written as dx
=Y
dY
dt = -4dt)Y
-
Applying the generalized Newton-Raphson formula, we obtain x:+1
= Y,+l
Yh+l = -41WYfi+1
- 2&)XnXn+1
+ 42(t)X:
with boundary conditions
%+do) = 1 Y n + d f f )= 0
Equations (4) and ( 5 ) represent a linear two-point boundary-value problem. Two different approaches have been used to solve this problem in Chapter 3. T h e first one is based on the superposition principle; and the second one is the finite-difference method. Both of these approaches have limitations. Because the rapid-access memory in current computers is limited, the finite-difference method is limited to problems with small number of variables. T h e problem of ill-conditioning in solving the algebraic equations can cause considerable difficulties when the first approach is used. Invariant imbedding can be used to form the third approach. This approach avoids the numerical solution of algebraic
4. LINEAR
225
BOUNDARY-VALUE PROBLEMS
equations. However, both this and the first approaches may have stability problems caused by the marching integration technique. T o use the invariant imbedding approach, consider again the problem with the more general boundary conditions (64
%+,(a) = c
YTZ,l(tf) = 0 (6b) An expression for the missing initial condition for the system represented by (4)and (6) can be obtained by considering a family of processes with different starting points, a. Let Yn+1(4 = rTZ+,(c, 4
(7)
+
where Y,+~ represents the missing initial condition for the (n 1)st quasilinearization iteration. We can obtain the invariant imbedding equation from Eq. (12.5) of Chapter 6
=
-q1(4rn+1(c,
4 - 2qz(a)xn(a)c + 4 2 ( a ) g ( 4
(8)
remembering that the variables with subscript n are known variables. T h e final condition for Eq. (8) is rn+,(c, 5) = 0
(9)
Since Eq. (4)is linear, the function r,+, can be expressed as rn+,(c, 4 = c 5 n + 1 ( 4
+
rln+1(4
(10)
Substituting (10) into (8) and equating terms involving like powers of c, we obtain
Equation (11) could have been obtained directly from Eq. (13.5) of Chapter 6. From Eqs. (9) and (lo), we see that 5n+1(tr) = 0
.ITZ+dtr) = 0
(12)
Equations (11) and (12) constitute a nonlinear initial-value problem. T h e missing initial condition for the original systems (4) and (6) can be obtained by integrating these equations.
226
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
Equations (1) and (2), which are now represented by the recurrence relations (4),( 5 ) , (1l), and (12), can be solved in the following iterative manner. With an assumed initial approximation for xnSo(t)= xnz0(u), Eq. (1 1) can be integrated backward starting at the final condition (12). Thus, we have the functions {,+,,,(a) and ~ ~ + , , ~ (for a ) all the processes with different starting points, a. We can then obtain the missing initial condition for the original process with a = 0 from the equation
remembering that the value of c for the original process is x,+,(O) = 1. Now, Eqs. (4),(5a), and (13) with n 1 = 1 become an initial-value problem and the results of the first quasilinearization iteration can be obtained by integrating these equations with the assumed initial approximation xo(t). With xl(t) and yl(t) known, Eq. (1 1) can be integrated with n = 1 to obtain the functions {n+l=2(a) and ~ , + ~ , ~ ( aand ) hence the missing initial condition, ~ , + , = ~ ( 0Thus, ) . the results for the second iteration can be obtained. This procedure can be continued until no more improvements on the values of x and y can be obtained. Notice that there are two independent variables in the above equations. T h e independent variable in systems (4)and (5) is t, while in the invariant imbedding equations (1 1) and (12) it is a. Notice also that the missing initial condition is obtained by integrating (1 1) backward starting with (12), which obeys the given final condition. Then, on the basis of this newly obtained missing initial condition, the solutions of (4)and ( 5 ) are obtained by integrating (4)forward. Thus, the calculated value of the given final condition, y,+,( tr), depends on two integration processes which are performed in series. Consequently, the present procedure is very sensitive to the accuracy of the integration processes. Obviously, this procedure is not suited to problems which are unstable when the marching integration technique is used. However, since no system of algebraic equations is solved, it seems to be more suited to problems with ill-conditioned algebraic systems.
+
5. Numerical Results
T h e problem formulated above is solved with the following numerical values: t , = 1.0 4 d t ) = -6 (1) q z ( t ) = -12 At da = 0.01
5.
227
NUMERICAL RESULTS
where d t and da are the integration step sizes for Eqs. (4.4) and (4.1 l), respectively. With the initial approximation equal to the given initial condition, X,Jtk)
=
1.0
Fz
=
0, 1, 2 ,..., N
(2)
TABLE 7.5
CONVERGENCE RATESWITH INVARIANT IMBEDDING
Iteration 1
0 0.2 0.4 0.6 0.8 1 .O
1 .O 0.78880 0.66685 0.59664 0.55740 0.54222
-1.3723 -0.79251 -0.45732 -0.261 74 -0.13737 -0.1io10 x 10-4
-2.7445 -2.7442 -2.7409 -2.7083 -2.3933 0.0
1.3723 1.3721 1.3705 1.3541 1.1967 0.0
Iteration 2
1 .O
1 .O 0.77547 0.62838 0.52591 0.44961 0.38042
- 1.4091 -0.89141 -0.60511 -0.43401 -0.34280 -0.38498
0 0.2 0.4 0.6 0.8 1.o
1 .O 0.77537 0.62807 0.52608 0.45633 0.42546
- 1.4093
0 0.2 0.4 0.6 0.8 1 .O
1 .O 0.77537 0.62807 0.52608 0.45634 0.42559
- 1.4093
0 0.2 0.4 0.6 0.8
1 .O 0.77537 0.62807 0.52608 0.45634 0.42560
- 1.4093
0 0.2 0.4 0.6 0.8
-2.5155 -2.1346 - 1.9007 - 1.7281 - 1.4025 0.0
1.1064 0.76436 0.59321 0.49566 0.38497 0.0
Iteration 3 -0.89226 -0.60591 -0.42483 -0.26931 -0.17480
-2.5003 -2.0697 -1.7571 - 1.4901 - 1.0994 X
lo-'
0.0
1.0909 0.71251 0.49766 0.35904 0.23208 0.0
Iteration 4 -0.89226
-0.60591 -0.42481 -0.26921 0.38191 x lo-'
-2.5002 -2.0695 -1.7583 - 1.5015 - 1.1418 0.0
1.0908 0.71235 0.49838 0.36489 0.25083 0.0
Iteration 5
1 .O
-0.89226 -0.60590 -0.42481 -0.26919 0.38776 x lo-'
-2.5002 -2.0695 - 1.7583 -1.5016 -1.1419 0.0
1.0908 0.71235 0.49839 0.36490 0.25087 0.0
228
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
with t,,, - t, = A t , the results listed in Table 7.5 are obtained. The Runge-Kutta integration scheme is used for both Eqs. (4.4) and (4.11). Since y,(t) does not appear in Eq. (4.4), the initial approximation yn=O(tk) is not needed. Each iteration takes approximately less than half a second to compute on the IBM 7094 computer. Table 7.5 shows that the functions 5 and 7 at a = 0.8 have the slowest convergence rates. Furthermore, it shows that the accuracy of y ( t ) at t near tf is not very good. T h e value of y(tf) should be zero, but a value of 0.38776 x is obtained. Thus, it seems that the accuracy of the solution for x and y is not as good as the results obtained by the superposition principle. With the numerical values listed in Eq. ( l ) , this problem is the same as that solved in Section 3 of Chapter 6 or Section 2 of Chapter 7 for c = 1 and a = 0. Thus, Table 7.4 can be compared with Table 7.5. T h e accuracy can be improved by a reduction in the step size, A t , used. T h e influence of the initial approximation on the convergence rates is shown in Table 7.6. Only the missing initial and missing final conditions are shown. It can be seen that the convergence rates are almost independent of the initial approximations within the range tested. No convergence has been obtained with xo(t,) = 5, k = 0, 1, 2, ..., N . Except for the differences in the given initial conditions, the present problem is the same as that solved in Section 5 of Chapter 3. Thus, approximate comparisons also can be made with the results obtained there. That convergence has been obtained with xo(t,) = 10 in Chapter 3 seems to suggest that the present approach has a smaller interval of convergence. I n order to avoid interpolations, the two step sizes, A t and d a , should be equal. Equation (4.11) must be integrated backward. It can be TABLE 7.6 AS
CONVERGENCE RATESOF x(1) AND y ( 0 ) FUNCTIONS OF INITIAL APPROXIMATION, INVARIANT IMBEDDING
W0(tk) =
Iteration
x(1)
0.0 y(0)
XO(tk) =
0.01
XO(tk)
=
0.1
XO(tb) =
2.0
41)
Y(0)
41)
Y(0)
41)
Y(0)
0.96751 0.46717 0.38156 0.42560 0.42559
-0.03945 -1.3748 -1.4091 - 1.4093 - 1.4093
0.74046 0.14186 0.41168 0.42561 0.42559
-0.35713 -1.3919 -1.4093 - 1.4093 - 1.493
1.0 0.54222 0.38042 0.42546 0.42559
0.0 -1.3723 -1.4091 - 1.4093 -1.4093
~ _ _ _ _
~
1 2 3 4 5
1.0 0.54222 0.38042 0.42546 0.42559
0.0 -1.3723 -1.4091 -1.4093 -1.4093
6.
OPTIMUM TEMPERATURE PROFILES IN TUBULAR REACTORS
229
shown easily that the same Runge-Kutta formula can be used for backward integration,. provided that the integration step size d a is replaced by --da.
6. Optimum Temperature Profiles in Tubular Reactors
In order to test the combined quasilinearization and invariant imbedding procedure and also to compare it with other approaches, the problem solved in Section 3 of Chapter 5 will be solved by this approach. T h e equations for the state variables and Lagrange multipliers are (see Section 2 of Chapter 5 )
dhl = k,h,
dt
-
k,X,
with boundary conditions X,(O)
= x;
T h e control variable T ( t ) is obtained from the following algebraic equations: XlklEl(h1
-
A,)
+ x,k,E,X,
=0
(3)
and
Equations (1)-(4) represent the system to be solved. T h e above equations have been linearized in Section 2 of Chapter 5, by the generalized
230
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
Newton-Raphson formula. These linearized equations have been listed in (2.17) of Chapter 5 and can be represented symbolically as follows:
d xl.n+l 4 1 i . n ( t ) ~ i . n + i+ qi2.n(t)Xzmn+i + ~ i i . n ( t ) L n + i dt
dxz,n+l
-dt
+ -
qZl,n(t)Xl.n+l
+
+
+ hz,n(t)
+
vlz.n(t)&,n+l
4 .n+l - ~ z l . r r ( t ) ~ l . n + l -dt
+
+
(54
+
q Z ~ . n ( t ) ~ ~ . n + l~ l , n ( t ) h l , n + l
Pzz,n(t)hz,n+1
4 .n+l - Wll.n(t)Xl.n+l -dt
+
+h l M
P1z,n(t)hz.n+1
“zz.n(t)~z,n+l
wlZ.n(t)X*,nil
+
(5b)
+
Vll.n(t)l,.n+l
(5 4
~l,n(t)
Wzz.n(t)Xz,n+i
+ %(t)
+
Wz1,n(t)Ln+1
(54
with boundary conditions xl,n+l(o)= x;
(64
4
(6b)
XZ.,+l(O)
=
hl,n+l(G)
=0
(64
1 (64 The control variable T(t) has been eliminated from the differential equations by using Eqs. (3) and (4) before linearization. The coefficients q, p , w , v , h, and u are functions of the variables x ~ , x~ ~, , Al,n, ~ , A,,n, which are known functions and are obtained from the previous nth iteration. The exact forms of these coefficients can be obtained by comparing Eq. (5) with Eq. (2.17) of Chapter 5. The system represented by Eqs. (5)and(6)has been solved in Chapter 5 by the superposition principle. This system also can be solved by obtaining the missing initial conditions by invariant imbedding. Thus, the original boundary-value problem becomes an initial-value probiem that can be solved by marching integration techniques. T o use the invariant imbedding approach, the more general boundary conditions hz,n+l(tf)
=
Xl.ra+l(4
= c1
(74
Xz,n+1(4
= CZ
(7b)
6.
OPTIMUM TEMPERATURE PROFILES I N TUBULAR REACTORS
23 1
< <
with a t tf again can be considered. Equations ( 5 ) and (7) are essentially the same equations as those listed in (13.1) and (13.2) of Chapter 6 with M = 2. T h e h in Eq. (5) corresponds to y in Eq. (13.1) of Chapter 6. If 4*n+l(a) = T1.n+&1 > c2 > 4 (84 Az*n+1(4 = Tz.n+1(C1
3
c2
9
4
(8b)
the invariant imbedding equations are
Equations (9) and (10) are obtained from Eqs. (13.5) and (13.6) of Chapter 6. T h e functions 5 and 77 are defined by Tl*n+l(Cl
9
CZ >
TZ*n+l(Cl
t
CZ
4 = C1511.n+1(4 4 = C1521.n+1(4
+ +
C251z.n+1(4 C2522.n+1(4
+ +
%,@+1(4
(114
rlz,n+1(4
(Ilb)
With the nth iteration known, the missing initial conditions for the (n 1)st iteration, Al,n+l(0)and A2,n+l(0), can be obtained by integrating (9) backward starting at the final conditions, Eq. (10); and by using Eqs. (8) and (11) at a = 0. Once these two missing initial conditions are obtained, the results for the (n 1)st iteration can be obtained by integrating ( 5 ) numerically. This iterative process can be continued.
+
+
232
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
7. Numerical Results
Using the equations and procedures discussed in the preceding section, some numerical experiments have been performed. This procedure will not converge if the numerical values listed in Eqs. (3.1)-(3.3) of Chapter 5 are used with (1)
da = d t
The values of the functions <(a) and ~ ( a increase ) rapidly during the numerical integration of the first iteration and extremely large values for these functions have been obtained at a = 0. However, convergence can be obtained if the initial approximations are obtained by using the following assumed temperature profile: T,(t,) = 335°K
k
= 0,
1, 2 ,...,N
(2)
These initial approximations have been obtained by first integrating Eqs. (6.la) and (6.lb) with (6.2a) and (6.2b) as the initial conditions and with T(tk)given by (2). T h e values of the Lagrange multipliers are then obtained by integrating Eqs. ( 6 . 1 ~ and ) (6.ld) backward with (6.2~)and (6.2d) as the final conditions. With the initial approximations obtained in this manner and with the numerical values listed in Eqs. (3.1) and (3.2) of Chapter 5 , the results listed in Table 7.7 are obtained. A figure very similar to Fig. 5.2 would be obtained if the convergence rate of the temperature profile were plotted. T h e optimum values for the functions and 77 are shown in Figure 7.1. Table 7.7a should be compared with
<
a
FIG.7.1. Optimum values of 1 and 7 , tt
=
8.
7.
233
NUMERICAL RESULTS
Table 5.2. Notice that the convergence rates for these two approaches are almost the same. However, since no convergence was obtained when the given boundary conditions were used as the initial approximations, the present procedure seems to have a smaller interval of convergence than the procedure used in Section 3 of Chapter 5. TABLE 7.7a CONVERGENCE RATESWITH t, Iteration
0 1 2 3 4
XI@,)
0.21 120 0.17760 0.17058 0.17044 0.17043
X&f)
0.67248 0.67737 0.67935 0.67943 0.67943
=
8
AND
To&) = 335°K
4(O)
UO)
T(O)
T(tf)
0.55788 0.60293 0.60946 0.61000 0.61000
0.87629 0.84683 0.82873 0.82823 0.82823
335.00 342.55 340.90 340.82 340.82
335.00 337.06 335.94 335.89 335.89
TABLE 7.7b
CONVERGENCE RATES
1 2 3 4 5
0.11490 0.12861 0.13607 0.13590 0.13591
WITH
-0.15062 -0.16847 -0.16547 -0.16514 -0.16515
tf
=
8
AND
-0.15062 -0.16847 -0.16547 -0.16514 -0.16515
TO&)= 335°K
0.35652 0.22436 0.20126 0.20072 0.20072
0.60680 0.61374 0.60903 0.60898 0.60898
0.77336 0.82154 0.82939 0.82945 0.82945
Some numerical experiments also have been performed with the more severe reaction conditions listed in (3.5) of Chapter 5. No convergence was obtained when the initial approximations were sbtained from the assumed temperature profile listed in (3.7) of Chapter 5. I n fact, even with the better initial temperature profile obtained from Eqs. (3.8) and (3.9) of Chapter 5, convergence still cannot be obtained. To obtain convergence, a smaller integration step size is needed, and thus a better accuracy in integration is obtained. Instead of Eq. (3.5) of Chapter 5, the following numerical values are used: t,
=
10 minutes
x! = 0.95 mole/liter
x i = 0.05 mole/liter
A t = 0.01
(3)
234
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
T h e initially assumed temperature profile is obtained by using Eqs. (3.8) and (3.9) of Chapter 5 ; and the initial approximations are obtained from this initially assumed temperature profile by integration. Use of the reaction rate constant listed in (3.1) of Chapter 5 gives the results shown in Table 7.8. T h e optimum values of the functions 5 and 7 are shown in Figure 7.2. Table 7.8a should be compared with Table 5.5. T h e conver-
0
2
4
6
10
8
a
FIG.7.2. Optimum values of 1 and 7,tt
=
10.
gence rates for these two procedures, again, are quite similar. However, Table 7.8 was obtained with a much smaller integration step size than that used in obtaining Table 5.5. T h e obtained value of A,($), which is Since this value should be zero, not shown in Table 7.8, is 0.23 x an accuracy of 0.23 x only is obtained at this final point. A much TABLE 7.8a
CONVERGENCE RATESWITH t,
0 1 2 3
4 5
0.16542 0.17709 0.17253 0.17265 0.17264 0.17264
0.67805 0.68069 0.68025 0.68012 0.68012 0.68012
0.67630 0.67704 0.67846 0.67859 0.67858 0.67859
=
10
0.71152 0.70862 0.70938 0.70948 0.70948 0.70948
345.00 360.44 359.96 359.94 359.94 359.94
335.00 335.72 336.07 336.11 336.1 1 336.11
8.
235
DISCUSSION
TABLE 7.8b
CONVERGENCE RATESWITH
1 2 3 4 5
0.85653 0.13502 0.11586 0.11549 0.11551
t, = 10
x -0.82157 x -0.82159 X 0.45777 x lo-' -0.20195 x lo-' -0.20195 X lo-' 0.36437 x lo-' -0.20494 x lo-' -0.20494 x lo-' 0.36339 x lo-' -0.20458 x lo-' -0.20458 X 10-l 0.36325 x lo-' -0.20460 x lo-' -0.20460 X lo-' 0.36326
0.66894 0.67819 0.67851 0.67851 0.67851
0.68652 0.71035 0.71078 0.71076 0.71076
higher accuracy has been obtained by using the superposition principle even with a much larger integration step size (see Table 5.4). No convergence can be obtained if the following values are used to replace the values listed in Eq. (3.8) of Chapter 5: To(0)= 360°K
T o ( t f )= 335°K
(4)
The other numerical values used remain the same. Notice that Figs. 7.1 and 7.2 are very similar. In fact, if the right-hand portion of Fig 7.2 for a = 2 to a = 10 is compared with Fig. 7.1, we find that they are almost the same. This is not surprising. T h e initial conditions for the problem with tr = 8 are xIo = 0.53 and x20 = 0.43; while from the problem with tf = 10, we obtain x1 ( t = 2) = 0.55255 and x2 ( t = 2) = 0.41405 (see Table 5.4). Thus, the problem with tf = 8 is almost the same as the problem with tf = 10 except that the former problem starts 2 minutes later. T h e initial temperature profile listed in Eqs. (3.8) and (3.9) of Chapter 5 has been estimated from the optimum temperature profile of the problem with tr = 8. T h e Runge-Kutta integration scheme has been used for all the numerical integrations discussed in this section.
8. Discussion
T h e combined quasilinearization. and invariant imbedding procedure is very sensitive to the numerical errors of the integration procedure. Thus, a much smaller integration step size is needed when this procedure is used. When the difficulties in solving the problem lie in the numerical integration process, this procedure is not suitable. On the other hand, it should be used when the difficulties in solving the problem are connected with the numerical solution of the algebraic equations.
236
7. QUASILINEARIZATION
A N D INVARIANT IMBEDDING
9. Dynamic Programming and Quasilinearization-I
For optimization problems, dynamic programming also can be used to form various combinations of numerical procedures [l, 3, 61. Quasilinearization and dynamic programming can be combined in various ways. In the first place, dynamic programming can be used to replace invariant imbedding in the predictor-corrector formula discussed in Section 2. A search procedure must be added when dynamic programming is used. Second, the original problem can be divided into a series of subprobbms. Quasilinearization can be used to obtain the optima of these subproblems first. Then dynamic programming can be used to obtain an approximate solution of the original problem by combining these subproblems in an optimum manner. To illustrate the procedure, consider the minimization of the integral
J
=
J t ’ f ( x , x’) dt 0
with the terminal conditions x(0) = cg
(24
4 t f ) = Cf
(2b)
Instead of solving the above problem directly, let us consider the problems with the shorter duration
(3) with t, = 0 and tN = tf . In other words, we have divided the original problem into (N - 1) subproblems. Now, let us see what the terminal conditions of these subproblems are. The terminal condition at t, is given by (2a). However, the terminal condition at t, is unknown. A series of terminal conditions will be allowed at this point. Suppose the value of x at t, can be reasonably assumed to be within the range xa
< <
xb
(4)
Then we can assume m terminal conditions at t, with the ith terminal condition as q t * ) = x,
+- i
x b - xa
m-1
i
= 0,1,2,..., ( m - 1)
9.
DYNAMIC PROGRAMMING AND QUASILINEARIZATION-I
237
Thus, in the interval [tl = 0, t,], m optimization problems must be solved by the quasilinearization technique. Each problem is represented by Eq. (3) with k = 1. T h e initial conditions for these problems are given by (2a) and the final conditions are given by ( 5 ) with i = o , 1 , 2 ,..., ( m - 1 ) . T h e next subproblem with the interval [t, , t,] can be treated in the same manner except that now both the initial and final terminal condi) be tions are unknown and all the possible values for x ( t z ) and ~ ( t , must considered. This procedure can be continued until the last subinterval in which the initial terminal is unknown. But, the final terminal is known and is given by Eq. (2b). Once the optima for all the subproblems are obtained, the optimum of the original problem with the interval [0, tr] can be obtained easily by dynamic programming. Consider the first subinterval, From co , we can go to any of the m points at t, . However, only one of these points gives the minimum value of (1). From these points at t , , we can go to any of the points at t, . This process can be continued. T h e problem is to find the best way to go from co to ct . This is a special case of the general routing problem which is especially suited for dynamic programming and has been considered in detail in various publications [4, 51. T h e functional equation of dynamic programming can be established easily. Let hidtk)
=
I
the minimum of
Jk with terminal conditions x t i l ( t k )and ~ ‘ j ) ( t ~ + ~ )
and gi(‘?k)
=
+
1
+ +
the minimum of the sum h,(t,) hil(tk+l) ... h i j ( t ~ - l ) ) where the process begins at tk with the starting terminal diI(tk)j
Note that the subscripts i and j in the definition of g are optimum values at times tk+l , t,,, ,..., t N - , . At time t, , the value o f j is optimum, but the value of i is not. On the other hand, at time t N P l ,the subscript i is optimum, but j must correspond to the given terminal condition ct . T h e principle of optimality yields the functional equation gdtk) = mp[hij(tk)
fg j ( t k + l ) ]
=
2,***,( N - 2,
(6)
T h e optimum of the last subinterval can be obtained by considering a process with duration t = t,-l to t = t , , gi(tN-1)
= hidtN-1) = hicf(tiv-l)
(7)
238
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
Thus, the optimum of the original process can be obtained by minimizing Eq. (6) recursively starting at the last subinterval represented by Eq. (7). This is essentially an optimization problem with N - 1 stages. The optimum obtained is, in general, not the true optimum, but only an approximate one. This approximate result can be used as the initial approximation and the quasilinearization procedure, again, can be used to obtain better results. 10. Discussion
As has been shown in Section 12 of Chapter 2, there are two ways to ensure convergence of the quasilinearization procedure. The first is to obtain an initial approximation which is sufficiently close to the desired solution; and the second is to make the interval of interest of the independent variable [ t o ,tr] sufficiently small. The simplest way to make this interval small is to consider a process with shorter duration first, and then use the results of this process to solve a process with longer duration. As has been discussed in Section 7, this approach has been used to obtain the initial temperature represented by Eqs. (3.8) and (3.9) of Chapter 5 from the results of the process of shorter duration with tf = 8. The procedure discussed in the preceding section constitutes a more sophisticated way to reduce the interval of interest of the variable. It is true that this procedure is fairly time consuming when x is multidimensional. However, if only an approximate initial approximation is needed, fairly small values of N and m can be chosen and hence the time required can be reduced considerably. 11. Linear Differential Equations
We have seen in previous chapters that the quasilinearization procedure is a useful tool for solving nonlinear boundary value problems. In the next section it will be shown that this procedure also can be used to reduce the dimensionality of an optimization problem when dynamic programming is used. However, before we discuss this approach, some well-known results concerning the solution of linear differential equations will be summarized. Consider the linear differential equations with variable coefficients
12.
DYNAMIC PROGRAMMING AND QUASILINEARIZATION-I1
239
I n vector-matrix notation, Eq. (1) becomes
with initial conditions x(0) = xo
(3)
where Q ( t ) is an M x M matrix, and p(t) and x are M-dimensional vectors. T h e solution of (2) is [7, 81
+ 1 X(t)X-'(s)p(s) ds t
x(t) = X(t)xo
0
(4)
where X ( t ) is an M x M matrix and is the solution of the matrix equation dX dt
- = Q(t)X
with initial conditions X(0)
where I is a unit matrix. At t
=
=
I
tr , Eq. (4) becomes
where c
=
X(tf)XO
K(s)
=
X(tt)X-'(s)
(7)
and c is an M-dimensional constant vector and K(s) is an M x M matrix. Equation (6) can be rewritten as
12. Dynamic Programming and Quasilinearization-I1
Now we wish to show that by the combined use of dynamic programming and quasilinearization a reduction in the dimensionality of an optimization problem can be obtained. Consider the problem of maximizing the function +I(b),
XZ(tf)!..., x m ( 5 ) )
(1)
240
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
over the control variables z(t), which are related to the state variables x by means of the nonlinear differential equations dxi -=ffi(x,z) dt
i = 1 , 2,..., M
(2)
1, 2, ..., M
(3)
with initial conditions X i ( 0 ) = x:
i
=
<
where m < M , 0 \< t tf , and x and z are M-dimensional vectors. I n addition, the problem must satisfy the constraints zi.min
< z i ( t ) < zi,max
i
= 1, 2,.-, M
(4)
For simplicity, the number of control variables has been assumed to be equal to the number of state variables. Obviously, the procedure to be discussed can also be applied if this is not the case. If the functional equations of dynamic programming are used, this problem involves the computation and storage of functions of M variables. However, Bellman and Kalaba [9, 101 have shown that if Eq. (2) were linear, the above problem could be treated by dynamic programming involving sequences of functions of m variables. Since m is equal to one or two for a number of significant problems, this is a very important reduction in the dimensionality of the problem. This is especially true in view of the fact that Eq. (2) can be linearized by the quasilinearization procedure. Thus, an iterative scheme can be formed by the combined use of dynamic programming and quasilinearization. Applying the generalized Newton-Raphson formula, Eq. (2) can be linearized
I n vector-matrix notation, Eq. (6) becomes d%+l dt - f (xn 9 zn)
+ J(xn
zn)(xn+i - xn)
+ Jz(xn
>
zn)(zn+1
- zn) (7)
where Jz(x,, z), is the same Jacobi matrix as that in Eq. (17.10) of Chapter 2, except that the differentiation is with respect to z not x.
12.
24 1
DYNAMIC PROGRAMMING AND QUASILINEARIZATION-I1
I n obtaining Eq. (7), we have considered both x and z as unknown functions. Comparing Eqs. (11.2) and (7), we obtain Q(t>
=
J(xn
(8)
2,)
+
2,) - J ( x n G)X, JAxn > ~ n ) ( ~ n + l2,) (9) Using Eqs. (8) and (9), the solution of (7) can be represented by Eq. (1 1.8). Introducing the Lagrange multiplier A, the problem becomes the maximization of
~ ( t= ) f (xn
J
7
= + ( ~ l ( t f )x*~ t f ) , * - *Xvm ( t r > ) -
J”
tf 0
f ( z )dt
(10)
over all z ( t ) ,satisfying Eqs. (2)-(4). Using Eq. (11.8), Eq. (10) becomes
where z,+,(t) are the unknown control variables after n iterations. If we consider X as a known parameter, the maximum value of Eq. (11) depends only on c1 , c2 ,..., c, , and tj . Furthermore, by examining Eqs. (11.5), (11.7), and (8), it can be seen that c is independent of the control variables zn+l(t). Thus, we wish to imbed the original problem with particular values of cl, c2 ,..., c, and duration tf within a family of processes in which cl, c2 ,..., c, and t j are parameters. Notice that if the explicit solution, Eq. (1 1.8), were not used, the maximum value of Eq. (10) would depend on c l , c2 ,..., cM and tf . We have reduced the number of variables from M to m. If m is equal to one or two, a feasible computational procedure has been obtained. Following the approach used in Chapter 6, we define g(c,
i
the maximum value of J where the process begins
> ~2 >.*.>
cm
2
a) = at t
=
u with starting state c1
,..., c,
.
1
Since the process is nonstationary, we have fixed the final time tf . A family of processes with different starting points a will be considered. T h e new maximization problem is
242
7.
QUASILINEARIZATION A N D INVARIANT IMBEDDING
T h e maximization is excuted by choosing the proper values of z over the interval [a, tf]. Applying the principle of optimality, we obtain the desired recurrence relation
T h e terms under the integral sign may be approximated by
M
cm
+ c hmj(a)Pj(a) j=l
a
+A)]
(16)
T o obtain the final condition for Eq. (16), observe that if the process had zero duration or a = tf , then the maximum value of Eq. (11) would be equal to zero. Thus g(c, > cz ~
. * cm * ~
7
tj) = 0
(17)
Notice that we have assumed that the duration of the process tf is divided into small intervals of A width. Let tf = AN, then a = 0, d, 24, ..., N d . Thus, Eq. (16) can be solved in a backward recursive fashion starting with the known final condition, Eq. (17), at a = tf . T h e computational procedure can now be summarized as follows: (a) Estimate a reasonable control policy ~ ~ = ~satisfying ( t ) , Eqs. (4) and ( 5 ) .
13.
FURTHER REDUCTION I N DIMENSIONALITY
243
(b) Calculate xn=,(t)from Eqs. (2) and (3), using the newly obtained values of z,=,(t). (c) Obtain z,,,(t) by maximizing Eq. (1l), using the newly obtained values of z,=,(t) and ~ , = ~ ( tT) h. e values of z,,,(t) must satisfy Eq. (4)(d) Return to (b) with n = 1. Equation (11) can be maximized by using the recurrence relation, Eq. (16), remembering that the solution of the linearized equation, Eq. (7), has been used in obtaining this recurrence relation.
13. Further Reduction in Dimensionality
Owing to the limited rapid-access memory of current computers, the above algorithms cannot be used if m is larger than three. However, for a large number of problems, a further reduction of the dimensionality can be obtained. Consider the problem of maximizing the function
Let us introduce a new state variable, xM+,(t), defined by
Differentiating Eq. (2) with respect to t , we have
T h e initial condition is X M + m = H(x(0))
(4)
If we consider Eqs. (12.2) and (3) as the system of differential equations, the objective function, Eq. (12. l), becomes a function of one variable
4 = XM+l(tf)
(5)
If the algorithms obtained in the previous section are used, a problem with a dimensionality of one is obtained. This is a significant reduction in terms of computational requirements.
244
7.
QUASILINEARIZATION AND INVARIANT IMBEDDING
14. Discussion
T h e above procedure can be generalized easily to problems with more general objective function. For example, the problem of maximizing the integral
J
=
/”f(x, 0 z) dt
(1)
can be treated by the above procedure if we introduce a new state variable, ~ ~ + ~defined ( t ) , by
with initial condition %4+1(0) = 0
(3)
T h e problem now becomes the maximization of ~ ~ + ~A( variety t ~ ) .of other forms of objective functions also can be treated by the algorithms obtained in the previous section. More discussion can be found in the references listed at the end of the chapter [ll-131. Let us consider briefly how Eq. (12.16) may be solved. I n order to obtain the values of c and K(a), the homogeneous differential equation, Eq. (1 1.5), must be solved first. Since Q ( t ) does not contain the unknowns x , + ~and z,+~ , Eq. (1 1.5) can be solved easily. Thus, in actual computations c and K(a) in Eq. (12.16) can be considered known. T h e functions f and p contain the unknown control variables Z ~ + ~ ( U T ) . h e problem is to find Z,+~(U) so that the expression inside the square bracket on the right-hand side of Eq. (12.16) is maximized. Since dynamic programming is especially suited for discrete processes, the present approach in the reduction of dimensionality appears to be a promising tool for solving stagewise processes. REFERENCES 1. Bellman, R., Kagiwada, H. H., and Kalaba, R., Numerical studies of a two-point nonlinear boundary-value problem using dynamic programming, invariant imbedding, and quasilinearization. RM-4069-PR. RAND Corp., Santa Monica, California, March, 1964. 2. Bellman, R., and Kalaba, R., Dynamic programming, invariant imbedding and quasilinearization: Comparisons and interconnections. RM-4038-PR. RAND Corp., Santa Monica, California, March, 1964; see also Bellman and Kalaba in “Computing Methods in Optimization Problems.” Balakrishnan, A. V., and Neustadt, L. W., eds., Academic Press, New York, 1964.
REFERENCES
245
3. Bellman, R., and Kalaba, R., “Quasilinearization and Nonlinear Boundary-Value Problems.” American Elsevier, New York, 1965. 4. Bellman, R., and Dreyfus, S., “Applied Dynamic Programming.” Princeton Univ. Press, Princeton, New Jersey, 1962. 5. Kalaba, R., On some communication network problems. “Combinatorial Analysis.” Am. Math. SOC.,Providence, Rhode Island, 1960. 6. Kalaba, R., Graph theory and automatic control, in “Applied Combinatorial Mathematics” (E. F. Beckenbach, ed.). Wiley, New York, 1964. 7. Bellman, R., “Introduction to Matrix Analysis.” McGraw-Hill, New York, 1960. 8. Bellman, R., “Stability Theory of Differential Equations.” McGraw-Hill, New York,
1953. 9. Bellman, R., Some new techniques in the dynamic programming solution of variational problems, Quart. Appl. Math. 16, 295 (1958). 10. Bellman, R., and Kalaba, R., Reduction of dimensionality, dynamic programming, and control processes, J . Basic Eng. 83, 82 (1961). 11. Rozonoer, L. I., The maximum principle of L. S. Pontryagin in optimal system theory, Automatika Telemekhhanika 20, 1320, 1441, 1561 (1959), [English transl. Automation Remote Control 20, 1288, 1405, 1517 (1960)l. 12. Katz, S., Best operating point for staged systems, Ind. Eng. Chem. Fundamentals 1, 226 (1962). 13. Fan, L. T., “The Continuous Maximum Principle.” Wiley, New York, 1966.