I n d . Eng. C h e m . Res. 1996,34,536-544
536
PROCESS DESIGN AND CONTROL Bracketing Method for On-Line Solution of Low-Dimensional Nonlinear Algebraic Equations Yash P. Gupta Department of Chemical Engineering, Technical University of Nova Scotia, Halifax, Nova Scotia, Canada B3J 2x4
On-line solution of system of nonlinear equations is required in several areas of chemical engineering. Failure of currently available numerical methods to converge t o a solution or convergence to a n infeasible solution creates problems. For these applications, methods are needed that would find a feasible solution within specified bounds without breaking down. Since bracketing methods offer a high probability of convergence, a n extension of the well-known bracketing method is presented in this paper. Its performance is demonstrated on 12 example problems from the chemical engineering literature. The high probability of convergence obtainable makes the proposed method a viable candidate for consideration.
1. Introduction Solution of nonlinear algebraic equations (NLE) is required in many areas of chemical engineering, such as process design, simulation, optimization, and control. A large number of numerical methods have been developed for solution of these equations. A comprehensive treatment of these methods is given by Ortega and Rheinboldt (1970) in their book. The book by Dennis and Schnabel(1983) covers most Newton-based methods. Numerical solution of constrained NLE is discussed by Shacham (1986). Mapped continuation methods for computing all solutions of NLE have been proposed by Seader et al. (1990). For polynomial equations, sofiware packages such as HOMPACK (Watson et al., 1987) and CONSOL (Morgan, 1987) are available. Reviews of the available software have been made by Hiebert (19821, Chen and Stadtherr (1981),and More and Cosnard (1979). Shacham (1985) has compared five softwares on chemical engineering problems. An increasing number of applications are requiring an o n - l i n e solution of these equations. For instance, t o handle the nonlinearities of real-life processes, the online optimization and control calculations generally require solution of nonlinear algebraic equations as a part of the overall calculations. Since numerical methods fail occasionally to converge to a solution, this causes a problem for these applications. The problem is enhanced especially when the initial guess is not close to the solution. Nonlinear equations often have infeasible solutions, and convergence to one of these solutions still causes a problem because these solutions cannot be employed. The human intervention t o fix the problem causes an unacceptable amount of delay. For such applications, methods are needed that offer a high probability of convergence to a feasible solution. The methods need to be capable of handling nonnegativity and simple bounds on the variables to obtain feasible solutions within desired bounds. On the few occasions when the methods still fail, some backup plans, such as solution of linearized equations, can be used. Bracketing methods, although somewhat inefficient computationally, offer a high probability of convergence. These 0888-588519512634-0536$09.00/0
may be employed either t o begin with or as a backup when other methods fail. The methods of bisection and regula-falsi are well-known bracketing methods for onedimensional problems (single equations). The purpose of this paper is to extend the bracketing method for single equations t o simultaneous equations and to show that it converges from a wide range of initial guesses. To keep the solution feasible and in the desired range, the method allows bounds on variables without requiring extra equations. An additional advantage is that the derivatives of the functions are not required. Although the method is recommended for lowdimensional problems, it is useful because of the number of applications that require solution of such equations. Moreover, at times it is possible to reduce the number of equations through elimination of variables. The handling of additional constraints, introduced by the elimination procedure is discussed in this paper. The high probability of convergence offered by the proposed method is useful in general. It is particularly useful for on-line applications and applications where the NLE solver is to be embedded in other large programs because the failure of the embedded program can be very costly.
2. Method for Two-Dimensional Problems Let the system of equations to be solved be represented by
Like the method for one-dimensional problems, the proposed method is simple in concept and is explained below: Let the equations f 1 = 0 and f 2 = 0, when drawn on a x 2 versus XI plot, be as shown in Figure 1. The points on these lines represent values of x 1 and x2 where f 1 or f 2 is zero. The solution of the simultaneous equations is given by the intersection of the two lines. Let us assume that f 2 is positive in the region to the
0 1995 American Chemical Society
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 537
I
f,=O
f,
f, = 0
=o
J x,
I
1‘
f,
I
N
XlP
=o
+ XI
Figure 1. Functions where the bracket on the solution can be maintained.
Figure 2. Functions where the bracket on the solution can be lost.
right of the line representing f2 = 0 and negative in the region to the left of this line. For this situation, the solution can be reached through the following steps: Step 1. Start with a guessed value of x2. Step 2. For the current value of x2, solve eq 1using a one-dimensional bracketing method. Step 3. Evaluate f2 using the current value of x2 and the value of x1 found in step 2. Step 4. Designate the value of x2 as x2p if f2 is positive or as X2N if f2 is negative. If f2 is zero, then the solution t o eqs 1 and 2 has been found. Step 5. Now repeat steps 2 and 3 starting with different values of x2 until f2 has a sign that is opposite of the one found in step 4. Designate this value of 2 2 as X2N or x2p as appropriate. Note that to find the required x2, one needs to search both sides of the starting guess for 2 2 . Step 6. Now since the required solution has been bracketed between x2p and XZN, it can be found by reducing the interval between x2p and X2N using the bisection method or linear interpolation (regula falsi) method. Note that steps 2, 3, and 4 need to be carried out for the value of x2 selected between x2p and X2N at each iteration. In the above setup, x2 and x1 may be regarded as being in the outer and inner loops, respectively. It may be noted that the above procedure will not always maintain a bracket on the solution. Other strategies need to be incorporated in the above procedure t o ensure convergence and to recover from infeasible solutions. Therefore, the following strategies were included in the proposed method: 1. Consider the situation shown in Figure 2. In this case as the interval between x2p and X2N is reduced, the bracket on the solution can be lost. The search may stop around x2s where the interval between x2p and X2N has reached a small value but the values ofx1, corresponding t o x2p and XZN, are quite different. For the situation shown in Figure 2, there are other values of x2s where the search may stop in this manner. When this happens, the bracket on the solution can be reestablished by switching x1 to the outer loop and x2 to the inner loop, that is, by reducing the interval between x1p and X1N where eq 1 is now solved for x2. This switching of variables was performed whenever the interval between the outer variables reached a specified low limit and the correspondinginner variables were quite different. A limit on the number of times
the switching was allowed was included in the program to avoid an endless cycle. A value for this limit equal to the number of equations being solved worked well on the example problems tested. 2. Since it is possible for eq 1to have no solution for certain values of x2 (or XI) being tried, provisions were made in the method to select a different x2 (or X I ) after a certain number of attempts were made. A new x2 (or X I ) was also needed whenever the solution of eq 1 did not result in the required sign for f2. In such cases, both sides of x2 (or XI) were searched between the specified minimum and maximum values. In order to find the required x2 (or XI) quickly, different searches were initiated in “parallel”, some of which moved outward to the limits of the interval and some that moved inward from the outer limits. The intervals were searched twice, using first a coarse grid and then a finer grid. 3. Provisions were made to switch the sequence of equations if the initial sequence did not work. For example instead of solving the first equation for a given value of x2, the second equation was solved. 4. An option was provided to automatically recover from infeasible solutions. When this option was selected and an infeasible solution was encountered, the total interval between the specified minimum and maximum values for the outer variable (x2) was divided into two regions around the infeasible point. The search for a feasible solution was then restarted in the two new regions. 5. The following three options were implemented in the program and studied for reducing the interval whenever the bracketing x values had been found: (i) bisection of interval; (ii) linear interpolation (regula-falsi method); (iii) quadratic interpolation. The linear interpolation generally worked the best for execution times. All the results reported in this paper were generated using this option. 6. Other commonly used strategies were implemented to prevent the program from hanging up and to speed up convergence. For example, the denominators were checked before the divisions were performed. If found to be zero, then an appropriate step was taken to keep the search going. In the linear and quadratic interpolation methods, if the new point found was very close to one of the limits of the interval, it was shifted toward the center.
538 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
3. Method for Three-DimensionalProblems
The methodology presented above can be extended to three-dimensional problems. For example consider the following equations:
S t a r t w i t h guessed x , . x, and x,
I
‘For Current x,, solve Eqns ( 3 ) & (4) for x, and x,
1
~
~
-~
of 2 - D problem in each case
Step 1. Start with a guessed value O f x 3 . Step 2. For the current value of x 3 , find the solution of eqs 1 and 2 as described in the previous section. Step 3 . Evaluate f 3 using the current value of x 3 and the values of x1 and x2 found in step 2. Step 4. Designate the current value of x 3 as x 3 p if f 3 is positive or as X 3 N if f 3 is negative. If f 3 is zero, the solution t o eqs 3-5 has been found. Step 5. Now repeat steps 2 and 3 starting with different values of x 3 until f 3 has a sign that is opposite of the one found in step 4. Designate this value of x 3 as X 3 N or x 3 p as appropriate. Step 6. Now a solution to the three equations has been bracketed. It may be found by reducing the interval between x 3 p and X 3 N using the bisection method, linear interpolation, or quadratic interpolation. Other strategies as mentioned in the previous section need to be incorporated to ensure convergence. 4. Method for Higher-DimensionalProblems
The methodology can be applied in a similar way to higher-dimensional problems. However, since the computational effort increases sharply as the number of equations being solved increases, its suitability depends upon the time that is available for finding a solution. In general, an effort should be made, whenever possible, to reduce the number of equations through elimination of variables. Note that the constraints on the eliminated variables still need to be maintained. One way to do this is to translate these to constraints on remaining variables. The other is to ignore them and then check if the solution found violates constraints on eliminated variables. If it does, the solution is infeasible and the search needs to be restarted. The first approach is better and should be used whenever possible.
5. Computer Program Since the overall pattern of calculations is the same irrespective of the number of equations involved, a program in BASIC using recursive subprograms was developed. Please note that a recursive subprogram can call itself and can call other subprograms that in turn call the first subprogram before the first subprogram has completed execution. Because of this, the size of the program could be kept small and the same program could be used irrespective of the dimension of the problem. An outline of the computer program is shown in Figure 3 through a flow chart. This flow chart is drawn for a three-dimensional problem and shows the outer loop. Essentially the same flow chart is involved when solving one equation or two equations in the inner loops for the three-dimensional problem. For the inner loops, whenever the block marked “Unsuccessful Exit” is
solution of 2-D problem in each case
Successful Exit
Unsuccessful Exit
-
Figure 3. Outer loop for a three-dimensional problem.
encountered, the program continues with a different value of the adjacent outer loop variable.
6. Testing The performance of the proposed method is presented on 12 example problems from the chemical engineering literature. Most of these problems have been used by other researchers for evaluating the performance of NLE solvers. For every case considered for each of the problems, a “rectangular” grid was chosen around the solution. The method was started using one of the grid points as the initial guess and its performance was recorded. This procedure was repeated for each of the points of the grid. This enabled us to test the ability of the method to converge from different initial guesses. The percentage of the grid points from which the method converged to a solution was noted. The bounds on the variables were the same as the bounds of the grid. To put the results obtained in perspective, the performance of CONLES and Powell-Hybrid methods on the problems was also recorded. CONLES method was chosen because it performed generally the best among the five methods compared by Shacham (1985) on chemical engineering problems. CONLES is a combination of Newton-Raphson, Marquardt-Levenberg, and continuation methods. A computer program for this method is given by Shacham (1986) and was used in this study. The Powell-Hybrid method was chosen because it is one of the standard methods for solving NLE and was available on our university mainframe. The runs for the proposed bracketing method and the CONLES method were performed on an IBM-compatible 486DX 33 MHz personal computer. The total time for the grid was noted and the average time for a run was calculated. The times for the Powell-Hybrid runs are not reported because these runs were made on the mainframe and run times were not available.
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 539 Table 1. Comparison of Convergence to Feasible and Infeasible Solutions for Example Problem Set 1 convergence (%I problem parameters
grid size
N/A
0 5 x1 5 100 0 5 x2 5 10 0 5 x3 5 100 50 5 X I 5 100 5 5 x2 5 10 50 5 x3 5 100 75 5 X I 5 85 65x258 55 5 x3 5 65 05x152
1
2
ki=31.240 krl = 2.062 k2 = 0.272 k,2 = 0.020
grid size
N/A
0 5 x 1 5 100 0 5 x 2 5 10 0 5 x3 5 100
1
34.7
16.2
100
53.7
25.5
100
69.9
48.1
100
81.0
82.6
100
= 17.721 = 3.483 kz = 0.118 k,z = 0.033
05
= 505.051 k i = 17.721 k,1 = 6.966 k2 = 0.118 krz = 333.333
0.89 5 ~2 5 1 0 5 21 5 2
5
2
96.7
2 89.3
= 17.721 3.483 kz = 0.118 kn = 0.033
100
9.7
100
21.8
17.6
100
57.4
16.7
100
38.8
43.8
100
70.2
86.0
100
34.7
51.2
100
97.5
86.0
86.0
100
0 5 ~2 5 2 0.03 5 X I 5 2 0 5 ~1 5 2
k,i=
05 52 0.05 5 ~1 5 2
0 522 5 2 0.05 5 5 2
100 = 505.051 0.89 5 x z 5 1 0 5 21 5 2 k l = 17.721 k3
100
100
100
kri
0 522 5 2
= 6.966
k2 = 0.118 k,2 = 333.333
0.05 5 21 5 2 100
100
100 k3
3 100
88.4
98.3
100
100
100
= 505.051 0.89 5 ~2 5 1
NIA
05x151
100
1
100
0 5 xz 5 2 0.05 5 5 2
100
05x151 0 5 x 4 5
7.9
= 303.030 0.79 5 2 2 5 1
ki
75.2
5 5 x2 5 10 50 5 x3 5 100 75 5 21 5 85 6 5 22 5 8 55 5 x3 5 65 0 5 XI 5 2
100
= 505.051 0.89 5 zz5 1
N/A
= 31.240 k,1 = 2.062 k2 = 0.272 krz = 0.020 ki
k3
90.9 k3
Powell CONLES bracketing
50 5 X I 5 100
0 5 22 5 2 0.03 5 X I 5 2
kri
3
problem parameters
= 303.030 0.79 5 x2 5 1
ki
k3
convergence (%)
~ _ _ _ _
Powell CONLES bracketing
85.1 k3
Table 2. Convergence to Feasible Solutions for Example Problem Set 1
66.1
43.8
100
71.1
54.5
100
05x451 0.0 5 X I 5 0.5
0.0 5 2 1 5 0.5
100
95.9
100
0.5 5 x4 5 1.0
The 12 problems are grouped into four sets to reduce the number of tables and thus to conserve space. For the two-dimensional problems, the rectangular grids consisted of eleven equally spaced values for each of the two variables. Thus each grid contained 112 = 121 initial guesses. For the three-dimensional problems, the grids consisted of six equally spaced values for each of the three variables. These grids contained 63 = 216 initial guesses. Generally, two grids of different sizes and/or two different sets of parameters are considered for each of the problems. 6.1. Example Problem Set 1. In this example set, the performance of the methods is presented on three reactor problems that have many infeasible solutions. Since the computer programs for the Powell and CONLES methods had not been designed to recover from infeasible solutions, we first compare the performance of the methods for a case where both feasible and infeasible solutions are accepted. To keep this comparison fair, the option in the bracketing method to recover from infeasible solutions was not used. This comparison is shown in Table 1. Please note that the infeasible solutions found for the reduced system of equations do satisfy these equations. These are considered infeasible because these violate constraints on the variables in the equations or on the eliminated variables. Since the infeasible solutions reached are essentially failures from a practical standpoint, one needs to recover from them. To do so, it is not enough to restart the search from another initial guess because the same infeasible solution may again be found. Note that we are already evaluating the performance of the methods when they are started from different initial guesses by scanning the grids mentioned earlier. For this recovery one needs to introduce limits in which the restarted
0.5 5 ~4
5
1.0
search must remain. This feature of limiting the search region is inherent in the bracketing method. In Table 2 we report the percentage of runs that reached feasible solutions by using the automatic recovery option in the bracketing method. The percentage of runs that reached feasible solutions using Powell and CONLES methods are also included for reference. Please note from the title of Table 2 that we are not making a comparison here. A comparison of the methods for finding feasible solutions is made later on in the discussion section by assuming perfect recovery from infeasible solutions by the Powell and CONLES methods. Problem 1. This problem considers the 1985 AIChE Student Contest Problem where the esterification reaction is carried out in two CSTR in series. The sevenequation model for this problem has been studied by Seader et al. (1990). This problem has 44 solutions and only 1 of them is feasible. The seven equations were reduced to the following three equations through elimination of variables.
5.5(0.098&1 - X2 - Xg - o . o k : , ) ( X 5
+ Xs) -
X#8
=0 (8)
where x1 is inlet flow rate (lb-mom),x2 is water in offgas from first reactor (lb-molh), x3 is BuOH in off-gas from first reactor (lb-mom), and
540 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
.x5 = 1 2 . 1 2 8 2 5 6 ~-~x2 - 1 2 h 4 - 1.752, - 75.85948
x6 = 61.177 - 8 . 7 6 1 4 ~-~X 2 X8
+ 9 h 4- + X5
= 1 . 0 9 8 6 ~-~X2 - 9X4 - X5 - x6
+
time (s)
X7
X7
The feasible solution is given by xT = (80.580, 6.976, 61.154). The methods were applied to the reduced three-dimensional problem, and their performance is shown in Tables 1 and 2. Problem 2. This problem considers the Cutlip sixequation model to represent the rate of different steps in a chemical reaction. This problem has many infeasible solutions and was presented by Shacham (1986). The six equations were reduced to the following two equations through elimination of variables. k1x1x6 - krlx4 - k#4x5 = 0
(9)
1.5(k+96 - kr95) - k g 4 ~ 5 =O
(10)
where x1 and x2 represent moles of two species and x4 = (XI
+ k1x1x6 - l)/krl
= (X2
+k996 - 114.2
X5
x3 = 2kg4X5
The reduced two-dimensional problem was solved for three sets of parameters considered by Shacham. The feasible solutions found for the first, second, and third sets of parameters are (XI= 0.97442, x 2 = 0.982831, (XI = 0.97104, x2 = 0.98037), and (XI = 0.97895, xz = 0.98605), respectively. Two grids were considered for each set of parameters. The limits of the wider grids were same as given in the problem statement. The limits of the narrower grids were derived t o keep the eliminated variables, x3 through X6, nonnegative. These limits in terms of the parameters are as follows: x1 L 141
+ kl),
x2 2 141
+ k2)
The performance of the methods is shown in Tables 1 and 2. Problem 3. This problem considers the production of sythetic gas in an adiabatic reactor. The seven equations involved were reduced to the following two equations by Kuno and Seader (1988):
+
- O . O O ~ X ~ X , ~ 0.16517~: - 0.046354~1~4 -
+
0.2070~ -~o.262O9&cl2 4- 0 . 1 7 6 9 6 7 ~ ~ 0.043528 = 0 (11)
0.h:
Table 3. Average Execution Run Time for Example Problem Set 1
+
-0.185963~~ 0 .~0~2 1 3 8 7 5 ~-~
0.086512~ -~ 0~. 0 1 6 4 1 5 ~=~0 (12)
where XI and x4 represent mole fractions. The feasible solution is X I = 0.32287 and x4 = 0.61817. The performance of the methods on the reduced twodimensional problem is shown for two grids in Tables 1 and 2.
Droblem
CONLES
bracketing
1 2 3
0.16 0.03 0.08
1.10 0.55 0.04
Discussion of Results for Example Problem Set 1. When both feasible and infeasible solutions are accepted, Table 1 shows that the performance of the Powell and CONLES methods is good except for problem 1. The convergence improves as the initial guesses are chosen closer to the solution. The bracketing method found a solution starting from all of the initial guesses. Now to make a comparison of the methods in finding feasible solutions, let us assume that a feature was available in the programs for the Powell and CONLES methods that always provided recovery from an infeasible solution to a feasible solution. In this case, the percent convergence of these methods to feasible solutions would increase from the values shown in Table 2. However, these values will only increase to the values shown in Table 1 because the recovery feature would take the infeasible solutions to feasible solutions. The percentage of runs, where neither feasible nor infeasible solutions were found, would remain the same. Comparing Table 1values for the Powell and CONLES methods with Table 2 values for the bracketing method, we find that the bracketing method provides an equal or better convergence to feasible solutions even when we consider a perfect recovery from the infeasible solutions for the other two methods. For runs in Table 2, although infeasible solutions were encountered by the bracketing method, it recovered from them. For certain runs, as many as four infeasible solutions were encountered. Since the narrowest grids for problems 1 and 2 did not contain any infeasible solutions, the bracketing method completely avoided them by staying within the grid ranges. This is the advantage of translating the constraints on the eliminated variables to constraints on the remaining variables. The average execution run times for the CONLES and bracketing methods are given in Table 3. It should be pointed out that the CONLES method allows lower limits on variables. The upper limits can be imposed indirectly by adding extra equations in the set whereby the dimensionality of the problem increases. 6.2. Example Problem Set 2. In this example set, the performance of the methods is presented on three reactor problems that have been used in control studies. For problems in this and the following sets, infeasible solutions were not encountered within the bounds of the grids. Therefore, the recovery feature of the bracketing method was not utilized. Problem 1. This problem considers an exothermic reversible reaction, A t B, in a CSTR. The dynamic model for this problem was presented by Li and Biegler (1988) and was used by Patwardhan and Madhavan (1993) to study nonlinear predictive control. The steady state equations are given by
+
- O . ~ ~ X ~ X ~ -kl(l ~ U~ x) - k p , = 0
+
(13)
O . l 6 u , ~ p ~-- O ~.l63~p~-~u~ 5(k1(1 - xl) - k g l ) = 0 (14)
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 541 0 . 1 6~~ 0 . 4 Q 5= o
(15)
Table 4. Comparison of Convergence to Solutions for Example Problem Set 2 convergence (70)
where
k, = 3 x k, = 6
problem parameters
io5 exp(-5000/x,)
u1=l ~z=300
1
05x15 1 0 5 ~2 5 1000
Powell CONLES bracketing
19.0
54.6
100
28.7
78.7
100
76.9
69.9
100
600 83.3 05x351 Q=5.1x1O6O5C~51 0 5 CB 5 1 90.3
83.3
100
82.9
100
05x3~1
x lo7 exp(-7500/x2)
and x3 represent conversion, temperature, and liquid level; and u1 and u~ represent inlet flow and inlet flow temperature, respectively. This three-dimensional problem was solved for XI, X Z , and x3 for two specified sets of u1 and u2. The solutions for the two sets of u are (XI= 0.017, x2 = 300.086, x3 = 0.160) and (XI = 0.551, x 2 = 402.757, x3 = 0.640). The performance of the methods for two grids considered for each set of u is shown in Table 4. Problem 2. This problem considers consecutive reactions, A B C, in a CSTR. The dynamic model for this problem was presented by Bequette (1991) t o study nonlinear predictive contol of temperature and concentration. The steady state component and energy balances are given by
05x151 200 5 ~2 5 600
O5z35l
XI, x2,
u1= 2 uz=400
(16)
05x15
1
0
1000
5x2 5
05x351 05x151 200 5 XI, 5
2
0 5 T 5 1000
05C*Sl 0 5 C B 1~ 100 5-T 5 200 Q = 5.2 x lo6 0 5 CA5 1 0 5 CB 5 1
- -
0.1( 1 - CA) - k1CA2 = 0
grid size
05
3
T5
93.1
100
86.6
100
82.9
100
1000
u1=l
O5C*51 0 5 C B 1~ 100 5 T 5 200 0 5 5 10
u2 = 0
0 5 x2 5 10
91.7
100
100
2.5
73.6
100
05x151
15.7 0.5
05x251 0 5 X i 5 10
u2 = 1.0
0 5 x2 5 10
U1=
100
3.3
91.7
100
05x151
where
+ 273.16))) k 2 = 172.2 exp(-34833/(8,314(T + 273.16)))
97.5
CAand CBrespresent concentration of A and B, respectively. This three-dimensional problem was solved for CA,CB, and T for two specified values of Q. The solutions for two sets of Q are (CA= 0.158, CB = 0.771, T = 153.03) and (CA = 0.157, CB = 0.763, T = 158.8). The performance of the methods for two grids considered for each value of Q are shown in Table 4. Problem 3. This problem considers consecutive reactions A B C, in an isothermal CSTR. Here the reaction A B is second order and the reaction B C is half order. The dynamic model for this problem was used by Ray (1981) in different control studies. The steady state equations are as follows:
-- -
-xl -xl 2 + u , x12 - x,
+
=o
- zx21’2 u2 = 0
(19)
(20) where x1 and x2 are normalized concentrations of A and B in reactor and u1 and u2 are normalized inlet concentrations of A and B. This two-dimensional problem was solved for XI and x 2 for two specified sets of u1 and u2. The solutions for the two sets of u are xT = (0.618, 0.031) and (0.366, 0.212). The performance of the methods for two grids considered for each set of u is shown in Table 4. Discussion of Results for Example Problem Set 2. Although the performance of Powell and CONLES methods is generally good, especially when the grids are narrowed around the solution, they fail to converge to a solution from several grid points. The bracketing
100
100
05x251
k, = 11 exp(-4180/(8.314(T
-
100
Table 5. Average Execution Run Time for Example Problem Set 2 time problem
1 2 3
(9)
CONLES 0.10
bracketing
0.06
0.08 0.02
0.03
0.16
method converged to the solution from all of the grid points. The average run times for the CONLES and bracketing method are given in Table 5. 6.3. Example Problem Set 3. In this example set, the performance of the methods is presented on three problems that were used by Shacham (1985) for comparing software for solving nonlinear algebraic equations. Problem 1. This problem represents a steady state material and energy balance on a chemical reactor. The equations are 120~ - 75k(l - X ) = 0
(21)
~ ( 8 7 3- 2‘) - 11(T - 300) = 0
(22)
where
k = 0.12 exp(1258UT - 298)/(29811) T and x represent temperature and conversion, respectively. The solution to this problem is x = 0.96387 and T = 346.16. As pointed out by Shacham (1985),this problem has a stationary point at x = 0.79 and T = 338 which makes it difficult to solve if the initial guess is below this point. The performance of the methods on this problem is presented in Table 6 for two different grids.
542 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 Table 6. Comparison of Convergence to Solutions for Example Problem Set 3 convergence (%I problem parameters grid size Powell CONLES bracketing 1
N/A
0 5 x 5 1
20.7
49.6
100
0 5 T 5 1000 0 5 x 5 1
2
57.9
76.0
100
37.2
77.7
100
54.6
98.3
100
300 5 T 5 500 0 5 x 5 1
NIA
0 5 T 5 3000 0 5 x 5 1 3
x1
Table 7. Average Execution Run Time for Example Problem Set 3 time (SI problem
CONLES
bracketing
1
0.06 0.06
0.04
2
3
0.29
Table 8. Comparison of Convergence to Solutions for Example Problem Set 4 convergence (%I problem parameters a i d size Powell CONLES bracketing
1000 5 T 5 2000 -20 5 51 5 20
=0.1
42.1
72.7 a = 0.1
-20 5
t2 5
20
-10 5
51
5
10 48.8
70.2
a = 0.4
-20 5 -10 5
-10
‘12
47.9
21.5
100
66.1
28.9
100
5 20
tl 5
5 52 5
19.8
100
14.9
95.9
100
31.4
99.2
100
100
5 t2 5 10 -20 5 51 5 20
= 0.5
2.5 100
-10 XI
0.06 0.10
100
10
100
100
10
Problem 2. The following equations represent material and energy balances for a plug flow reactor and need to be solved to find equilibrium conversion and temperature.
2.8
58.8
100
3.7
76.4
100
0
44.4
100
2.8
65.3
100
(0.91 - 0.5~)/(9.1 - 0 . 5 ~) x2/((1 - xI2Kp)= 0 (23)
+
T ( 1 . 8 4 ~ 77.3) - 4 3 2 6 0 ~- 105128 = 0
(24)
where
Kp = exp(42300/T - 24.2 + 0.17 ln(T)) x is conversion and T is temperature. The solution to this problem is x = 0.53334 and T = 1637.7. The results obtained for this problem are shown in Table 6 for two different grids. Problem 3. The following equations are involved in correlating liquid-liquid equilibrium data.
+
1/(x1x2)- 2~,G,~/(x, x2G2I32z1G1/(x2 X ~ G ,=) ~ 0 (25)
+
+
(xl - X , ) / ( X , ~ ~ , ~ ) 6 ~ , G , ~ ( l G2)/(xl
+
+ x2G2I4+
solution from all of the grid points. The average run times for the CONLES and bracketing method are given in Table 7. 6.4. Example Problem Set 4. In this example set, the performance of the methods is presented on three problems that are badly scaled andor require variables to stay nonnegative during the search. Problem 1. This example considers the Wayburn and Seader’s (1987) problem which is as follows:
+
x12 x: = 17
(27)
+ x21’2
(28)
6~1G2~(G1 - 1)/(x2 x,G1I4 = 0 (26) 23c11’3
where Gi = exp(-azi) x1 and x2 are mole fractions and xz = 1 - x1.
These equations need to be solved for unknown parameters zl and t2 for specified values of x1 and a and are extremely nonlinear. This problem was solved for two sets of x1 and a, and the performance of the methods is shown in Table 6. The solutions for the two sets of parameters are ( z ~= -5.926, z2 = 8.690) and (z1 = z2 = 1.044). Discussion of Results for Example Problem Set 3. The problems in this example set were quite difficult to solve, as the Powell and CONLES methods did not converge to a solution from almost half of the grid points. The bracketing method again converged to a
=4
The variables, x1 and x 2 , must stay nonnegative during the search so these can be raised to the fractional powers involved. The solutions to this problem are xT = (1, 4) and (4.07, 0.65). The performance of the methods on this problem is given in Table 8. Problem 2. This example considers the Hiebert (1982) chemical equilibrium problem which consists of six equations. These equations are very badly scaled, and all six variables must be nonnegative at the end of the search because these represent amount of chemicals. The six equations were reduced to the following two equations through elimination of variables. (0.1/1.1)(0.001- x4) - 104(x4f
+
“6)3t4
=0
(55 - x 6 ) - 55 x iOl4(x4 x6)xs = 0
(29) (30)
Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995 543 Table 9. Average Execution Run Time for Example Problem Set 4
where x1 = (0.1/1.1)(0*001- x4) 22
time (e.) problem 1 2 3
= (1/1.1)(0.001 - x4)
CONLES
bracketing
0.06 0.03 0.04
0.30 0.10 3.88
x9 = 2.155 x 10-4xlx30~5~~5/x4
To keep the eliminated variables nonnegative, the limits on the variables are x4 I 0.001 and X6 I 55. The solution to this problem is x4 = 90.909 x and X6 = 1.099 x 10-lo. Because the magnitude of the variables is very small, x4 and X6 in eqs 29 and 30 were replaced by larger variables x4' and X6', respectively, using the following definitions. x i = x4 x
lo6 and x i = x6 x lo6
The two equations were solved for x4' and X6', and the performance of the methods is shown in Table 8. Problem 3. This example considers the Hiebert (1982) chemical equilibrium problem which consists of 10 equations. There are scaling problems, and all variables must be nonnegative at the end of the search. In addition, some of the variables must be maintained nonnegative throughout the search because their square roots need to be calculated. These equations were reduced to the following three equations through elimination of variables. 2xl + x 2
+ x4 + x 7 + x S + xg + 2xlo - R = 0
(31)
x1x5 - 1.93 x ~ O - ' X $ ~ = 0
(32)
(33) i=l
The retained variables in the above equations were x4, x5, and x10. These were selected so that the nonnega-
tivity constraints on the eliminated variables could be easily maintained. The other variables appearing in the above equations were calculated from the following expressions whenever needed.
T = ~ ~ d ( 3 . 8 4x 6 x1 = 3 - x4 x3 = 2R - x5/2
x7 = 3.448 x 10-3x10.5p.5
xs = (1.799/3.846)xcp,dx4
To keep the eliminated variables nonnegative, the limits on the variables are x4 I3, x5 5 4, and x10 I0.5R. The reduced three-dimensional problem was solved for two values of parameter R that were considered by Shacham (1985). The solutions for R = 10 and R = 40 are (x4 = 0.119895, x5 = 0.031741, x10 = 0.001044) and (x4 = 0.002364, x5 = 0.000604, x10 = 0.0044991, respectively. The performance of the methods is shown in Table 8. Discussion of Results for Example Problem Set 4. In this example set, the Powell method did not work well on problems 1 and 3. The CONLES method did well especially from the narrower grids around the solution. The bracketing method again converged to a solution from all of the grid points. The low convergences of the Powell method for problems 1 and 3 were mainly due to the variables getting negative which could not be raised to the fractional powers involved. The CONLES and bracketing methods did not encounter this problem as they handled nonnegativity constraints on the variables. The average run times for the CONLES and bracketing method are given in Table 9. 7. Other Comments a. Since the bracketing method searches for the solution within the specified bounds, these bounds must contain the desired solution; otherwise the solution cannot be found. The solutions lying at the boundary can be found by including the boundary points in the search. b. As in the case of single equations, the bracketing method does not find singular solutions unless one of the search points happens t o fall on the solution. In those cases where no solution is found and other cases where the equations have no solution, the program exits after making a certain number of attempts that are coded in the program. c. The same program for the bracketing method was used for all the results reported in this paper; that is, its tuning parameters did not require adjustment from one problem to another. However, it is desirable t o adjust the tuning, whenever possible, to suit the type of problem to be solved. This can reduce the computational time. d. The comparisons with very good methods, such as Powell and CONLES, show that the example problems used in this study were not trivial because of the difficulties experienced by these methods. e. As may be seen from the results, the proposed method performed very well on a number of problems converging in all cases starting from a wide range of initial guesses. It allows bounds on variables without requiring additional equations. Moreover, it does not require the derivatives of the functions. f. For on-line applications, whenever the NLE solver fails t o find the desired solution, a backup plan such as a solution of linearized equations needs to be used which
544 Ind. Eng. Chem. Res., Vol. 34, No. 2, 1995
lowers the performance. The high probability of convergence offered by the proposed method reduces the need for backup plans and thus allows a higher overall performance.
8. Conclusions a. It is possible to extend the bracketing method for single equations to higher dimensional problems and to obtain a high probability of convergence to feasible solutions within specified bounds. b. The proposed method is worth considering for solution of low-dimensional systems of nonlinear equations for on-line applications and other applications where a high probability of convergence is desired.
Acknowledgment The financial support provided by the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. The assistance provided by Mr. Phillip Christie, Ms. Rhonda Crews, and Mr. Sean Jenkins in running the computer programs is appreciated.
Literature Cited Bequette, B. W. Nonlinear Predictive Control Using Multi-rate Sampling. Can. J. Chem. Eng. 1991, 69, 136-143. Chen, H. S.; Stadtherr, M. A. A Modification of Powell's Dogleg Method for Solving Systems of Nonlinear Equations. Comput. Chem. Eng. 1981,5, 143-150. Dennis, J. E.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; PrenticeHall: Englewood Cliffs, NJ, 1983. Hiebert, K. L. An Evaluation of Mathematical Software that Solves Systems of Nonlinear Equations. ACM Trans. Math. Softw. 1982, 8, 5-11. Kuno, M.; Seader, J. D. Computing All Real Solutions to Systems of Nonlinear Equations with a Global Fixed-point Homotopy. Znd. Eng. Chem. Res. 1988,27, 1320-1329.
Li, W. C.; Biegler, L. T. Process Control Strategies for Constrained Nonlinear Systems. Znd. Eng. Chem. Res. 1988,27,1421-1433. More, J. J.; Cosnard, M. Y. Numerical Solution of Nonlinear Equations. ACM Trans. Math. Softw.1979,5, 64-73. Morgan, A. P. Solving Polynomial Systems using Continuation for Engineering and Scientific Problems; Prentice-Hall: Englewood Cliffs, NJ, 1987. Ortega, J . M.; Rheinboldt, W. C. Iterative Solutions of Nonlinear Equations i n Several Variables; Academic Press: New York, 1970. Patwardhan, S. C.; Madhavan, K. P. Nonlinear Model Predictive Control Using Second-Order Model Approximation. Znd. Eng. Chem. Res. 1993,32, 334-344. Seader, J. D.; Kuno, M.; Lin, W. J.; Johnson, S. A.; Unsworth, K.; Wiskin, J. W. Mapped Continuation Methods for Computing All Solutions to General Systems of Nonlinear Equations. Comput. Chem. Eng. 1990,14, 71-85. Shacham, M. Comparing Software to the Solution of Systems of Nonlinear Algebraic Equations Arising in Chemical Engineering. Comput. Chem. Eng. 1985, 9 , 103-112. Shacham, M. Numerical Solution of Constrained Nonlinear Algebraic Equations. Znt. J . Numer. Methods Eng. 1986,23, 14551481. Ray, W. H. Advanced Process Control; McGraw-Hill: New York, 1981; p 119. Watson, L. T.; Billups, S. C.; Morgan, A. P. Algorithm 652, HOMPACK A Suite of Codes for Globally Convergent HOMOTOPY ALGORITHMS. ACM Trans. Math. Softw. 1987, 13, 281-310. Wayburn, T. L.; Seader, J. D. Homotopy Continuation Methods for Computer-Aided Process Design. Comput. Chem. Eng. 1987, 11, 7-25.
Received for review May 18, 1994 Revised manuscript received September 28, 1994 Accepted October 5, 1994@ IE940320P Abstract published in Advance ACS Abstracts, December 15, 1994. @
