- Email: [email protected]
Simulation and analysis of complex distillation columns using continuation methods
SIMULATION AND ANALYSIS OF COMPLEX DISTILLATION COLUMNS USING CONTINUAnON METHODS Kyung Taek Han*, Kang Iu Lee, En Sup Yoon *Hanyang Chemicals Central Research Center,Taejeon Yuseong Sinseong 6,KOREA Department of Chemical Engineering, Seoul National University, Seoul 151-742, KOREA
ABSTRACT This paper presents an efficient homotopy continuation method to solve a nonlinear eqaution system which describes chemical engineering problems . The characteristics of the Newton homotopy function is studied and parametric continuation is proposed to study steady state multiplicities for a constrained system. An equation based distillation simulator employing the proposed algorithm was implemented. The proposed algorithm and simulator were tested on many numerical problems and complex distillation columns. KEYWORDS Homotopy continuation, parametric continuation, solution, bifurcation
complex distillation column,
multiple
INTRODUCTION Recently, homotopy continuation methods have been used successfully to solve nonlinear equations due to their global convergent characteristics and their ability to find multiple solutions from one starting point. This paper shows how to improve the computational efficiency and capability of tracking homotopy path. This was achieved mainly by reducing the number of function evaluations and by avoiding segment jumping. The proposed algorithm had been applied to many complex distillation columns. A few workers have used the homotopy continuation method to get multiple steady state solutions for linked distillation columns(Chavez et al. ,1986, Lin et al., 1987). Parametric continuation followed by homotopy continuation was used to investigate the steady state multiplicities of the complex column and it was confirmed that this parametric continuation is quite useful in study design and operability limitations of the system. Features of the algorithm include: 1. The Euler prediction and Newton correction to track the path , 2. Step size control using determinant monitoring and curvature, 3. Flexible convergence criteria , 4. Capability of parametric continuation followed by homotopy continuation . ALGORITHM Newton Homotopy Function 5479
European Symposium on Computer Aided Process Engineering-2
S480
For solving the nonlinear equation system, f(x)=0, the Newton homotopy function h(x,t) is defined as (1)
h(x,t) = f(x) - (I-t)f(xo),
where x ERn, t e [0, I] and x = xo at t = O. If we differentiate h(x,t) by arclength, s, and introduce the definition of the arclength, we get the following IVP(initial value problem),
[
2.!!... (ws) aWT dw (Sk) ds
]
dw (Sk) ds
= e n+ 1
•
(2)
where w = [x,t], k = continuation step number and en+ 1 is a column vector of which the n+ It h column equals one and the others are zero. Any types of predictor-corrector schemes can be used to solve the IVP. simple Euler predictor and Newton corrector in this work.
We use the
Step Size Control In this study, Allgower's algorithm (1980) based on curvature and error contraction and Choi's (1990) determinant monitoring scheme were used. The determinant of the augmented matrix approachs zero as the solution does to another path(Allgower et a l. , 1980). The negative determinant value means that the path has already passed the bifurcation point and is tracking another path. At first, (k+ l)-th prediction is carried out with a step size calculated as l\sk+1 = 9id/K, where £lid is the ideal turning angle and K is the approximated curvature. Then, determinant is checked at this predicted point by equation (3). If this determinant doesn't satisfy equation (3), the proposed step size is reduced by a ratio of 0.5. This step size control strategy is essential in tracking the path by preventing segment jumping.
I detAk+I - detAk I Max ( IdetAlk+1 , IdetAlk ) < £det where A is Jacobian of the augmented matrix and
(3)
e det is tolerance for determinant
moni toring.
Convergence Criteria Each corrected point on the path is not affected by the error of the previous step (Allgower et al.,1980). To improve computational efficiency, we can use larger convergence tolerance during the intermediate path tracking than near target point(Le. t= I) or turning point. The following criteria provide flexible convergence tolerance and so increase computational efficiency.
II h(Wi+lk+ 1 ) II < Min(£mu,£min)
(4)
£mu = £tol/(det A(wo k+ I»lIn
(5)
If ( I det A(wO k+ l) I
(6)
Else
£min
=
< I ) £min = £cor I det A(wo k+ l) I
£cor I det A(wo k+ l)l/n I
European Symposium on Computer Aided Process Engineering-2
S481
5r-;;:=================::::::;'1 1- (3.5.3.0) - (3.5.2.7) _._- (3.5.2 .4) Solution I 0
3
·1
·3
·3
3
·1
5
X1
Fig. 1. Variation of homotopy path with different starting points
Parametric Continuation To solve the nonlinear equation expressed by f(x,p) = 0, we first set the parameter in order to satisfy degree of freedom. We can get a solution curve by changing the parameter p. This method is refered as parametric continuation and the difference from the homotopy continuation is that all points on the path tracked are solutions of the nonlinear system. To start the parametric continuation we need a solution at an initial point. The starting point for parametric continuation can be easily obtained by homotopy continuation . In order to perform parametric continuation on a selected parameter, pe, following Newton homotopy function defined in Rn~2 space can be obtained as, h(x,pc ,t)
=
[f(X,ps) ] pc - C
O
- (l - t) [f(X ,ps) ] pcD - C
=
0
(7)
where ps is a known parameter and c is a starting value of parameter pe. We can find a solution at t= I by applying homotopy continuation to the above equations . Parametric continuation on pc is proceeded by equation (8). It is regarded as homotopy continuation where pc is pameterized by t as given in equation (9). h(x,pc,t) t
=
where
=
[f(X,ps) ] pc - Ct
- (l - t) [
0
C - Ct
]
(c - pc)/(c - ce) , Ct
=
0
(8) (9)
is the given target value of parameter pc.
Han(1992) proposed a general Jacobian structure for complex column to carry out the parametric continuation on any parameter with the same Jacobian structure used by homotopy continuation.
EXAMPLES OF THE ALGORITHM Starting Point Problems and Mapping Function Consider the following nonlinear equations known as Himmelblau's function (Lin et ai, ,1987). CAn 11 S""pl-ff
5482
European Symposium on Computer Aided Process Engineering-2
4,-----------$
4,---------,--
--,
--.
Solution
All + Direction
·4.----------------1 ·2 2
.4.-1:"4---------11--------1
o
Xl
Fig. 2. Role of mapping function 2X 13 +2X1XZ - 21XI + Xzz - 7 XI Z + 2X1XZ + 2XZ3
4
Xl
-
Fig. 3. Path tracked in both directions
=
13Xz - 11
0
(10)
=0
(11)
If we solve the above equations by using the proposed Newton homotopy, a total of nine roots are obtained as shown in Fig. 1. It shows the different paths with different starting points. We can obtain only five roots when we start with (3.5, 2.4). As far as we know, there is no general starting point selection criteria for multi-dimensional problems. It is thought that the best way to find many solutions as possible is to try many different starting points.
Most problems in chemical engineering have constraints. To overcome these problems, a mapping function may be introduced. Lin (Lin et al . , 1987) proposed a mapping function having a type of yZ = x where x z O. Figure 2 shows the results of Himmelb1au's function with constraints, xl ~ 0 and x2 ~ 0, when Lin's mapping function was applied. As shown, a real positive solution on the opposite side of the starting point was found by tracking through the negative region. However, when we tracked both the positive and negative direction of t without the mapping, two other solutions that Lin's technique couldn't detect were found(Fig.3) . There may be a starting point that leads to a path connecting all solutions without violating constraints . A path connecting all four positive real roots were actually found with starting points of either (3.5, 1.3), (3.5, 1.0) or (3.5, 0.7). From this point, the parametric continuation whose path travels feasible region was proposed to implement steady state multiplicities for distillation problems. Performance Test Table 1 shows performance compansions of the proposed algorithm with other workers'(Rion and Van Brunt, 1990) for Himmelblau's function and Ho1odniok's one. As shown in Table 1, the number of steps and corre ctions were slightly larger but the number of Jacobian and function evaluations were much smaller than other's. Considering that the Jacobian evaluation is quite time-consuming (about more than 70% of the total computational load), the proposed algorithm would be more efficient for large nonlinear systems such as complex distillation columns. The Jacobian was evaluated by the finite difference method and a maximum of six iterations were allowed for each continuation step . The maximum allowable turning angle was 300.ideal turning angle was 300, and e dot was 0.5.
5483
European Symposium on Computer Aided Process Engineering-2
Table 1. Performance comparisons of the proposed algorithm
Holodniok's function
Himmelblau's function
NSTEP NCORR NFAIL NFE NJE
Frantz' Brunt
Rion" Brunt
Proposed MNR u NR
40 127 50 961 222
34 103 19 1147 106
45 100 6 612 152
47 139 10 524 103
Franz' Brunt
99 275 137 2994 603
Rion' Brunt
89 234 108 2961 517
Proposed MNRu NR
133 279 65 1774 431
140 316 64 1633 370
NSTEP=no. of continuation steps; NCORR=no. of corrections; NFAIL=no. of failures ; NFE=no. of function evaluations; NJE=no. of Jacobian evaluations * Data from Reference of Rion et al., 1990 ** MNR denotes modified Newton method.
Complex Distillation Column Problems An equation oriented column simulator has developed based upon the proposed homotopy and parametric continuation algorithm. To evaluate the performance of the developed system, more than 20 test examples were solved. All examples were extracted from literatures and most of them are difficut to solve having multiple steady state solutions. Thirteen types of column configuration were selected and three of them were heat-integrated column types. Due to space problems, only one experiment is summarized in this paper. The column configuration and basic operation data are shown in Fig. 4. The parametric continuation on distillate rate and reflux rate was studied. As shown in Fig. 5, under the composition specification of each component, the operating range of the distillate rate is very small(approximately from 19 to 21) and multiple solutions are found in this range. It can be seen that it is impossible to attain the design specification below the reflux rate of 180 Kgmollh. The same answer was obtained from Fig. 6 and the connection stream V and L have two solutions for a specified reflux rate. CONCLUSIONS An improved Newton homotopy continuation algorithm based on determinant monitoring and flexible convergence criteria is presented . There is a critical starting point above and below which the homotopy path shows completely different behavior. For a constrained system, such as most of chemical engineering problems, the parametric continuation will be robust without any violations of constraints. Also, this parametric continuation method is quite useful in the area of design and operability studies. Multiple solutions were found for most of the linked columns, especially those having concentration sepcification. It is noted here that not all variables have multiple solutions but some of variables may have multiple solutions.
S484
European Symposium on Computer Aided Process Engineering-2
Configuration N=20,1O NRl=5 Total condenser, P= 1.013 bar Feed Conditions Stage Temp. Rate(Kmollh) 100 10 383.8K Composition Benzene 0.2 Toluene 0.4 Xylene 0.4
F
""""'"I/Vl'-p
Specifications DBen
I---:/l/~-- B
Pre i BXyl
= 0.95
= 0.95 = 0.95
Fig. 4. Example: column with sidestripper
230,----------------. 180
I'l
130
J 80
lSOI1:l:9-------.-------~
20
Distillale(Kgmollh)
Fig. 5. Multiple solutions for distillate
21
I-,---vl
3
////
..-................•... ~
(~....
_
_-_ •..••.......•.............••...
200 250 Reflux Rete(Kgmollh)
300
Fig. 6. Multiple solutions for linking flowrate
REFERENCES Allgower, E. and Georg K.(1980). Simplical and continuation methods for approximating fixed points and solutions to systems of equations. SIAH Review,22,28 Chavez, C.R, J.D. Seader, and T.L. Wayburn (1986). Multiple steady state solutions for interlinked separation systems. Ind. Eng. Chem. FUlldam. ,25,566 Choi, S.H.(1990). Ph.D. Dissertation, Univ. of Missourri-Rolla Han, K.T. (1992). Ph.D. Dissertation, Seoul National University Lin, W.J.,J.D. Seader and T.L. Wayburn (1987).Computing multiple solutions to systems of interlinked separation column. AIChE J. 33,6,886 Rion W.L. and V.V. Brunt (1990). Differential geometry based homotopy continuation. COli/put. Chem. Eng. ,14,889
ABSTRACT This paper presents an efficient homotopy continuation method to solve a nonlinear eqaution system which describes chemical engineering problems . The characteristics of the Newton homotopy function is studied and parametric continuation is proposed to study steady state multiplicities for a constrained system. An equation based distillation simulator employing the proposed algorithm was implemented. The proposed algorithm and simulator were tested on many numerical problems and complex distillation columns. KEYWORDS Homotopy continuation, parametric continuation, solution, bifurcation
complex distillation column,
multiple
INTRODUCTION Recently, homotopy continuation methods have been used successfully to solve nonlinear equations due to their global convergent characteristics and their ability to find multiple solutions from one starting point. This paper shows how to improve the computational efficiency and capability of tracking homotopy path. This was achieved mainly by reducing the number of function evaluations and by avoiding segment jumping. The proposed algorithm had been applied to many complex distillation columns. A few workers have used the homotopy continuation method to get multiple steady state solutions for linked distillation columns(Chavez et al. ,1986, Lin et al., 1987). Parametric continuation followed by homotopy continuation was used to investigate the steady state multiplicities of the complex column and it was confirmed that this parametric continuation is quite useful in study design and operability limitations of the system. Features of the algorithm include: 1. The Euler prediction and Newton correction to track the path , 2. Step size control using determinant monitoring and curvature, 3. Flexible convergence criteria , 4. Capability of parametric continuation followed by homotopy continuation . ALGORITHM Newton Homotopy Function 5479
European Symposium on Computer Aided Process Engineering-2
S480
For solving the nonlinear equation system, f(x)=0, the Newton homotopy function h(x,t) is defined as (1)
h(x,t) = f(x) - (I-t)f(xo),
where x ERn, t e [0, I] and x = xo at t = O. If we differentiate h(x,t) by arclength, s, and introduce the definition of the arclength, we get the following IVP(initial value problem),
[
2.!!... (ws) aWT dw (Sk) ds
]
dw (Sk) ds
= e n+ 1
•
(2)
where w = [x,t], k = continuation step number and en+ 1 is a column vector of which the n+ It h column equals one and the others are zero. Any types of predictor-corrector schemes can be used to solve the IVP. simple Euler predictor and Newton corrector in this work.
We use the
Step Size Control In this study, Allgower's algorithm (1980) based on curvature and error contraction and Choi's (1990) determinant monitoring scheme were used. The determinant of the augmented matrix approachs zero as the solution does to another path(Allgower et a l. , 1980). The negative determinant value means that the path has already passed the bifurcation point and is tracking another path. At first, (k+ l)-th prediction is carried out with a step size calculated as l\sk+1 = 9id/K, where £lid is the ideal turning angle and K is the approximated curvature. Then, determinant is checked at this predicted point by equation (3). If this determinant doesn't satisfy equation (3), the proposed step size is reduced by a ratio of 0.5. This step size control strategy is essential in tracking the path by preventing segment jumping.
I detAk+I - detAk I Max ( IdetAlk+1 , IdetAlk ) < £det where A is Jacobian of the augmented matrix and
(3)
e det is tolerance for determinant
moni toring.
Convergence Criteria Each corrected point on the path is not affected by the error of the previous step (Allgower et al.,1980). To improve computational efficiency, we can use larger convergence tolerance during the intermediate path tracking than near target point(Le. t= I) or turning point. The following criteria provide flexible convergence tolerance and so increase computational efficiency.
II h(Wi+lk+ 1 ) II < Min(£mu,£min)
(4)
£mu = £tol/(det A(wo k+ I»lIn
(5)
If ( I det A(wO k+ l) I
(6)
Else
£min
=
< I ) £min = £cor I det A(wo k+ l) I
£cor I det A(wo k+ l)l/n I
European Symposium on Computer Aided Process Engineering-2
S481
5r-;;:=================::::::;'1 1- (3.5.3.0) - (3.5.2.7) _._- (3.5.2 .4) Solution I 0
3
·1
·3
·3
3
·1
5
X1
Fig. 1. Variation of homotopy path with different starting points
Parametric Continuation To solve the nonlinear equation expressed by f(x,p) = 0, we first set the parameter in order to satisfy degree of freedom. We can get a solution curve by changing the parameter p. This method is refered as parametric continuation and the difference from the homotopy continuation is that all points on the path tracked are solutions of the nonlinear system. To start the parametric continuation we need a solution at an initial point. The starting point for parametric continuation can be easily obtained by homotopy continuation . In order to perform parametric continuation on a selected parameter, pe, following Newton homotopy function defined in Rn~2 space can be obtained as, h(x,pc ,t)
=
[f(X,ps) ] pc - C
O
- (l - t) [f(X ,ps) ] pcD - C
=
0
(7)
where ps is a known parameter and c is a starting value of parameter pe. We can find a solution at t= I by applying homotopy continuation to the above equations . Parametric continuation on pc is proceeded by equation (8). It is regarded as homotopy continuation where pc is pameterized by t as given in equation (9). h(x,pc,t) t
=
where
=
[f(X,ps) ] pc - Ct
- (l - t) [
0
C - Ct
]
(c - pc)/(c - ce) , Ct
=
0
(8) (9)
is the given target value of parameter pc.
Han(1992) proposed a general Jacobian structure for complex column to carry out the parametric continuation on any parameter with the same Jacobian structure used by homotopy continuation.
EXAMPLES OF THE ALGORITHM Starting Point Problems and Mapping Function Consider the following nonlinear equations known as Himmelblau's function (Lin et ai, ,1987). CAn 11 S""pl-ff
5482
European Symposium on Computer Aided Process Engineering-2
4,-----------$
4,---------,--
--,
--.
Solution
All + Direction
·4.----------------1 ·2 2
.4.-1:"4---------11--------1
o
Xl
Fig. 2. Role of mapping function 2X 13 +2X1XZ - 21XI + Xzz - 7 XI Z + 2X1XZ + 2XZ3
4
Xl
-
Fig. 3. Path tracked in both directions
=
13Xz - 11
0
(10)
=0
(11)
If we solve the above equations by using the proposed Newton homotopy, a total of nine roots are obtained as shown in Fig. 1. It shows the different paths with different starting points. We can obtain only five roots when we start with (3.5, 2.4). As far as we know, there is no general starting point selection criteria for multi-dimensional problems. It is thought that the best way to find many solutions as possible is to try many different starting points.
Most problems in chemical engineering have constraints. To overcome these problems, a mapping function may be introduced. Lin (Lin et al . , 1987) proposed a mapping function having a type of yZ = x where x z O. Figure 2 shows the results of Himmelb1au's function with constraints, xl ~ 0 and x2 ~ 0, when Lin's mapping function was applied. As shown, a real positive solution on the opposite side of the starting point was found by tracking through the negative region. However, when we tracked both the positive and negative direction of t without the mapping, two other solutions that Lin's technique couldn't detect were found(Fig.3) . There may be a starting point that leads to a path connecting all solutions without violating constraints . A path connecting all four positive real roots were actually found with starting points of either (3.5, 1.3), (3.5, 1.0) or (3.5, 0.7). From this point, the parametric continuation whose path travels feasible region was proposed to implement steady state multiplicities for distillation problems. Performance Test Table 1 shows performance compansions of the proposed algorithm with other workers'(Rion and Van Brunt, 1990) for Himmelblau's function and Ho1odniok's one. As shown in Table 1, the number of steps and corre ctions were slightly larger but the number of Jacobian and function evaluations were much smaller than other's. Considering that the Jacobian evaluation is quite time-consuming (about more than 70% of the total computational load), the proposed algorithm would be more efficient for large nonlinear systems such as complex distillation columns. The Jacobian was evaluated by the finite difference method and a maximum of six iterations were allowed for each continuation step . The maximum allowable turning angle was 300.ideal turning angle was 300, and e dot was 0.5.
5483
European Symposium on Computer Aided Process Engineering-2
Table 1. Performance comparisons of the proposed algorithm
Holodniok's function
Himmelblau's function
NSTEP NCORR NFAIL NFE NJE
Frantz' Brunt
Rion" Brunt
Proposed MNR u NR
40 127 50 961 222
34 103 19 1147 106
45 100 6 612 152
47 139 10 524 103
Franz' Brunt
99 275 137 2994 603
Rion' Brunt
89 234 108 2961 517
Proposed MNRu NR
133 279 65 1774 431
140 316 64 1633 370
NSTEP=no. of continuation steps; NCORR=no. of corrections; NFAIL=no. of failures ; NFE=no. of function evaluations; NJE=no. of Jacobian evaluations * Data from Reference of Rion et al., 1990 ** MNR denotes modified Newton method.
Complex Distillation Column Problems An equation oriented column simulator has developed based upon the proposed homotopy and parametric continuation algorithm. To evaluate the performance of the developed system, more than 20 test examples were solved. All examples were extracted from literatures and most of them are difficut to solve having multiple steady state solutions. Thirteen types of column configuration were selected and three of them were heat-integrated column types. Due to space problems, only one experiment is summarized in this paper. The column configuration and basic operation data are shown in Fig. 4. The parametric continuation on distillate rate and reflux rate was studied. As shown in Fig. 5, under the composition specification of each component, the operating range of the distillate rate is very small(approximately from 19 to 21) and multiple solutions are found in this range. It can be seen that it is impossible to attain the design specification below the reflux rate of 180 Kgmollh. The same answer was obtained from Fig. 6 and the connection stream V and L have two solutions for a specified reflux rate. CONCLUSIONS An improved Newton homotopy continuation algorithm based on determinant monitoring and flexible convergence criteria is presented . There is a critical starting point above and below which the homotopy path shows completely different behavior. For a constrained system, such as most of chemical engineering problems, the parametric continuation will be robust without any violations of constraints. Also, this parametric continuation method is quite useful in the area of design and operability studies. Multiple solutions were found for most of the linked columns, especially those having concentration sepcification. It is noted here that not all variables have multiple solutions but some of variables may have multiple solutions.
S484
European Symposium on Computer Aided Process Engineering-2
Configuration N=20,1O NRl=5 Total condenser, P= 1.013 bar Feed Conditions Stage Temp. Rate(Kmollh) 100 10 383.8K Composition Benzene 0.2 Toluene 0.4 Xylene 0.4
F
""""'"I/Vl'-p
Specifications DBen
I---:/l/~-- B
Pre i BXyl
= 0.95
= 0.95 = 0.95
Fig. 4. Example: column with sidestripper
230,----------------. 180
I'l
130
J 80
lSOI1:l:9-------.-------~
20
Distillale(Kgmollh)
Fig. 5. Multiple solutions for distillate
21
I-,---vl
3
////
..-................•... ~
(~....
_
_-_ •..••.......•.............••...
200 250 Reflux Rete(Kgmollh)
300
Fig. 6. Multiple solutions for linking flowrate
REFERENCES Allgower, E. and Georg K.(1980). Simplical and continuation methods for approximating fixed points and solutions to systems of equations. SIAH Review,22,28 Chavez, C.R, J.D. Seader, and T.L. Wayburn (1986). Multiple steady state solutions for interlinked separation systems. Ind. Eng. Chem. FUlldam. ,25,566 Choi, S.H.(1990). Ph.D. Dissertation, Univ. of Missourri-Rolla Han, K.T. (1992). Ph.D. Dissertation, Seoul National University Lin, W.J.,J.D. Seader and T.L. Wayburn (1987).Computing multiple solutions to systems of interlinked separation column. AIChE J. 33,6,886 Rion W.L. and V.V. Brunt (1990). Differential geometry based homotopy continuation. COli/put. Chem. Eng. ,14,889