2320
Ind. Eng. Chem. Res. 2008, 47, 2320-2330
PROCESS DESIGN AND CONTROL Nested Tabu Search (TS) and Sequential Quadratic Programming (SQP) Method, Combined with Adaptive Model Reformulation for Heat Exchanger Network Synthesis (HENS) Xi Chen,* Zhaohua Li, Jing Yang, and Zhijiang Shao State Key Lab of Industrial Control Technology, Department of Control Science and Engineering, Zhejiang UniVersity, Hangzhou 310027, People’s Republic of China
Lingyu Zhu College of Chemical Engineering and Materials Science, Zhejiang UniVersity of Technology, Hangzhou 310014, People’s Republic of China
A nested method combining tabu search (TS) and sequential quadratic programming (SQP) for the mixedinteger nonlinear programming (MINLP) is proposed. It involves MINLP in a nested inner and outer loop, where TS in the outer loop involves the integer variables, and SQP in the inner loop solves the nonlinear programming (NLP) with integers specified from the outer loop. Good performance has been demonstrated by several mathematic examples and a pump system example, but problems with regard to heat exchanger network synthesis (HENS) are still encountered. Further analysis shows that even given the optimal structure of a heat exchanger network (HEN), it is difficult to solve the inner loop NLP. Based on the Yee and Grossmann model, an adaptive model reformulation method is later proposed, to bridge the inner loop and the outer loop. Given the integer variables by the outer loop, it automatically simplifies the NLP model, in regard to both dimension and complexity, before delivering it to the inner loop SQP method. The results are presented and discussed through examples. I. Introduction Mixed-integer nonlinear programming (MINLP) involves both continuous and integer variables; many chemical process design and operational problems are MINLP problems. Until now, several techniques have been developed for convex MINLP problems, such as the Branch and Bound algorithm (BB),1 the Generalized Benders Decomposition algorithm (GBD),2 and the Outer Approximation algorithm (OA).3 However, global optimality cannot be guaranteed for many nonconvex cases. Thus, heuristic search methods, such as simulated annealing (SA),4 genetic algorithm (GA),5 and tabu search (TS),6,7 are also broadly applied to solve the nonconvex MINLP problems, because of their good performance in regard to global optimality. Heat exchanger network synthesis (HENS) has received significant attention during the past 30 years for the large savings that are achieved in energy costs. In the mathematic programming framework, HENS is generally formulated as MINLP superstructure models.8-10 Previous studies have revealed that solving HENS either sequentially or with a stagewisesimultaneous formulation is NP-hard.11 The deterministic methods such as GBD and OA have been successfully applied to HENS. However, their computation time changes exponentially with the problem size. Sometimes, the global optimality cannot be guaranteed. A few studies have reported the applications of
heuristic search methods to HENS, such as SA,12 GA,13,14 and TS.15 Although the heuristic search methods can find global optimal solutions in many cases, they are normally slower than the deterministic methods. Besides, as most heuristic search methods such as TS are originally developed for combinatorial optimization, they may not perform at their best in regard to dealing with continuous variable optimizations, especially those with constraints. Considering the pros and cons of heuristic and deterministic methods, a nested method that combines TS and sequential quadratic programming (SQP) is proposed in this paper. It involves the MINLPs in nested inner and outer loops, whereas TS, which is a typical heuristic search method, involves integer variables in the outer loop, and SQP, which is a deterministic method, solves the sub-NLP with integers that are specified from the outer loop in the inner loop. When it is applied to HENS, an adaptive model reformulation technique is also proposed, to greatly simplify the solution of the inner loop NLP. The structure of this paper is organized as follows. In next section, the idea of combining the TS algorithm with SQP is described. The performance is shown by several case studies; and the unsuccessful applications to HENS also are presented. In section 3, the mechanism of adaptive model reformulation for HENS is described. The results and discussions on HENS are given in section 4. The last section concludes the paper. II. Nested TS & SQP Approach for MINLP
* To whom correspondence should be addressed. Tel.: 86-57187951970, Ext. 401. Fax: 86-571-87951068. E-mail address: xichen@ iipc.zju.edu.cn.
1. Structure of Nested TS & SQP. Generally, the MINLPs can be defined as follows:
10.1021/ie071245o CCC: $40.75 © 2008 American Chemical Society Published on Web 03/04/2008
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008 2321
min f(x,y) x,y
(1)
s.t. h(x,y) ) 0 g(x,y) e 0 x ∈ X ⊂ Rn y ∈ Y integer where x is a vector of n continuous variables, y is a vector of m integer variables (normally binary variables), h denotes equality constraints, g denotes inequality constraints, X and Y are sets, and f is the objective function. Cases with integer variables other than binary ones can also be transformed to binary forms.16 Moreover, if some of the constraints relate only to integer variables, the formulation can be represented as
min f(x,y) x,y
(2)
s.t. h(x,y) ) 0 g(x,y) e 0 p(y) ) 0 q(y) e 0 x ∈ X ⊂ Rn y ∈ {0,1}m The optimal solution to the aforementioned problem is also the same as the following form:
min min f(x,y)
y∈{0,1}m x⊂Rn
(3)
which is a nested optimization with two loops. The outer loop is an integer optimization that involves binary variables only as
min V(y) y
(4)
s.t. p(y) ) 0 q(y) e 0 y ∈ {0,1}m The inner loop is an NLP with y′ specified by the outer loop as follows:
V(y′) ) min f(x,y′) x
(5)
s.t. h(x,y′) ) 0 g(x,y′) e 0 x ∈ X ⊂ Rn In this paper, we propose a combined TS and SQP method to examine MINLPs in a nested inner loop and outer loop. Figure 1 shows a schematic description of this method. In the outer loop, TS is used to address the integer optimization. It generates the feasible binary variables y′ for the inner NLP. SQP, which is an efficient algorithm for NLP, addresses the subsequent NLP problems in the inner loop. The optimal solution to the inner NLP then is returned to the outer loop TS, to generate a new solution for the next iteration. This procedure continues until the termination criterion of the outer loop is satisfied.
Figure 1. Structure of the nested TS & SQP approach.
2. Tabu Search Algorithm. Tabu search (TS) is a metaheuristic method originally developed for combinatorial optimization by Glover.6,7 It uses adaptive memory and responsive exploration to search the solution space effectively and prevent it from being trapped in a local optimum. Good performances have been demonstrated by many applications, such as scheduling problems,17 traveling salesman problems,18 and HENS problems.15 The algorithm used in this paper has been developed by the authors, and the details are presented as follows. 2.1. Neighbor Generation. In the outer loop, an integer optimization is conducted by TS. First, at each iteration, a group of initial neighbor solutions is generated; the best among the neighbors is used to update the current solution. A set Y(k) ini is defined as follows:
Y(k) ini ) {y′
|
m
|y′i - yi| ) k} ∑ i)1
(6)
where m is the dimension of the binary variables. It means that the elements of Y(k) ini are exactly k bits different from the current solution, y. The initial neighbor set (Yini) then can be constructed as a union of sets as follows: [m/2]
Yini ) ∪ Y(k) m
(7)
k)1
where Y(k) m is constructed by randomly selecting m points from each set of Y(k) ini , for k values from 1 to the nearest integer of m/2 toward infinity. Note that the merit of this type of definition is that we can limit the number of neighbors with an upper bound and keep it evenly distributed. After the initial neighbors are generated, these points are checked by the pure integer constraints; the violated ones are removed from Yini. If none of the solutions in Yini meets the pure integer constraints, then a new Yini is randomly generated, according to the previously described rules. 2.2. Tabu and Aspiration Criteria. In this paper, 2-fold tabu lists are constructed, to improve the efficiency of TS. Tabu List One (TL1) saves the information of the points that have been visited in the past L1 iterations. If there is any new point y′ satisfying the conditions m
∃ y ∈ TL1 and
|y′ - y | e 1 ∑ j)1 j
j
(8)
it is considered to be a taboo point. This tabu area covers the points recently visited and those with only one bit different from any solution in TL1. The length of TL1 is set as L1. Tabu List Two (TL2) saves the information of local optimal points. If the objective function has not been improved for k iterations, the corresponding point is treated as a local optimum. The length
2322
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008
of TL2 is set as L2. If there is any point y′ satisfying the conditions m
∃ y ∈ TL2 and
|y′j - yj| e 2 ∑ j)1
(9)
it is also tabooed. This taboos the points with no more than two bits different from any solution in TL2. Through this way, it reduces the probability of avoiding being trapped into a local optimum. The algorithm also gives a probability p to aspire a tabu point as a neighbor candidate of the current point. A random value is generated before the tabu checking of each initial neighbor. If the value is less than a preset parameter p, this point is set as a neighbor, regardless of whether it has or has not been tabooed previously. The value of parameter p determines the degree of aspiration. It can be adjusted according to the size of the MINLP. 2.3. Termination Criterion. The procedure is designed to terminate if the maximum number of iterations (Imax) is reached or if the objective function has not been improved for a given number of iterations (M). The maximum number of iterations should increase with the dimension of the binary variables. In this paper, the default maximum number of iterations (Imax) is defined as
Imax ) round
1 )2 ] [(20m m
(10)
which rounds the elements up to the nearest integer. The range for M in this paper is between 0.2(Imax) and 0.5(Imax). 3. SQP Algorithm.19 With integers specified by the outer loop, the inner loop becomes a typical constrained NLP. SQP has been demonstrated with good performance, in regard to constrained optimization problems. Therefore, it is selected in the inner loop of our nested method. SQP is a class of optimization methods through solving a quadratic programming (QP) sub-problem at each iteration. Each QP sub-problem minimizes a quadratic model of a certain modified Lagrangian function subject to linear constraints. The solution to the QP sub-problem is used to form a search direction for a line search procedure. The merit function is reduced along each search direction, to ensure convergence from any starting point. Given the NLP problems described by expression 5, the QP subproblem can be defined by linearizing both the inequality and equality constraints as follows:
1 min dTHkd + ∇f(xk)Td d∈Rn 2
(11)
s.t. ∇h(xk)Td + h(xk) ) 0 ∇g(xk)Td + g(xk) e 0 The sub-problem is solved by the active-set methods. The solution is used to form a new iterate:
xk+1 ) xk + Rkdk
(12)
where the step length parameter (Rk) is determined by an appropriate line search procedure, so that a sufficient decrease in a merit function is obtained. The matrix, Hk, is a positive definite approximation of the Hessian matrix of the Lagrangian function; it can be updated by any of the quasi-Newton methods. The iteration continues until it converges.
4. Basic Procedure of Nested TS & SQP Algorithm. As shown in Figure 2, the basic procedure of the algorithm is given as follows: Step 1: Initialize the control parameters, such as the maximum number of iterations (Imax), the aspiration probability (p), and the length of the TL1 and TL2, etc. Then randomly select a point as the current solution with empty initial tabu lists. Step 2: Construct an initial neighbor solution set Yini of the current solution y and select the first element of Yini as the candidate. Step 3: Judge if the candidate satisfies the pure integer constraints; if it is satisfied, go to Step 4; otherwise, go to Step 6. Step 4: Evaluate the aspiration criteria on the selected candidate; if it is satisfied, add the candidate into the new neighbor set Ynew; otherwise, go to Step 5. Step 5: Judge the tabu criteria; if the candidate is tabooed, go to Step 6; otherwise, add the candidate into Ynew. Step 6: Judge if all the points in the neighbor solution set Yini have been evaluated. If yes, go to Step 7; otherwise, select the next point in Yini as the candidate and return to Step 3. Step 7: Judge if Ynew is empty; if it is empty, return to Step 2; otherwise, go to Step 8. Step 8: For each element in Ynew, apply SQP to solve the inner NLP. Step 9: Select the best solution in Ynew as the current solution and update the tabu lists. Step 10: Judge the termination criterion; if yes, stop the program and output the best-known solution; otherwise, return to Step 2 and repeat the procedure. 5. Case Studies. The tests in this section, as presented in Appendix A, include all the MINLP examples from the handbook.20 Simulations are conducted in MATLAB 7.1 with a Pentium IV CPU (1.7 GHz and 512M RAM). The inner loop SQP algorithm uses the standard function of the MATLAB optimization toolbox with default parameter settings. The other programs, including the TS algorithm, are coded by the authors. Each test is ran 20 times from different random initials. The parameter settings and results are summarized in Tables 1 and 2, respectively. The slight difference in the parameter settings is due to the scale difference among the cases. In Table 2, the second column lists the successful run numbers to reach the global optimal solution among a total of 20 runs. We can see that this method effectively avoids the possibility of being trapped in local minima. The first five cases (100%) reach the global optimal solutions. Even for case 6, 18 of the 20 runs successful converge to the global optimal solution. The third column lists the total time cost for the 20 runs. Note that cases 5 and 6 are not binary variable cases; we first convert them to binary variables, following the method presented by Floudas.16 After this transformation, there are 15 binary variables in case 6. The increase in dimension also increases the complexity and the time cost. The successful rate of case 6 can be further improved by enlarging the parameters in Table 1; this would consequently cost more time for the TS algorithm to locate the global optimal solution as a tradeoff. In the last column of Table 2, “AvgTabuProb” means the average taboo probability in the procedure of neighbor generation, which is defined as the ratio of the number of the tabooed neighbor candidates over the number of all initial candidates. The average taboo probability reveals the efficiency of the TS algorithm: if it is too small, it means the taboo action is at low effect; if it is too large, it is difficult for the program to jump out of local optimums. The average tabu probability of each case is ∼10%-25%, which indicates that the TS mech-
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008 2323
Figure 2. Flowchart of the nested TS & SQP algorithm. Table 1. Parameter Settings for Cases 1-6 case
Imax
M
k
p
L1
L2
1 2 3 4 5 6
5 2 5 30 5 100
3 1 3 15 3 20
5 1 3 5 3 8
0.08 0.08 0.08 0.08 0.08 0.08
3 1 3 10 3 30
3 1 3 5 3 10
Table 2. Optimization Results of Cases 1-6 case
Global
Time (s)
AvgTabuProb
1 2 3 4 5 6
20 20 20 20 20 18
8.34 0.81 15.4 35.3 2.82 256.3
0.226 0.2 0.149 0.1545 0.13 0.1692
anism used in the outer loop prevents the solution from revisiting the points that previously have been visited effectively. 6. HENS Cases. In this section, the nested TS & SQP method is applied to three HENS cases. They are all based on the model proposed by Yee and Grossmann.8,9 The model and case descriptions are presented in Appendix B. In our method, the HEN structure is represented by a vector Y ) [Zijk,Zcu,i,Zhu,j]. For each case, the program is ran 20 times from different random initials. However, none of them reaches the global optimal solutions. Further study shows that this failure is because of the fact that the inner loop SQP cannot locate the
optimal solution. Even given the optimal HEN structure, the inner NLP cannot be successfully solved. SQP first converts the original NLP problem to a QP problem; the QP problem then is solved using the active-set method to find a decreasing search direction, dk, to determine the next solution, as shown in eq 12. For HEN cases, however, the active-set method cannot find a feasible starting point, because of the ill condition of the original NLP problem. Under this situation, the direction dk generated by the sub-QP cannot lead the main SQP to the global optimal solution. Note that the Yee and Grossmann model was constructed as a general superstructure model, which introduced some auxiliary decisive variables and constraints. With specified binary variables, the inner NLP is ill-conditioned and, thus, difficult to solve. The constraints of the original problem are also dependent. For example, equalities B.4 and B.7 can lead to equality B.3 (see Appendix B). To enable the application of our nested method to HENS, we further propose a mechanism, called the adaptive model reformulation method, between the inner loop and the outer loop to adaptively reformulate the inner loop NLP model after the binary variables have been specified by the outer loop TS method. This model reformulation is anticipated to facilitate the solution of the inner NLP for HENS. III. Adaptive Modeling for HEN Synthesis Given the binary variables by the outer-loop TS method, the HEN structure is also fixed. As shown in Appendix B, the inner-
2324
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008
Figure 3. Structure of adaptive inner NLP.
Figure 4. Static and dynamic components of the adaptive variables, bounds, and constraints.
loop optimization that is based on the Yee and Grossmann model becomes a typical NLP with the following linear constraints:
min f(X)
(13)
s.t. AX e B AeqX ) Beq Xlb e X e Xub Here, f is the objective function; A, B, Aeq, and Beq are the coefficients of the linear inequalities and equalities; Xlb and Xub are the respective lower and upper bounds of the continuous variables X. As the original MINLP model is generalized for a superstructure, it becomes dependent after fixing the binary variables by the outer loop. Therefore, we propose adding a block to adaptively reformulate a new inner NLP model between the inner loop and the outer loop. As shown in Figure 3, given the parameters Y, HP, CP, SI, and Data, the reformulation includes the generation of reduced dimension variables and their bounds (X*, X*lb, X*ub), the generation of constraints (A*, B*, A* eq, B* eq), and the rebuilding of the object function (f*). Details are shown as follows. 1. Generation of Reduced Dimension Variables and Their Bounds. The new set of the continuous variables are denoted by a vector X*, which stores information regarding the temperature and heat load. It can be divided into two parts: a static part and a dynamic part, as shown in Figure 4. The static part, which is denoted as Xs, includes temperatures of the hot and cold streams at each stage. These variables always exist, regardless of the outer loop match. The other continuous variables in the original superstructure model may have zero values or be physically meaningless, relative to variation of the HEN structure. For example, Qijk only exists when Zijk ) 1; this would also hold for Qcu,i and Qhu,j. In addition, ∆Tijk denotes the temperature difference between Tik and Tjk at stage k only if Zijk ) 1; it is physically meaningless if Zijk ) 0. Therefore, the original ∆Tijk are excluded in the new model, because it can be substituted by ∆Tijk ) Tik - Tjk when Zijk ) 1; this also holds for ∆Tcu,i and ∆Thu,j. We also introduce a vector Xd, which represents the dynamic variables used to store the information of the heat exchange of existing matches. It varies dimensionally with the outer loop solution. In accordance with the different
variable types, the bounds are also divided into the static and dynamic parts, as X*lb ) [Xlbs; Xlbd] and X*ub ) [Xubs; Xubd]. In this way, the dimension of the original model can be greatly reduced by removing unnecessary variables. Consequently, the time cost can be reduced as the scale of the inner NLP is reduced. In addition, it helps to eliminate the illness of the inner loop NLP after removing the zero-value variables and physically nonexistent variables. Details are shown in the next section. 2. Generation of Constraints. Given the HEN structure, the constraints are also adaptively rebuilt, because most of them involve integer variables. Following the original idea of the Yee and Grossmann model, the constraints can be categorized as follows. 2.1. Overall Heat Balance for Each Stream. Because each hot stream must exchange heat with a cold stream or cold utility, and each cold stream must exchange heat with a hot stream or hot utility, equalities B.3.1 and B.3.2 in Appendix B always activate and are included in the reformulated model. 2.2. Heat Balance at Each Stage. If ∑i∈HP Zijk g1, for j ∈ CP, k ∈ SI (i.e., a match between cold stream j and one hot stream at stage k exists), then equality B.4.1 in Appendix B activates; otherwise, this constraint is deactivated and excluded in the reformulated model. A similar analogy is observed for equality B.4.2. 2.3. Assignment of Superstructure Inlet Temperatures. Because equalities B.5.1 and B.5.2 in Appendix B always activate, regardless of the binary variables, they are always included in the reformulated model. 2.4. Feasibility of Temperatures. If ∑j∈CP Zijk g1, for i ∈ HP, k ∈ SI (i.e., the match between hot stream i and one cold stream at stage k exists), then inequality B.6.1 in Appendix B holds in its original form; otherwise, it is reformulated as an equality constraint, Ti,k ) Ti,k+1. If ∑i∈HP Zijk g1, for j ∈ CP, k ∈ SI (i.e., the match between cold stream j and one hot stream at stage k exists), then inequality B.6.2 in Appendix B holds in its original form; otherwise, it is reformulated as Tj,k ) Tj,k+1. If Zcu,i ) 1, for i ∈ HP (i.e., the match between hot stream i and cold utility exists), then inequality B.6.3 in Appendix B holds in its original form; otherwise, it is reformulated as Tout,i ) Ti,NS+1. If Zhu,j ) 1, for j ∈ CP (i.e., the match between cold stream j and hot utility exists), then inequality B.6.4 in Appendix B holds in its original form; otherwise, it can be reformulated as Tout,j ) Tj,1. 2.5. Hot and Cold Utility Loads. If Zhu,j ) 1, j ∈ CP (i.e., the match between cold stream j and hot utility exists), then equality B.7.1 in Appendix B holds in its original form; otherwise, it is simplified as Tout,i ) Ti,NS+1. If Zcu,i ) 1, i ∈ HP (i.e., the match between hot stream i and cold utility exists), then equality B.7.2 in Appendix B holds in its original form; otherwise, it is simplified as Tout,j ) Tj,1. 2.6. Logical Constraints. Similarly, inequalities B.8.1, B.8.2, and B.8.3 in Appendix B hold only when Zijk ) 1, Zhu,j ) 1,
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008 2325
and Zcu,i ) 1, respectively; otherwise, they are removed from the reformulated model. 2.7. Minimum Approach Temperature Constraints. If Zijk ) 1, ∆Tijk is regard as an auxiliary variable as: ∆Tijk ) Tik Tjk, and inequality B.9 in Appendix B is rebuilt as Tik - Tjk g EMAT; when Zijk ) 0, inequality B.9.1 in Appendix B is deactivated. A similar analogy is observed for inequalities B.9.2, B.9.3, and B.9.4 in Appendix B; inequality B.9.5 in Appendix B is excluded in the new model, because it has been substituted into the new constraints. As a summary, the constraints can also be divided into two types: the static one, which is independent of the HEN structure, and the dynamic one, which is dependent on the HEN structure. As shown in Figure 4, the constraints that do not change with the outer-loop structure are represented by the following coefficients: As, Aeqs, Bs, and Beqs. The dynamic constraints, on the other hand, are represented by the following coefficients: Ad, Aeqd, Bd, and Beqd. Following the steps previously described, the constraints can be automatically developed as follows:
A*X* e B*
(13a)
A*eqX* ) B*eq
(13b) Figure 5. Flowchart of the nested TS & SQP method, combined with adaptive model reformulation for HENS.
where
( ) ( )
X X* ) Xs d B B* ) Bs d
( ) ( ) As 0 0 Ad
A* )
B B*eq ) Beqs eqd A*eq )
(
Aeqs 0 Aeqd 0
)
3. Generation of the Objective Function. Similarly, the objective function can also be automatically rebuilt in accordance with the given HEN structure. The original investment cost terms CFijZijk and Cij(Qijk/U(LMTD)i,j)β exist only if Zijk ) 1; otherwise, they are ignored in the reformulated objective function. This also holds for cases with Zcu,i and Zhu,j. The operation cost terms (CcuQcu,i and ChuQhu,j) exist only if Zcu,i and Zhu,j, respectively, are equal to 1. 4. Overall Structure and Flowchart. As shown in Figures 3 and 4, the model can be adaptively reformulated in accordance with the HEN structure specified by the outer loop. The variables, bounds, and constraint coefficients are divided into two types: the static one, which is independent of integer variables, and the dynamic one, which is dependent on integer variables. Based on the mechanism of the adaptive modeling previously described, the HEN synthesis can be implemented
using the nested TS & SQP method with adaptive model reformulation, as shown in Figure 5. The TS algorithm involves integer variables in the outer loop. After initiation, the TS algorithm first generates feasible neighbors for the current HEN structure. Then, for each neighbor, adaptive model reformulation will automatically generate a new NLP, according to the HEN structure specified by the integers. SQP is later used to solve each sub-NLP. The optimal solution among those of all the subNLPs is selected to update the current HEN structure; tabu lists are also updated. The program repeats until the termination criterion in the outer loop is satisfied. In comparison to the original nested TS & SQP method for HENS, the major difference lies in the adaptive model reformulation, which is used to bridge the inner loop and the outer loop. The performance is presented and discussed in the next section. IV. Results and Discussion To illustrate how the adaptive model reformulation improves the inner loop optimization, we first assume that the HEN structures that are generated by the outer loop are already at their optimal structures. The tests are conducted to run the inner NLP 100 times from different random initials. The results are presented in Table 3, which compares the success times to locate the global optimal solutions (Global), the number of continuous variables (Nx), the number of inequality constraints (Nneq), the number of equality constraints (Neq), the condition number (CN), and the time cost (Time) for the inner NLP before and after the model reformulation. As described in section II.6, SQP cannot solve the inner loop NLP, even given the optimal
Table 3. Comparison between Original NLP and Reformulated NLP of Cases 7, 8, and 9 Case 7
Case 8
Case 9
NLP
original
reformulation
original
reformulation
original
reformulation
Global Nx Nneq Neq CN Time (s)
0 52 60 24 5.1 × 1016 30.063
100 21 25 20 358.459 7.743
0 72 84 30 1.5 × 1016 55.89
100 25 27 25 14.876 3.26
0 175 200 50 9.26 × 1016 282.34
100 40 48 40 126.518 8.08
2326
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008
Table 4. Iterative Process for Case 7 from a Fixed Start Point Iter 1 2 3 4 5 6 7 8
Y
NX
1111111111111111 1000010101001010 1000010100001010 1000100100001010 1010100100001010 1010100100001110 1010100100000110 1010100100000100
Nneq
32 22 21 21 22 23 22 21
...
Neq
Obj Fun
feasible?
44 26 24 24 26 27 26 25
20 20 20 20 20 20 20 20
1.0 × 1.0108 × 105 0.9505 × 105 0.9505 × 105 0.8726 × 105 0.9234 × 105 0.8096 × 105 0.7471 × 105
no yes yes yes yes yes yes yes
Nneq
Neq
Obj Fun
feasible?
20 20 20 20 20 20 20
1.0 × 1.0 × 1010 1.0 × 1010 0.9357 × 105 0.8721 × 105 0.8096 × 105 0.7471 × 105
no no no yes yes yes yes
1010
...
Table 5. Iterative Process for Case 7 from Another Random Start Point Iter 1 2 3 4 5 6 7
Y
NX
1111010111001110 1000110001000010 1000100001000110 1010100101001010 1010100101000110 1010100101000100 1010100100000100
27 21 21 23 23 22 21
...
35 25 24 28 28 27 25
1010
...
Table 6. Iterative Process for Case 8 from a Fixed Start Point Iter 1 2 3 4 5 6 7 8 ...
Y 11111111111111111111111 01000000011000100001001 01000000011000100011001 01000000011000100011101 00000000010000100011101 00000000000000100011101 00000100000000110011001 00000000000000110011001
NX
Nneq
43 26 27 28 26 25 26 25
61 30 31 32 28 26 29 27
Neq
Obj Fun
feasible?
25 25 25 25 25 25 25 25
1.0 × 1.0 × 1010 1.0 × 1010 1.032 × 105 0.9408 × 105 0.8419 × 105 0.8592 × 105 0.8204 × 105
no no no yes yes yes yes yes
1010
...
structure. As a comparison, the NLP after the model reformulation can reach the optimal solution at a 100% success rate. Taking case 7, for example, the number of continuous variables is dramatically decreased from 52 to 21; the numbers of inequality and equality constraints are also decreased, from 60 to 25 and from 24 to 20, respectively. The condition number, as a consequence, is obviously improved from 5.1 × 1016 to 358.459 after the model reformulation, which means the problem is efficiently improved from an ill-posed one to a well-posed one. This is very suitable for SQP to solve it. The total time cost of the 100 times of the inner NLP is 7.743 s. The effect can be obvious for larger-scale problems. As shown by the results of case 9 in Table 3, although the scale is much larger than those of the other two cases, the time cost of the inner NLP has the same order of magnitude after the model reformulation. As a conclusion, we can see that the adaptive model reformulation greatly reduces the inner-loop optimization scale and improves its solvability. We also test how the nested TS & SQP method iteratively leads the HENS to the optimal solution. Initializing all integer variables as 1, the process iterates and soon converges to the optimal structure. As shown in Table 4, after only eight iterations in the outer loop, case 7 reaches its best structure, as presented
by the second column in the table. The third, fourth, and fifth columns in the table respectively list the numbers of the variables, the inequalities, and the equalities, after the model reformulation under each structure. Among them, the number of equalities is independent of the HEN structure, whereas the other two parameters can be greatly reduced as the structure solution approaches the optimal one. Moreover, note that the first solution is infeasible. It does not satisfy the condition of no stream splitting (i.e., the pure integer constraint). To avoid having the iteration go back to those solutions, we assign a large penalty value, 1.0 × 1010, to them. After the eighth iteration, the problem cannot be improved and terminates after the next M iterations. The optimal objective function is $74 711. Figure 6 illustrates the optimal solution for case 7, which is same as the best-known solution. Later, we repeat the test for case 7 from another random initial, [1 1 1 1 0 1 0 1 1 1 0 0 1 1 1 0]. As shown in Table 5, the procedure reaches the same optimal solution at the seventh iteration. Table 6 lists the iterative results for case 8. Also starting from an initial with values that all are 1, it reaches the best-known solution at the eighth iteration. Figure 7 illustrates its optimal solution.
Figure 6. Optimal structure of HEN case 7.
Figure 7. Optimal structure of HEN case 8.
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008 2327 Table 7. Parameter Settings for Cases 7, 8, and 9 case
Imax
M
k
P
L1
L2
7 8 9
50 100 200
30 30 70
20 20 20
0.08 0.08 0.08
50 60 100
30 40 70
Table 8. Optimization Results of Cases 7, 8, and 9 case
NY
Global
average time (s)
7 8 9
16 23 60
27 26 17
63.8 129.7 3556.8
Because the TS algorithm is a heuristic search method, we further test these three cases for 30 runs from different random initials. The parameter settings of the outer loop TS are listed in Table 7; the parameters of the inner loop SQP are set at the default values of the Matlab function. The statistical results are shown in Table 8. The second column in Table 8 presents the dimension of the integer variables of each case. The third column in Table 8 lists the success times to reach the global optimal solution among the 30 runs. The last column in Table 8 shows the average time cost of each run. We can see that, for cases 7 and 8, most runs can reach the global optimal solutions very fast. Taking case 7, for example, the number of the binary variables is 16 (i.e., there are 216 possible HEN structures). The global optimal solution can be obtained only if the optimal HEN structure represented by the binary variables has been experienced in the outer loop. The probability of locating the global optimal solution can be further improved if we enlarge the number of neighbors, and the same can be observed for the parameters in Table 7. Of course, as a tradeoff, the time cost will also increase. For case 9, there are 60 binary variables (i.e., there are 260 possible HEN structures). Therefore, it is relatively difficult to obtain the global optimal structure, but still 17 of the 30 runs succeed with the best-known solution, as shown in Figure 8. Moreover, in the other 13 runs, the entire program stops at a local solution, the structure of which is illustrated in Figure 9. Its structure is very different from that of the best-known solution in Figure 8. However, the objective function of this sub-optimal solution is $43 331, which is very similar to the best-known one ($43 330). We have not found other literature that has reported on this sub-optimal solution. The slight difference in objective function makes the global optima very difficult to locate. On the other hand, the slight improvement has little significance in regard to engineering. We may also accept it as a global optimal solution.
Figure 9. Inferior optimal structure of HEN case 9.
mixed-integer optimization problems. The TS algorithm is used in the outer loop to address integer variables. The SQP algorithm is used in the inner loop to solve the sub-NLPs after the integer part is given by the TS algorithm. Based on the Yee and Grossmann model,8,9 an adaptive model reformulation method is further proposed to combine with the nested TS & SQP approach for HENS. Given integer variables by the outer loop, it automatically simplifies the NLP model, in regard to both dimension and complexity, before delivering it to the inner loop SQP method. The good performance of the proposed method is also demonstrated through examples. In addition, this method also has the potential to be used for other process synthesis problems that are represented by mixed-integer nonlinear programming (MINLP). Appendix A. MINLP Examples from the Handbook of Test Problems in Local and Global Optimization20 Case 1.
min 2x1 + 3x2 + 1.5y1 + 2y2 - 0.5y3 x,y
s.t. x12 + y1 ) 1.25 x21.5 + 1.5y2 ) 3 x1 + y1 e 1.6 1.333x2 + y2 e 3 -y1 - y2 + y3 e 0 x1,x2 g 0
V. Conclusion This paper applies the tabu search (TS) algorithm and the sequential quadratic programming (SQP) algorithm to solve
y ∈ {0,1}3 The optimal objective function is 7.66726. Case 2.
min -0.7y + 5(x1 - 0.5)2 + 0.8 x,y
s.t. -exp(x1 - 0.2) - x2 e 0 x2 + 1.1y e -1 x1 - 1.2y e 0.2 0.2 e x1 e 1 -2.22554 e x1 e -1 y ∈ {0,1} Figure 8. Optimal structure of HEN case 9.
The optimal objective function is 0.1.
2328
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008
Case 3.
3
∑ (Ci + Ci′Pi)NpiNsizi x,V˘,ω,P,∆p,Np ,Ns ,z i)1 min
min (y1 - 1)2 + (y2 - 2)2 + (y3 - 1)2 x,y
ln(y4 + 1) + (x1 - 1)2 + (x2 - 2)2 + (x3 - 3)2
s.t. P1 - 19.9
s.t. y1 + y2 + y3 + x1 + x2 + x3 e 5 2
2
i i
( )
ω1 3 ω1 2 - 0.1610 V˘ + ωmax ωmax 1 0.000561
y3 + x1 + x2 + x3 e 5.5 2
( )
i
2
y1 + x1 e 1.2 y2 + x2 e 1.8
P2 - 1.21
( )
( )
ω2 3 ω2 2 - 0.0644 V˘ + ωmax ωmax 2
y3 + x3 e 2.5
0.000564
y4 + x1 e 1.2 y2 + x2 e 1.64 2
2
P3 - 6.52
( )
( )
ω3 3 ω3 2 - 0.1020 V˘ + ωmax ωmax 3
y32 + x32 e 4.25 y22 + x32 e 4.64 x1, x2, x3 g 0 y ∈ {0, 1}4 The optimal objective function is 4.55796. Case 4.
min -x1x2x3
0.000232
( ) ( ) ( )
( )
ω2 V˘ 2 ) 0 ωmax 2
( )
ω3 V˘ 2 ) 0 ωmax 3
ω1 2 ω1 - 0.696 V˘ + 0.0116V˘12 ) 0 ωmax ωmax 1
∆p2 - 215
ω2 2 ω2 - 2.950 V˘ + 0.115V˘22 ) 0 ωmax ωmax 2
∆p3 - 316
ω3 2 ω3 - 0.530 V˘ + 0.00946V˘32 ) 0 ωmax ωmax 3 x 1 + x2 + x3 ) 1
s.t. x1 + 0.1y10.2y20.15y6 ) 1
∆Ptotzi - ∆PiNsi ) 0,
x2 + 0.05y40.2y50.15y3 ) 1
i ) 1, 2, 3
V˘iNpi - xiVtot ) 0,
x3 + 0.02y70.08y8 ) 1
i ) 1, 2, 3
ωi - ωmaxzi e 0,
-y1 - y2 - y3 e -1
Pi - Pmax,izi e 0,
i ) 1, 2, 3
i ) 1, 2, 3
-y7 - y8 e -1
V˘i - Vvotzi e 0,
3y1 + y2 + 2y3 + 3y4 + 2y5 + y6 + 3y7 + 2y8 e 10
Npi ∈ {0,1,2,3},
i ) 1, 2, 3
0 e x1,x2,x3 e 1
Nsi ∈ {0,1,2,3},
i ) 1, 2, 3
yi ⊂ {0,1}
ω1 V˘ 2 ) 0 ωmax 1
∆p1 - 629
x,y
-y4 - y5 - y6 e -1
( ) ( ) ( )
( )
i ) 1, 2, 3
8
The optimal objective function is -0.94347. Case 5.
min 3y - 5x x,y
s.t. 2y2 - 2y0.5 - 2x0.5y2 + 11y + 8x e 39 -y + x e 3 2y + 3x e 24 1 e x e 10 y ∈ {1,2,3,4,5,6} The optimal objective function is -17. Case 6. This example is a three-level pump system.21-23 The objective is to determine the most effective cost configuration of a pump system that can achieve a specified increase in pressure. The nonconvex participation of integer variables leads to the presence of many local minima.20 The model is formulated as follows:
zi ∈ {0,1},
i ) 1, 2, 3
The optimal objective function is 122894. Appendix B. Model and Case Descriptions of the Nested TS and SQP Method, Applied to Three HENS Cases, Based on the Model Proposed by Yee and Grossmann8,9 The Yee and Grossmann model, reported in 1990,8,9 can be described as follows:
∑ ∑ ∑ CFi,jZi,j,k + i∈HP ∑ CFcu,iZcu,i + ∑ CFhu,jZhu,j + i∈HP ∑ CcuQcu,i + i∈CP ∑ ChuQhu,j + j∈CP
min
T,Q,∆T,Z i∈HP j∈CP k∈SI
(
(
∑ ∑ ∑ Ci,j U i∈HP j∈CP k∈SI Qcu,i
Ui,cu(LMTD)cu,i
)
)
β
Qi,j,k
i,j(LMTD)i,j
β
+
∑
j∈CP
(
Cj,hu
+
∑ Ci,cu ×
i∈HP
Qhu,j
Uj,hu(LMTD)hu,j
)
β
(B.1)
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008 2329 Table B1. Stream Data for HEN Case 7
∑ ∑ Qi,j,k + Qcu,i ) (Tin,i - Tout,i)Fcpi,
i ∈ HP (B.3.1)
∑ ∑ Qi,j,k + Qhu,j ) (Tout,j - Tin,j)Fcpj,
j ∈ CP (B.3.2)
stream
Tin (K)
Tout (K)
FCP (kw/K)
k∈SI j∈CP
H1 H2 C1 C2 steam water
443 423 293 353 450 293
333 303 408 413 450 313
30 15 20 40
k∈SI i∈HP
(2) Heat balance at each stage:
(Ti,k - Ti,k+1)Fcpi -
∑ Qi,j,k ) 0,
k ∈ SI, i ∈ HP (B.4.1)
∑ Qi,j,k ) 0,
k ∈ SI, j ∈ CP (B.4.2)
j∈CP
Table B2. Parameters for HEN Case 7 parameter CFi,j Ci,j Ui,j CHU
value
parameter
$6250/yr $83.26/yr 0.8 kw/ (m2 K) $80/(kw/yr)
CFcu,i Ci,cu Ui,cu CCU
value $6250/yr $83.26/yr 0.8 kw/ (m2 K) $20/(kw/yr)
parameter
value
CFhu,j Cj,hu Uj,hu
$6250/yr $99.91/yr 1.2 kw/ (m2 K) 1K
∆Tmapp
(Tj,k - Tj,k+1)Fcpj -
i∈HP
(3) Assignment of superstructure inlet temperatures:
Ti,1 ) Tin,i,
i ∈ HP
Tj,NS+1 ) Tin,j,
(B.5.1)
j ∈ CP
(B.5.2)
Table B3. Stream Data for HEN Case 8 stream
Tin (K)
Tout (K)
FCP (kw/K)
h (kw/(m2 K))
H1 H2 H3 C1 C2 steam water
159 267 343 26 118 300 20
77 80 90 127 265 300 60
2.285 0.204 0.538 0.933 1.961
0.1 0.04 0.5 0.01 0.5 0.05 0.2
Ti,k - Ti,k+1 g 0, k ∈ SI, i ∈ HP Tj,k - Tj,k+1 g 0,
Tout,j g Tj,1,
(B.6.2)
i ∈ CP
(B.6.3)
j ∈ CP
(B.6.4)
(5) Hot and cold utility load:
parameter
value
parameter
value
parameter
value
CFi,j Ci,j CHU
$7400/yr $80/yr $110/(kw/yr)
CFcu,i Ci,cu CCU
$7400/yr $80/yr $10/(kw/yr)
CFhu,j Cj,hu ∆Tmapp
$7400/yr $80/yr 1K
FcpiTi,NS+1 - Qcu,i ) FcpiTout,i, FcpjTj,1 + Qhu,j ) FcpjTout,j,
i ∈ HP (B.7.1) j ∈ CP
(B.7.2)
(6) Logical constraints:
Table B5. Stream Data for HEN Case 9 stream
Tin (K)
Tout (K)
FCP (kw/K)
H1 H2 H3 H4 H5 C1 C2 C3 C4 C5 steam water
433 522 544 500 472 355 366 311 333 389 509 311
366 411 422 339 339 450 478 494 433 495 509 355
8.79 10.55 12.56 14.77 17.73 17.28 13.90 8.44 7.62 6.08
max Zi,j,k e 0, Qi,j,k - Qi,j,k
k ∈ SI, i ∈ HP, j ∈ CP (B.8.1)
max Zcu,i e 0, Qcu,i - Qcu,i
i ∈ HP
(B.8.2)
max Zhu,j e 0, Qhu,j - Qhu,j
j ∈ CP
(B.8.3)
(7) Minimum approach temperature constraints: max max Zi,j,k e ∆Ti,j , ∆Ti,j,k - Ti,k + Tj,k + ∆Ti,j k ∈ SI, i ∈ HP, j ∈ CP (B.9.1) max max ∆Ti,j,k+1 - Ti,k+1 + Tj,k+1 + ∆Ti,j Zi,j,k e ∆Ti,j , k ∈ SI, i ∈ HP, j ∈ CP (B.9.2)
where LMTD is approximated using the Chen approximation,
[
(B.6.1)
k ∈ SI, j ∈ CP
Tout,i e Ti,NS+1,
Table B4. Parameters for HEN Case 8
LMTD )
(4) Feasibility of temperatures:
]
∆T1∆T2(∆T1 + ∆T2) 2
1/3
(B.2)
max max ∆Tcu,i - Ti,NS+1 + ∆Tcu,i Zcu,i e ∆Tcu,i - Tout,cu, i ∈ HP (B.9.3) max max ∆Thu,j + Tj,1 + ∆Thu,j Zhu,j e ∆Thu,j + Tout,hu,
∆Ti,j,k e EMAT, and the constraints are defined as follows. (1) Overall heat balance for each stream:
j ∈ CP (B.9.4)
k ∈ SI, i ∈ HP, j ∈ CP
(B.9.5)
where EMAT is the minimum approach temperature.
Table B6. Parameters for HEN Case 9 parameter CFi,j Ci,j Ui,j CHU
value 0 $145.63/yr 0.852 kw/(m2 K) 0
parameter CFcu,i Ci,cu Ui,cu CCU
value 0 $145.63/yr 0.852 kw/(m2 K) $18.12/kw/yr
parameter CFhu,j Cj,hu Uj,hu ∆Tmapp
value 0 $145.63/yr 0.852 kw/(m2 K) 10 K
2330
Ind. Eng. Chem. Res., Vol. 47, No. 7, 2008
Case 7. This problem20 involves two hot streams, two cold streams, one hot utility, and one cold utility, with three stages and no stream splitting. Equality B.2 is replaced by LMTD ) (∆T1 + ∆T2)/2 and the value of β in expression B.1 is 1. Its stream data and parameters are presented in Tables B.1 and B.2. The global optimal solution to this problem consists of five heat exchangers (see Figure 6). The annual cost is $74 711. Case 8. This problem24 involves three hot streams, two cold streams, one hot utility, and one cold utility with three stages and no stream splitting. Equality B.2 is replaced by LMTD ) (∆T1 + ∆T2)/2 and the value of β in expression B.1 is 1. Its stream data and parameters are presented in Tables B.3 and B.4. The overall heat transfer coefficient (U) are computed with film heat transfer coefficient (h):
Ui,j )
(
)
1 1 + hi hj
-1
(B.10)
The global optimal solution to this problem consists of five heat exchangers (see Figure 7). The annual cost is $82 043. Case 9. This problem25 involves five hot streams, five cold streams, one hot utility, and one cold utility with three stages. The value of β in expression B.1 is 0.6. Its stream data are presented in Tables B.5 and B.6. The global optimal solution to this problem consists of ten heat exchangers (see Figure 8). The annual cost is $43 330. Nomenclature Data ) vector storing information about the stream data and parameters of a HEN CP ) set of cold process streams HP ) set of cold process streams NS ) number of stages Zijk ) existence of the match between hot stream i and cold stream j at stage k Zcu,j ) existence of the match between hot stream i and a cold utility Zhu,j ) existence of the match between cold stream j and a hot utility Qijk ) heat exchanged between hot stream i and cold stream j Qcu,i ) heat exchanged between hot stream i and a cold utility Qhu,j ) heat exchanged between cold stream j and a hot utility Tik ) temperature of hot stream i at stage k Tjk ) temperature of cold stream j at stage k ∆Tijk ) temperature approach for the match of hot stream i and cold stream j at stage k ∆Tcu,i ) temperature approach for the match of hot stream i and a cold utility ∆Thu,j ) the temperature approach for the match of cold stream j and a hot utility Nx ) number of continuous variables Nneq ) number of inequality constraints Neq ) number of equality constraints CN ) model condition number Acknowledgment We gratefully acknowledge the financial support of National Key Basic Research and Development Program (No. 2002CB312200), National Natural Science Foundation of China
(No. 60704029), and Zhejiang Provincial Natural Science Foundation of China (No. Y406080). Literature Cited (1) Borchers, B.; Mitchell, J. E. An improved branch and bound algorithm for mixed-integer nonlinear programming. Comput. Oper. Res. 1994, 32, 359-367. (2) Geoffrion, A. M. Generalized benders decomposition. J. Optim. Theory Appl. 1972, 10 (4), 237-260. (3) Duran, M. A.; Grossmann, I. E. An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math. Prog. 1986, 36, 307. (4) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by simulated annealing. Science 1983, 220, 671-680. (5) Holland, J. H. Adaptations in Natural and Artificial Systems; University of Michigan Press: Ann Arbor, MI, 1975. (6) Glover, F. Tabu SearchsPart I. INFORMS J. Comput. 1989, 1 (3), 190-206. (7) Glover, F. Tabu SearchsPart II. INFORMS J. Comput. 1990, 2 (1), 4-32. (8) Yee, T. F.; Grossmann, I. E. Simultaneous optimization models for heat integrationsI. Area and energy targeting and modeling of multi-stream exchangers. Comput. Chem. Eng. 1990, 14 (10), 1151-1164. (9) Yee, T. F.; Grossmann, I. E. Simultaneous optimization models for heat integrationsII Heat exchanger network synthesis. Comput. Chem. Eng. 1990, 14 (10), 1165-1184. (10) Ciric, A. R.; Floudas, C. A. Heat exchanger network synthesis without decomposition. Comput. Chem. Eng. 1991, 15, 358. (11) Furman, K. C.; Sahinidis, N. V. Computational complexity of heat exchanger network synthesis. Comput. Chem. Eng. 2001, 25 (9-10), 13711390. (12) Athier, G.; Floquet, P.; Pibouleau, L.; Domenech, S. Synthesis of heat exchanger network by simulated annealing and NLP procedures. AIChE J. 1997, 43 (11), 3007-3019. (13) Lewin, R. A Generalized method for HEN synthesis using stochastic optimization. II. Synthesis of cost-optimal networks. Comput. Chem. Eng. 1998, 22 (10), 1387-1405. (14) Lewin, R.; Wang, H.; Shalev, O. A Generalized method for HEN synthesis using stochastic optimization: General framework and MER optimal synthesis. Comput. Chem. Eng. 1998, 22 (10), 1503-1513. (15) Lin, B.; Miller, D. C. Solving heat exchanger network synthesis problems with Tabu Search. Comput. Chem. Eng. 2004, 28, 1451-1464. (16) Floudas, C. A. Nonlinear and Mixed-Integer Optimization:Fundamentals and Applications; Oxford University Press: New York, 1995. (17) Dowsland, K. A. Nurse scheduling with Tabu Search and strategic oscillation. Eur. J. Oper. Res. 1998, 106, 393-407. (18) Gendreau, M.; Laporte, G.; Semet, F. A Tabu Search heuristic for the undirected selective traveling salesman problem. Eur. J. Oper. Res. 1998, 106, 539-545. (19) Optimization Toolbox for Use with MATLAB, User’s Guide. Version 3; The Mathworks, Inc.: Natick, MA, 2004. (20) Floudas, C. A. et al. Handbook of Test Problems in Local and Global Optimization; Kluwer Academic Publishers: London, 1999. (21) Westerlund, T.; Petersson, F.; Grossmann, I. R. Optimization of pump configurations as a MINLP problem. Comput. Chem. Eng. 1994, 18, 845-858. (22) Hangandi, V.; Nikolaou, M. A hybrid approach to global optimization using a clustering algorithm in a genetic search framework. Comput. Chem. Eng. 1998, 22, 1913-1925. (23) Adjiman, C. S.; Androulakis, I. P.; Floudas, C. A. Global optimization of mixed-integer nonlinear problems. AIChE J. 2000, 46 (9), 17691798. (24) Zamora, J. M.; Grossmann, I. E. A global MINLP optimization algorithm for the synthesis of heat exchanger networks with no stream splits. Comput. Chem. Eng. 1998, 22 (3), 367-384. (25) Linhoff, B.; Flower, J. R. Synthesis of heat exchanger networks: Systematic generation of energy optimal networks. AIChE J. 1978, 24 (4), 633-654.
ReceiVed for reView September 13, 2007 ReVised manuscript receiVed December 18, 2007 Accepted January 9, 2008 IE071245O
