Object-Oriented Disjunctive Programming with a Nested Heuristic and

Ind. Eng. Chem. Res. 2010, 49, 1779–1791

1779

Object-Oriented Disjunctive Programming with a Nested Heuristic and Gradient-Based Solver for Chemical Process Synthesis Daqing Tian,† Lingyu Zhu,‡ Xi Chen,*,† Zhijiang Shao,† and Jixin Qian† State Key Lab of Industrial Control Technology, Department of Control Science and Engineering, Zhejiang UniVersity, Hangzhou 310027, China, and College of Chemical Engineering and Materials Science, Zhejiang UniVersity of Technology, Hangzhou 310014, China

The generalized disjunctive programming (GDP) model has been proposed and applied in the past decade as an alternative to the mixed integer nonlinear programming (MINLP) model because it has the advantages of being straightforward when used in conditional modeling and being able to reduce the complexity in the sub-NLP. In this paper, we introduced an improved variant of the traditional GDP model, otherwise known as the object-oriented disjunctive programming (ODP) model. Such a method helps generate well-posed subNLP, thereby improving the solving process. A nested method combining the heuristic algorithm and gradientbased optimizer is also proposed to solve the GDP and ODP. It is a two-layer method, wherein a heuristic algorithm performs master iterations in the outer-loop when dealing with the Boolean variables, and a gradientbased NLP solver is applied in the inner-loop when dealing with the sub-NLP. Excellent performance has been demonstrated by applying the modeling and solving methods into the process synthesis of heat exchanger networks (HENs) and water networks (WNs). 1. Introduction Optimization-based approaches have been widely used as systematic methods for chemical process synthesis. The common strategy is to formulate the synthesis problem as a mathematical model, usually with the aid of a superstructure, and then solve it with coded optimization solvers. When discrete decisions have to be made, the mixed integer nonlinear programming (MINLP) model is required to address the optimization problem, wherein a set of integer variables (usually binary ones) is incorporated to represent the discrete values. In recent years, researchers introduced the idea of conditional modeling into the chemical process synthesis and built generalized disjunctive programming (GDP) models1 to comprise an alternative to the MINLP. The main characteristic of GDP is the adoption of a set of logic variables, which not only represents the discrete decisions but also accomplishes the conditional modeling of equations depending on the specific values they take. A variety of chemical process problems have been modeled using GDPs.2,3 It was pointed out that GDP models are more robust and computationally efficient compared to MINLP.4 Some deterministic solvers of MINLPs have been modified to solve GDPs. For instance, logic-based OA1 has been applied to solve a number of process network problems modeled in GDP. It requires the solution of the inner-layer sub-NLP generated by fixing the logic values and of outer-layer disjunctive OA master problems that provide new values for the Boolean variables. Similar to the problems in MINLP, the computation time of these methods increases exponentially with the size of GDPs. Moreover, the global optimality cannot be guaranteed when the model is nonconvex. Various techniques have been proposed to maintain global optimality, such as the convex underestimators.5 On the other hand, research and applications of heuristic algorithms (also known as stochastic algorithms) have been widely studied in past decades. Methods that can be categorized into heuristic algorithms include simulated annealing,6 Taboo search,7,8 and genetic algorithm,9 * To whom correspondence should be addressed. Tel.: +86-57187953068. Fax: +86-571-87951068. E-mail: [email protected]. † Zhejiang University. ‡ Zhejiang University of Technology.

among others. In contrast to the deterministic methods, the performance of heuristic algorithms is proven to be robust for NP-hard and even discontinuous problems. Moreover, these algorithms have the advantage of avoiding entrappment in a local optimum and reaching the global optimal solution or a solution close to it. A number of papers have been published on their application in process synthesis.10-13 Compared with the deterministic methods, however, the convergence speed of heuristic methods is normally slow. In this paper, a hybrid method is proposed to solve the GDP models. It adopts a heuristic algorithm in the outer-layer searching for the best Boolean values and a deterministic gradient-based solver in the inner-layer which can optimize reduced sub-NLP. However, before the algorithm implementation, a novel objectoriented disjunctive programming (ODP) model, as a variant of the original GDP, is proposed for the modeling of the process synthesis. This model helps generate more concise and well-posed sub-NLP which facilitates the solving of the problems. Finally, both the modeling and the nested solving method are applied to two typical types of chemical process synthesis, the heat exchanger network, and water network syntheses. 2. Object-Oriented Disjunctive Programming Model 2.1. Formulation of the GDP Model. Under some circumstances, disjunctions can involve more than two terms. In this paper, we focus on the disjunctions involving two terms only, as that is normally enough for the chemical process synthesis. Therefore, the general form of a GDP model is the following: min Z )

∑c

k

+ f(x)

k

s.t.

[

][

]

g(x) e 0 Yk ¬Yk h1k(x) e 0 ∨ h2k(x) e 0 k ∈ D ck ) γk ck ) 0 Ω(Y) ) true x ∈ Rn, c ∈ Rm, Y ∈ {true, false}m

10.1021/ie901010a CCC: $40.75  2010 American Chemical Society Published on Web 12/28/2009

(1)

1780

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010

where the objective function is usually linear with the variables ck representing the costs of relevant units. Meanwhile, g(x) is a group of functions which holds the mass and enthalpy balance equations and design specifications irrespective of the Boolean variables’ values. The second type of constraints is described in the disjunctive form in which Y is a set of Boolean variables denoting the existence of units, selection of operating conditions, and other parameters. These Boolean variables are the key operators of conditional modeling. If Yk is true, the corresponding first block of equations is activated; otherwise, the second block of equations is activated. In the chemical process synthesis, the second block of equations, as shown later in the example, is usually in the form where all the continuous variables are set to zero when the relevant unit does not exist. Finally, the term Ω(Y) ) true is the logical relation which must be satisfied by the Boolean variables. 2.2. Sub-NLP of the GDP Model. The GDP model is normally based on a superstructure. Once the Boolean variables are specified, a deterministic subnetwork is generated from the postulated superstructure. The GDP is also branched into a reduced NLP problem. An illustrative example is adopted from Tu¨rkay and Grossmann’s paper2 with a minor modification. It is an optimization problem of a simple process network with minimized capital and operating cost. The superstructure is illustrated in Figure 1. The raw material A can be processed in units Uab1 and Uab2 to produce B; or B can be purchased from the market. B is needed in units Ubc3 and Ubc4 for the production of C. The flow rates of each stream are shown in Figure 1. The GDP formulation of this example is described as min Z ) c1 + c2 + c3 + c4 + 7Fb - 11FC s.t. FA - Fin - Fin ) 0 1 2 FB - Fout1 - Fout2 - Fb ) 0 FB - Fin3 - Fin4 ) 0 FC - Fout3 - Fout4 ) 0 Fout2 e 5 Fout3 e 1 Uab1 ¬Uab1 Fout1 ) ln(1 + Fin1) ∨ Fin1 ) Fout1 ) 0 c1 ) 1 + 1.8Fin1 + Fout1 c1 ) 0 Uab2 ¬Uab2 Fout2 ) 1.2 ln(1 + Fin2) ∨ Fin2 ) Fout2 ) 0 c2 ) 1.5 + 1.8Fin2 + 1.2Fout2 c2 ) 0 Ubc3 ¬Ubc3 Fout3 ) 0.9Fin3 ∨ Fin3 ) Fout3 ) 0 c3 ) 3.5 c3 ) 0 Ubc4 ¬Ubc4 Fout4 ) 0.7Fin4 ∨ Fin4 ) Fout4 ) 0 c4 ) 2.8 c4 ) 0

[ [

[ [

][ ][

c1, c2, c3, c4 g 0;

][ ] ][ ]

Figure 1. Superstructure of the illustrative example.

Figure 2. Substructure of the illustrative example.

The sub-NLP generated from the GDP becomes min Z ) c1 + c2 + c3 + c4 + 7Fb - 11FC s.t. FA - Fin - Fin ) 0 1 2 FB - Fout1 - Fout2 - Fb ) 0 FB - Fin3 - Fin4 ) 0 FC - Fout3 - Fout4 ) 0 Fout2 e 5 Fout3 e 1 Fout1 ) 0 Fout1 ) 0 c1 ) 0 Fout2 ) 1.2 ln(1 + Fin2) c2 ) 1.5 + 1.8 Fin2 + 1.2Fout2 Fout3 ) 0.9Fin3 c3 ) 3.5 Fin2 ) 0 Fout4 ) 0 c4 ) 0 c1, c2, c3, c4 g 0;

(3)

Fink, Foutk g 0, k ) 1, ..., 4; FA, FB, FC, Fb g 0

] ]

Fink, Foutk g 0, k ) 1, ..., 4;

FA, FB, FC, Fb g 0 Uabi, Ubcj ∈ {true, false}, i ) 1, 2, j ) 3, 4

(2) When the Boolean variables are specified, e.g., Uab2 and Ubc3 are true and Uab1 and Ubc4 are false, the corresponding substructure is determined as shown in Figure 2.

The sub-NLP includes the same number of continuous variables as the GDP model. The variables belonging to the nonexisting units still exist in the form of zero terms. Compared with the ideal case where these zero variables should be eliminated, this NLP does not reduce the dimension of the continuous variables. Worse, in cases where these zero variables are present in bilinear or other nonlinear terms, the sub-NLP problem becomes difficult or even impossible to solve using specific solvers. Considering this drawback, we therefore propose a variant of the traditional GDP. 2.3. Introduction of the ODP Model. Generally, a GDP model consists of two types of variables, the Boolean and the continuous variables. In the chemical process synthesis, a Boolean variable represents the existence of a unit, and the continuous variables represent the operating conditions such as heat and mass. Most of the continuous variables can be categorized as one unit; therefore, a Boolean variable, combined with some continuous variables, can constitute an object. Aside from the variables, certain


equations (relations of these continuous variables) can also be viewed as members of this object. This presents great similarity to the object-oriented programming language wherein a class possesses two types of members: data members and member functions. When an object of a certain class is created, all the data members belonging to this object will be assigned a memory in the computer, and the member functions can be performed on the data members. Similarly, in the process synthesis, when one process object is selected (i.e., the corresponding Boolean variable is true), the relevant continuous variables belonging to this object will be present in the sub-NLP including their correlations (equations); otherwise, these continuous variables and their correlations should disappear. Therefore, the Boolean variables should not only decide the selection of the equations, but also the existence of some continuous variables. We propose a variant of the traditional GDP, the socalled object-oriented disjunctive programming (ODP) model, which facilitates the generation of concise sub-NLPs by eliminating the redundant variables and equations. Taking the example above for illustration, its ODP formulation can be written as:

min Z )

∑

i,Uabi.Y)true

s.t. FA -

∑

Uabi.c +

FA -

[

Uab1.Y ) true Uab1.Fout ) ln(1 + Uab1.Fin) Uab1.c ) 1 + 1.8Uab1.Fin + Uab1.Fout Uab1.Fin, Uab1.Fout, Uab1.c g 0

FB - Fb - Uab2.Fout ) 0 FB - Ubc3.Fin ) 0 FC - Ubc3.Fout ) 0 Uab2.Fout ) 1.2 ln(1 + Uab2.Fin) Uab2.c ) 1.5 + 1.8Uab2.Fin + 1.2Uab2.Fout Uab2.Fout e 5 Uab2.Fin, Uab2.Fout, Uab2.c g 0 Ubc3.Fout ) 0.9Ubc3.Fin Ubc3.c ) 3.5 Fout3 e 1 Ubc3.Fin, Ubc3.Fout, Ubc3.c g 0 FA,FB, Fb, FC g 0

Uabi.Fout ) 0

∑

Ubcj.Fin ) 0

∑

Ubcj.Fout ) 0

j,Ubcj.Y)true

[

FC -

j,Ubcj.Y)true

]

Uab1.Y ) true Uab1.Fout ) ln(1 + Uab1.Fin) Uab1.c ) 1 + 1.8Uab1.Fin + Uab1.Fout Uab1.Fin, Uab1.Fout, Uab1.c g 0 Uab2.Y ) true Uab2.Fout ) 1.2ln(1 + Uab2.Fin) Uab2.c ) 1.5 + 1.8Uab2.Fin + 1.2Uab2.Fout Uab2.Fout e 5 Uab2.Fin, Uab2.Fout, Uab2.c g 0 Ubc3.Y ) true Ubc3.Fout ) 0.9Ubc3.Fin Ubc3.c ) 3.5 Fout3 e 1 Ubc3.Fin, Ubc3.Fout, Ubc.3c g 0 Ubc4.Y ) true Ubc4.Fout ) 0.7Ubc4.Fin Ubc4.c ) 2.8 Ubc4.Fin, Ubc4.Fout, Ubc4.c g 0

[ [

[

]

]

(6)

min Z ) Uab2.c + Ubc3.c + 7Fb - 11FC s.t. FA - Uab2.Fin ) 0

Ubcj.c + 7Fb - 11FC

i,Uabi.Y)true

FB -

]

will be present in the sub-NLP if Uab1.Y ) true holds. When Uab2 and Ubc3 are assigned to be true and Uab1 and Ubc4 are false, the generated sub-NLP is as follows:

Uabi.Fin ) 0

∑

(5)

the variable denoting the inlet flow rate of unit Uabi will be present in the equation only if the condition “Uabi.Y ) true” is satisfied. In a superstructure-based model, the equations for a mixer or a splitter are usually in a summation form including the continuous variables of different units connected to this node. The sum-if notation implements this dynamic selection. Meanwhile, the original “if-else-end” conditional modeling form in GDP is transformed into an “if-end” mode. The equations

i,Uabi.Y)true

FB - Fb -

Uabi.Fin ) 0

i,Uabi.Y)true

j,Ubcj.Y)true

∑

∑

1781

(7)

Compared with the sub-NLP derived from GDP in (3), the above NLP is more concise since it only includes variables associated with the existing units, thereby reducing the model scale. This means that the dimension of the sub-NLP is dynamically varying according to the specified Boolean values. This is implemented by applying the object-oriented programming idea in the ODP model wherein the objects are constituted by variables and equations. When an object does not exist, the relevant variables disappear. It should also be noted that sometimes, there are static variables which do not belong to any objects, such as FA in this example. Figure 3 illustrates a comparison of the variable definitions in GDP and ODP. Extended from the traditional GDP model, the general formulation of the ODP model can be written as follows:

]

FA, FB, Fb, FC g 0, Uabi.Y, Ubcj.Y ∈ {true, false},

min Z )

i ) 1, 2, j ) 3, 4

(4) where Uab1, Uab2, Ubc3, and Ubc4 represent four objects associated with the four units. We introduce the operator “.” to indicate the “belongs-to” relationship between an object and its variables. The existence of the object Unitk is denoted by Unitk.Y. Here, we define a symbol, ∑k,Unitk.Y)true Unitk.x, for a “sum-if” operation. It is a summation operation of the continuous variable Unitk.x with respect to the index k wherein the continuous variable is added into the summation only if Unitk.Y ) true holds. For example, in the following equation taken from eq 4

∑

Unitk.c + f(xs)

k,Unitk.Y)true

s.t. g(

[

∑

k,Unitk.Y)true

Unitk.xd, xs) e 0

] [

]

Unitk.Y ) true Unitk.Y ) false ∨ ,k ∈ D hk(Unitk.xd, xs) e 0 hk(Unitk.xd, xs) e 0 Unitk.Y ) true hk(Unitk.xd, Unitk.c, xs) e 0 , k ∈ U Unitk.xd, Unitk.c g 0

[

]

Ω(Unitk.Y) ) true, Unitk.Y ∈ {true, false}m, k ∈ U

(8)

1782

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010 Table 2. Comparison between the GDP and ODP in the Global Optimality of Sub-NLPs GAMS case HENS1 HENS2 WNS2

Figure 3. Comparison of variables in GDP and ODP. Table 1. Comparison between the GDP and ODP in the Sub-NLP Scale case

subnetwork

GDP

ODP

HENS1

random optimal random optimal random optimal

79 79 201 201 88 88

35 31 57 53 51 49

HENS2 WNS2

where xs represents the static variables; Unitk.c, Unitk.xd, and Unitk.Y, belonging to a specific unit object (Unitk), denote the cost, the dynamic variables, and the existence variable of the unit, respectively. At the same time, U and D represent the unit object sets and D is usually a subset of U or sometimes a void set; while g is the compromised value of enthalpy and mass balance equilibrium equations which describe the relationship among different objects. The disjunctive parts, which may not be necessary in many cases of process synthesis, reveal the balance equations described in a disjunctive form. The objectoriented part describes the specific unit’s constraints which are included or excluded dynamically in the sub-NLPs according to the relevant unit’s existence. Finally, Ω(Unitk.Y) ) true describes the logical relations among the Boolean variables. 2.4. Sub-NLP Comparison of GDP and ODP. In the process of solving a GDP, an NLP solver will be frequently called to solve the generated sub-NLPs when the Boolean variables are specified in an outside loop. Compared with GDP, the advantage of the ODP lies in its ability to reduce the scale of both the equations and variables of the sub-NLP problems. Moreover, as the zero variables related to the nonexisting unit objects have been removed in the ODP, it avoids the singularities and facilitates the NLP solving process. Many NLP solvers are applicable for this case. On the other hand, the sub-NLPs generated by a GDP model are relatively difficult to solve and highly dependent on the performance of the NLP solver. A comparison between the sub-NLPs of the GDP and ODP models is conducted. Two examples of the heat exchanger network synthesis problem, HENS1 and HENS2, and one example of a water network synthesis problem, WNS2, were tested for the sub-NLP comparisons. The detailed descriptions of these examples can be found in sections 3 and 4. The comparison results are shown in Table 1. Two sub-NLPs were tested for each case, one generated from a random structure and another from the optimal structure. The last two columns of the table present the number of the continuous variables of the generated sub-NLPs. For the sub-NLPs generated from a different structure, the traditional GDP does not change in the variable scale. On the contrary, the ODP generates sub-NLPs with different scales (i.e., the variable number are dynamically varying with the structure). As the ODP is capable of removing the redundant variables in the generated sub-NLPs, it results in consistently smaller number of variables than those produced in the GDP. For a large scale case like HENS2, the reduction

MATLAB

subnetwork CONOPT MINOS FMINCON random optimal random optimal random optimal

100/100 100/100 98/100 100/100 99/100 87/88

85/99 100/100 0/100 36/100 20/89 11/99

91/100 92/99 57/93 100/100 44/76 35/58

time cost(s) 168.5/21.4 54.1/12.0 1382.5/42.4 223.6/19.6 440.1/173.6 366.5/139.9

can be dramatic: the number of variables in the ODP case is less than one-third of the GDP case. The next step involves the application of different NLP solvers including GAMS/CONOPT, GAMS/MINOS, and MATLAB/FMINCON, in solving the generated sub-NLPs. CONOPT implements the reduced gradient method with restoration while MINOS implements the reduced gradient method without restoration. FMINCON applies the sequential quadratic programming (SQP) method. The first two tests were conducted on GAMS14 (distribution 22.8), while the third was on MATLAB. The three tests were all processed on an Intel 3G HZ, 3.0 G memory machine with Windows XP. For each case, 100 runs of the sub-NLPs with different initial points were solved by the three NLP solvers separately. The results are shown in Table 2. As can be seen, the third, fourth, and fifth columns of the table show the global optimality results of different solvers. In each grid, the first number represents the run number among 100 runs to reach the global optimal point for the GDP, while the second number is the result of the ODP. In the first HEN case, as the scale is not large and the cost function of a heat exchanger is linear with its area (β ) 1), the performance of each solver on the sub-NLPs does not change much. In the second HEN case, however, the cost function of a heat exchanger is nonlinear with its area (β ) 0.6), resulting in the high nonlinearity in the sub-NLPs. In addition, the scale of this case is also large. For this reason, the performance of the different solvers was dramatically different for the GDP. On the other hand, the performance of different solvers in the ODP formulation was consistently good. The third case is a large scale water network synthesis example. With bilinear terms in constraints, the performance of the NLP solvers was worse than the HEN cases. Nevertheless, we can still observe the obvious advantage of the ODP. The figures in the last column in Table 2 present the overall time cost of 100 runs in MATLAB of the GDP and ODP. An obvious improvement is observed in the time cost of the sub-NLPs generated from the ODP compared to the time cost of the GDP. In conclusion, the proposed ODP generates concise and wellconditioned sub-NLPs, thereby improving the solving process of the sub-NLPs. When a superstructure-based model is highly dimensioned and complicated, the ODP model can be especially efficient. 3. Nested Heuristic and Gradient-Based Solver 3.1. Decomposition Strategy for the ODP. The logic-based OA method consists of a cycle of major iterations. A master disjunctive linear problem in the outer-layer was solved in each major iteration in order to produce a new set of Boolean variables, while a sub-NLP was also solved in the inner-layer. The algorithm will terminate if the upper and lower bound provided by the two layers will fall within a specified tolerance. Similarly, using the master and slave loop strategy, we propose a hybrid algorithm to overcome the local optimality disadvantage of the deterministic methods. The decomposition strategy is


1783

the next one. This iterative update is accomplished by selecting the best but not tabooed solution from neighbors of the current solution. Thus, the neighborhood generation plays an important role in the performance of the algorithm. A set, Yk, is defined as follows: M

∑ |y′

Yk ) {y′|

m

- ym | ) k}, k ∈ N+

(10)

m)1

where M denotes the dimensions of the Boolean variables, y is the current solution, y′ is one of its neighbors, and the subscript m stands for the mth element of the vector. Both the current solution and its neighbors have M elements denoted by 0 or 1. For the ODP, 1 denotes that an object is selected and 0 means otherwise. The definition of Yk means that its member, y′, is exactly k bits different from the current solution, y. The neighborhood of the current solution, YN, is composed of a number of such sets different from the index k, which controls the difference magnitude from the current solution: Figure 4. Decomposition strategy of the nested method for ODP.

YN ) {Yk |k ) 1, ..., K}

illustrated in Figure 4. In the outer-loop, a heuristic algorithm such as the Taboo search (TS) performs the master iterations, dealing only with the Boolean variables. The heuristic algorithm solves a combinatorial optimization as follows: min f(Y) s.t. Ω(Y) ) true

(9)

where f(Y) is an implicit function which calls for the solution of the inner-loop. It is responsible for iteratively locating the values of the Boolean variables to escape from being trapped in a local optimum. When a set of Boolean values, Y, denoting the existence of chemical objects is specified, the ODP model can be automatically converted into a reduced sub-NLP with no redundant equations or variables. The generated sub-NLP is illustrated in Figure 4, in which UT denotes the set of selected units specified by the outer-loop. In the inner-loop, a gradient-based NLP solver such as the SQP or the generalized reduced gradient (GRG) methods can be applied to solve the sub-NLPs. The inner-loop optimizer will return the objective function value to the outer-loop if it is successfully converged. Otherwise, a penalty will be given to the objective function in the outer-loop. The success of the nested method depends on two issues: (1) the ability of the heuristic algorithm to search for the globally optimal structure and (2) the convergence of the sub-NLPs. For the first issue, the heuristic algorithm can deal with nonconvex or even discontinuous problems at the cost of longer time consumed. As the process synthesis is normally conducted offline, a reasonable time cost can be accepted. As for the second one, the ODP model helps generate a well-posed sub-NLP which is concisely dimensioned, thereby efficiently avoiding singularity. The detailed implementation of the nested method is discussed below. 3.2. Heuristic Algorithm in the Outer-Layer. Taboo search has been given focused in this paper because of its ability to perform well in dealing with global optimization. It is a metaheuristic method which employs the neighborhood search and adaptive memory to explore the solution space and prevent it from being trapped in a local optimum. In this project, TS is coded in the C++ language. It is a general algorithm independent of the problem, although the parameters are initialized according to the dimensions of the Boolean variables. Some details of the developed algorithm are presented in the following sections. 3.2.1. Neighborhood Generation. The taboo search method seeks for the optimal solution through a sequence of moves; here, a move is an operation that changes from the current solution to

(11)

Generally, k starts from 1 to K, increasing incrementally by 1. The value of K is initialized according to the dimension of the Boolean variables, M, at the beginning of the execution. 3.2.2. Taboo and Aspiration Criteria. The neighborhood generation selects the best neighbors as the updated current solution. To avoid being trapped in a local optimum, a taboo operation was conducted to restrict the search around a local point. 2-Fold taboo lists were then built to perform this task. The first taboo list, TL1, keeps the information of the recently visited solutions. A new neighbor, y′, is tabooed if the following condition holds: M

∃y ∈ TL1,

∑ |y′

m

- ym | e 1

(12)

m)1

The second taboo list, TL2, stores the information of the “local optimal solutions.” If the objective function value has not been improved for L iterations, the best solution in the past L iterations will be considered as a local optimal solution and stored in TL2. Any y′ satisfying the following condition M

∃y ∈ TL2,

∑ |y′

m

- ym | e TP

(13)

m)1

is then tabooed. In this condition, TP is a penalty parameter for TL2, which is dependent on the dimensions of the Boolean variables. In addition to the taboo operation, an aspiration operation was also conducted to release the taboo point satisfying certain conditions. A new solution will be aspired if its objective value is better than the best solution known. Moreover, there is still a probability p for a taboo point to aspire as a candidate for the updated solution. 3.2.3. Termination Criterion and Parameter Specifications. The procedure will terminate if the maximum iteration count, L1, is reached or if the objective function has not been improved for L2 major iterations. In our code, almost all the parameters of the TS such as L1 and L2 as well as the length of two taboo lists, TL1 and TL2, were initialized according to the value of M. The parameter K in eq 11 decides the scale of the neighborhood. It is usually set specifically for each problem, since the logical constraints’ strictness of each ODP model can be very different. While the number of continuous variables and the complexity of the sub-NLP problems influence the time

1784


cost of the inner-loop, the parameters of TS generally determine the number of NLP calculations. 3.3. Implementation of the Nested Method. The program of the nested method was implemented using Microsoft Visual Studio 6.0 (the heuristic part) and GAMS (the inner NLP part). The choice of a specific heuristic and gradient-based NLP algorithm can vary, depending on the problem and the user’s preference. For example, the heuristic algorithm can be taboo search or genetic algorithm, and the NLP optimizer can be equipped with the CONOPT and MINOS methods. In this particular case, we only illustrated the basic procedure of the nested method, with TS as the outer heuristic method. The flowchart of the nested method is illustrated in Figure 5. The detailed procedures are as follows: Step 1: Initialization. At the beginning of the outer-loop, initialize the parameters of the heuristic algorithm based on the dimensions of the Boolean variables. Step 2: Generation of the initial “current solution”. The initial solution can be randomly and repeatedly generated until a feasible configuration is obtained. Another option is to follow a certain experience as described in the latter examples. Step 3: Neighborhood generation. A group of the current solution’s neighbors is produced following eqs 10 and 11. Afterward, each set of the Boolean variables, y′, will be handled sequentially. Step 4: Boolean constraints check. The term y′ is checked whether the Boolean constraint, Ω(y′) ) true, is satisfied or not. If yes, proceed to step 5; otherwise, go to step 9. Step 5: Taboo check. The neighbor y′ is checked whether it is tabooed or not. If not, proceed to step 7; otherwise, go to step 6. Step 6: Aspiration check. The neighbor y′ is further judged whether it can be aspired. If yes, proceed to step 7; otherwise, go to step 9. Step 7: Generation and solving of inner-loop sub-NLPs. Given the Boolean values, the ODP or GDP will be transformed into a sub-NLP problem and then solved by a gradient-based optimizer. Step 8: Updating. Once the results of the inner-loop are received, update the fitness value of y′. Step 9: Neighborhood examination. Examine if every neighbor has been handled. If not, select the next one as y′ and proceed to step 4; otherwise, go to step 10. Step 10: Solution selection. Choose the next current solution. Subsequently, the taboo lists will be updated. Step 11: Termination. Examine if the termination criteria are satisfied. If not, proceed to step 3; otherwise, exit. The proposed method was applied to solve two types of process synthesis problems, adopting TS in the outer-layer along the process. The performance between the GDP and ODP models with two NLP solvers, CONOPT and MNOS, in the inner-layer was also compared. All the executions of the program were conducted on an Intel 2.8G Hz, 3.0 G memory machine with Windows Server 2003. 4. HENS Applications The heat exchanger network synthesis (HENS) has been the most extensively studied problem in process synthesis for the past 40 years. A proposed simultaneous method to approach the HENS problem was to formulate it as an MINLP superstructure model.15-17 A GDP formulation18 was also proposed by Lee and Grossmann. In this GDP model, disjunctions were used only for heaters and coolers, leading to a hybrid MINLP/GDP model. In this paper,

we proposed an ODP model for HENS based on Yee and Grossmann’s MINLP model, coupled with the elimination of some redundant equations. The traditional GDP model with disjunctions for exchangers, coolers, and heaters can be also obtained easily by adaption of this ODP model. The ODP model for the HENS problem is presented as follows: HP

set of hot process streams set of cold process streams set of stages in the superstructure

Y

EMAT

Uniti,cu

the exchanger between hot stream i and cold stream j at stage k the cooler for hot stream i

Unitj,hu

the heater for cold stream j

Tj,k

CP TS Uniti,j,k

the existence of a unit the cost of a unit heat exchanged between two streams minimum approach temperature difference temperature of hot stream i at stage k temperature of cold stream j at stage k

C Q

Ti,k

(1) The objective function min

∑ ∑

∑

Uniti,j,k.C +

i∈HP j∈CP k∈TS,Uniti,j,k.Y)true

∑

∑

Uniti,cu.C +

i∈HP,Uniti,cu.Y)true

(14)

Unitj,hu.C

j∈CP,Unitj,hu.Y)true

(2) Heat balance equations (Ti,k - Ti,k+1)FCpi -

(Tj,k - Tj,k+1)FCpj -

∑

Uniti,j,k.Q ) 0,

∑

i ∈ HP, k ∈ TS Uniti,j,k.Q ) 0,

j∈CP,Uniti,j,k.Y)true

i∈HP,Uniti,j,k.Y)true

j ∈ CP, k ∈ TS

(15) (3) Assignment of the inlet temperatures Ti,1 ) Ti,in, i ∈ HP Tj,NS+1 ) Tj,in, j ∈ CP

(16)

(4) Disjunctive block for the heat balance Uniti,cu.Y ) true

[

FCpi(Ti,NS+1 - Ti,out) - Uniti,cu.Q ) 0

[

Uniti,cu.Y ) false

Ti,NS+1 - Ti,out ) 0 Unitj,hu.Y ) true

[

[

]

Unitj,hu.Y ) false Ti,out - Ti,1 ) 0

]

∨

, i ∈ HP

FCpj(Ti,out - Ti,1) - Unitj,hu.Q ) 0

[

]

]

(17) ∨

, j ∈ CP

(5) Object-oriented block for the exchangers Uniti,j,k Uniti,j,k.Y ) true max Uniti,j,k.Q e Qi,j,k

Ti,k - Tj,k g EMAT Ti,k+1 - Tj,k+1 g EMAT

(

Uniti,j,k.C ) CFi,j + Ci,j

Uniti,j,k.Q, Uniti,j,k.C g 0

[

Ui,j

Uniti,j,k.Q (Ti,k - Tj,k)(Ti,k+1 - Tj,k+1) 2

)

β

]

]

,

i ∈ HP, j ∈ CP, k ∈ TS

(18)

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010 -1

Ui,j ) (1/hi + 1/hj)

Table 3. Data of Case 1

1785

(23)

Stream Data 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

Uniti,j,k.Y ⇒ (¬Unitia,j,k.Y), i, ia ∈ HP, i * ia, j ∈ CP, k ∈ TS Uniti,j,k.Y ⇒ (¬Uniti,ja,k.Y), i ∈ HP, j, ja ∈ CP, j * ja, k ∈ TS

Parameters parameter

value

parameter

value

parameter

value

CFi,j Ci,j CHU

$7,400/y $80/y $110/kW/y

CFcu,i Ci,cu CCU

$7,400/y $80/y $10/kW/y

CFhu,j Cj,hu ∆Tmapp

$7,400/y $80/y 1K

[

(6) Object-oriented block for the coolers Uniti,cu Uniti,cu.Y ) true Uniti,cu.Q e Qmax i,cu Ti,NS+1 - Tcu,out g EMAT Uniti,cu.C ) CFi,cu + CCUQi,cu + Uniti,cu.Q Ci,cu (Ti,NS+1 - Tcu,out)(Ti,out - Tcu,in) Ui,cu 2 Uniti,cu.Q, Uniti,cu.C g 0

(

[

[

)

β

]

]

(7) Object-oriented block for the heaters Unitj,hu Unitj,hu.Y ) true Unitj,hu.Q e Qmax j,hu Thu,in - Tj,1 g EMAT Unitj,hu.C ) CFj,hu + CHUQj,hu + Unitj,hu.Q Cj,hu (Thu,in - Tj,1)(Thu,out - Tj,out) Uj,hu 2 Unitj,hu.Q, Unitj,hu.C g 0

(

[

(8) Variable bounds

)

β

]

]

,

,

i ∈ HP

(19)

j ∈ CP

(20)

Ti,k g 0, Tj,k g 0, Uniti,j,k.Y, Uniti,cu.Y, Unitj,hu.Y ∈ {true, false} i ∈ HP, j ∈ CP, k ∈ TS (21) In the above ODP model for HENS, the AMTD is used as the heat transfer driving force which can also be changed to Chen’s approximation: LMTD ) [∆T1∆T2(∆T1 + ∆T2)/2]1/3

It should be noted that the solution structure is constrained to be without any stream splits, which means the following set of logical constraints should be added to the ODP model:

(22)

In the code of the nested method in solving HENS, the initial configuration is generated by assuming that all process streams exchange heat with utilities. This means that the network consists of only coolers and heaters without any heat exchangers. However, this structure, while always feasible, is far from being the best solution. 4.1. Case 1: HENS1. The first case19 to be studied involves three hot and two cold streams, steam, and cooling water with three stages. The stream data and cost parameters are listed in Table 3. The value of β in the ODP model for HENS is 1 and the AMTD is used. The overall heat transfer coefficients (U) are computed with h:

(24)

As the two most important parameters in TS, L2 is set to be the number of objects, M, and K is half of M. The nested method obtained the global solution 82 042.91, which is the same as the best known solution. The corresponding configuration with the heat amount exchanged is shown in Figure 6. As TS is a stochastic method, a total of 60 runs have been tested. An iterative process of case 1 is shown in Table 4, which shows how the TS iterated in the outer-loop. In this run, the optimal structure represented by the Boolean variables was found at the eighth iteration. Table 5 shows the statistical results of case 1. Two different NLP solvers were applied for the innerloop as a comparison. The GDP and OPD formulations were also compared. The fourth column lists the successful instances where the global solution was reached among all the 60 runs. The fifth column presents an average count of the total NLP runs, which is highly relevant to the computation time. There are two numbers in the sixth column; the first denotes the average iteration number when the best solution is reached for the first time, and the second represents the total major iteration count. The last column gives the average time cost of all executions, including the file interface I/O as well as other executions. It can be observed that the nested method can reach global solution for both models. However, the ODP model is advantageous over ODP in the sense that it requires less computation time. The advantage of the ODP model will be more obvious in the following complex case. 4.2. Case 2: HENS2. The second case20 involves five hot streams, five cold streams, one hot utility, and one cold utility with two stages. The value of β in the model is 0.6, and Chen’s approximation for LMTD is adopted. The stream data and parameters are listed in Table 6. The TS parameters were set similarly as in case 1. L2 was set to be identical to M. K, the upper bound of k, was half of M, and when k is greater than 10, an increment of 2 other than 1 was adopted. For this case, the total number of objects is as large as 60, causing difficulty in solving the problem. The best known solution12 is 43 329 when the stage number is 2. Our nested method obtained the identical structure with the same operating conditions as shown in Figure 7. However, the objective function value is slightly different (43 330.50). Meanwhile, a suboptimal solution with an objective value of 43 331.37, which is very close to the global optimal solution, was frequently obtained in our tests. From the engineering viewpoint, the objective function values were so close that both were deemed acceptable. The structure corresponding to this suboptimal solution is shown in Figure 8. The statistical results of the nested method for solving case 2 are shown in Table 7. The fourth and fifth columns list the numbers among the total of 60 runs to reach the global solution in Figure 7 and the suboptimal solution in Figure 8, respectively. If the results are summed, we can see that our nested method has good performance in global optimality. When CONOPT was used as the NLP solver, the performance between GDP and ODP was closely similar. However, with MINOS as the NLP solver, the ODP performed much better than GDP.

1786


Figure 5. Flowchart of the nested method. Table 4. Iterative Process of Case 1 iteration no.

boolean variables (0/1)

objective value

feasibility

0 1 2 3 4 5 6 7 8

00000000000000000011111 00000000000010000011101 00000000000001000011101 00010000000001000011101 00001000010001000011101 00001000000001000011101 00001000000001010011101 00001000000001010011001 00000000000001010011001

95415.23 84189.18 84189.18 87639.28 92573.90 95039.28 92830.18 85431.14 82042.91

yes yes yes yes yes yes yes yes yes

Figure 6. Global optimum solution structure of case 1.

The results again demonstrate the advantage of the proposed ODP for its basic independence on the inner-loop NLP solver since it generates well-conditioned sub-NLPs which can be easily solved by most NLP solvers.

Table 5. Statistical Results of Case 1 no. of sub-NLP global average average major average objects solver model runs NLP count iteration time (s) 23

CONOPT

5. Water Network Synthesis Applications MINOS

The water network synthesis (WNS) problems were tested to illustrate the general application of the proposed nested method and the ODP modeling strategy. All these problems are complicated because of the presence of bilinear constraints which are different from the constraints of the HENS problems which are usually linear except for the cost functions.

GDP ODP GDP ODP

48/60 54/60 45/60 47/60

1227 1186 1287 1215

9/32 8/31 12/35 11/34

121 105 136 113

5.1. Case 3: WNS1. This example21 was originally an NLP problem which only minimized the fresh water consumption. The process data is presented in Table 8. Later, Lee and Grossmann adapted it by including the piping cost into the

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010 Table 6. Data of Case 2

1787

Table 8. Process Unit Data of Case 3 Stream Data

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

Parameters parameter

value

parameter

value

parameter

value

CFi,j Ci,j Ui,j CHU

0 $145.63/y 0.852 kW/m2 K 0

CFcu,i Ci,cu Ui,cu CCU

0 $145.63/y 0.852 kW/m2 K $18.12/kW/y

CFhu,j Cj,hu Uj,hu ∆Tmapp

0 $145.63/y 0.852 kW/m2 K 10 K

unit number

discharge load (kg/h)

maximum inlet concentration (ppm)

maximum outlet concentration (ppm)

1 2 3 4

2 5 30 4

0 50 50 400

100 100 800 800

objective function, after which they reformulated it as a GDP problem.22 The objective function was valued by the fresh water flow rate (t/h) and pipe costs (6.0 for each). In addition, a lower bound of 1 t/h was applied for each existing flow. The superstructure is shown in Figure 9. Lee and Grossmann solved this problem by introducing the tight convex under/over estimators for the bilinear terms, and a two-level branch and bound algorithm was employed. In this project, we formulated the same problem using the ODP and then solved it with the nested method. The ODP model is presented as follows:


Figure 8. Suboptimal solution structure of case 2. Table 7. Statistical Results of Case 2 no. of objects

sub-NLP solver

model

global runs

suboptimal runs

average NLP count

average major iteration

average time (s)

60

CONOPT

GDP ODP GDP ODP

25/60 32/60 3/60 18/60

27/60 27/60 8/60 27/60

35802 35543 51426 37804

29/89 32/92 43/103 40/100

3562 3466 6079 3665

MINOS

1788


(1) The objective function

∑

min Z )

(FUi.F + FUi.PC) +

i∈U,FUi.Y)true

∑

∑

∑

UUi,j.PC +

i∈U j∈U,UUi,j.Y)true

UDi.PC,

i*j

i∈U,UDi.Y)true

(25) (2) Mass balance constraints FiCini )

∑

UUj,i.FCoutj, i ∈ U, i * j

j∈U,UUj,i.Y)true

FiCouti ) FiCini + MLi,

[

(26)

]

(3) Disjunctive constraints for mass balance FUi.Y ) true

∑

Fi ) FUi.F +

[

j∈U,UUj,i.Y)true

[

FUi.Y ) false Fi )

∑

UDi.Y ) true

∑

Fi ) UDi.F +

[

j∈U,UUi,j.Y)true

UDi.Y ) false Fi )

∑

]

]

i, j ∈ U, i * j

]

(27)

UUi,j.F ∨

]

UUi,j.F ,

j∈U,UUi,j.Y)true

[ [ [

UUj,i.F ∨

UUj,i.F ,

j∈U,UUj,i.Y)true

i, j ∈ U, i * j

(4) Object-oriented constraints for each possible unit FUi.Y ) true FUi.F g LB , i ∈ U FUi.PC ) 6 UUi,j.Y ) true UUi,j.F g LB , i, j ∈ U, i * j UUi,j.PC ) 6 UDi.Y ) true UDi.F g LB , i ∈ U UDi.PC ) 6

] ]

FU UU UD Fi

flow rate of a process unit i

Y

TT UD

(29) set of process units set of pipes from source to process units set of pipes between process units set of pipes from process units to discharge

9. We can observe that the ODP have the advantage of less dependency on NLP solvers. 5.2. Case 4: WNS2. The model of case 3 is linear in the objective function. For this model, we tested another example with nonlinear objective functions. It is a water network synthesis problem consisting of three effluent streams with three contaminants and three treatment units. The data are from the papers of Wang23 and Gunaratnam.24 The superstructure is shown in Figure 11. This problem has 21 discrete variables, of which 6 variables denote the existence of discharge pipes for each unit, 9 variables are for the pipes from the 3 process units to 3 treatment units, and 6 variables are for the treatment unit interconnections. The 3 contaminants are hydrogen sulfide, oil, and suspended solids. The process stream data and removal ratio for each contaminant are listed in Table 10. The annual rate of return was set at 10%, and the operating time was 8600 h/y. The environmental limits for the three contaminants were as follows: 2 ppm (H2S,); 2 ppm (oil); and 5 ppm (suspended solids). PT

(28)

(5) Variable bounds Fi g 0, Cinimax g Cini g 0, Coutimax g Couti g 0 FUi.Y, UUi,j.Y, UDi.Y ∈ {true, false}, i, j ∈ U, i * j

U

Figure 9. Superstructure of case 3.

i∈U

set of pipes from process unit i to treatment unit j set of pipes between treatment units set of pipes from units to discharge

i

index of the process units

j

index of the treatment units

Y FR cin(i,c)

EL(c)

u

index of the set, UD

CP

PC

AR

annual rate of return

OP

F

flow rate of a pipe

H

Cin

input concentration into the a process unit output concentration out of a process unit

The initial configuration was obtained by assuming that each process unit had direct access to both the fresh water source and wastewater discharge. Furthermore, it was assumed that the network had no interconnections for any two process units, corresponding to the largest freshwater consumption. As for the parameters of TS, L2 was set to be M/3, and K was M/4. The nested method successfully reached the best known solution, 126.0. The global optimum structure is shown in Figure 10. The statistical results of a total of 100 runs are shown in Table

the input concentration of contaminant c for treatment unit i the output concentration of contaminant c for treatment unit i the environmental limit for contaminant c capital cost parameter operating cost parameter cost of a pipe

cout(i,c)

the existence of a unit cost of a pipe

Cout

the existence of a unit flow rate of a unit

operating hours

C

The ODP formulation for this problem is as follows: The objective function consists of: (a) the capital cost of the treatment units; (b) the operating cost of the treatment units; and (c) the cost of the three types of pipes. (1) Objective function min

∑

(AR)CP(j)FR(j)0.7 +

j∈TU

∑

HOP(j)FR(j) +

j∈TU

∑

∑

PTi,j.Y)true

PTi,j.C +

∑

UDu.C +

UDu.Y)true

TTja,j.C,

TTja,j.Y)true

i ∈ PU, j, ja ∈ TU, j * ja, u ∈ UD

(30)


1789

Figure 10. Global optimum solution structure of case 3. Table 9. Statistical Results of Case 3

Table 10. Data of Case 4

no. of sub-NLP global average average major average objects solver model runs NLP count iteration time (s) 20

CONOPT GDP ODP MINOS GDP ODP

90/100 91/100 63/100 77/100

940 944 934 950

6/12 7/13 6/12 7/13

Treatment Costs

100 97 110 102

units

type of costs

8600

TU1

capital cost ($) operating cost ($/h) capital cost ($) operating cost ($/h) capital cost ($) operating cost ($/h)

16800f 0.7 1.0f 12600f 0.7 0.0067f 4800f 0.7 0

TU2 TU3

The treatment cost was calculated as illustrated in Table 10. In the table, “f ” indicates the flow rate through the treatment pipe is calculated as follows: units. The cost of a pipe costi,j Ai,j ) 3600Fi,j /Vi,j costpipe i,j

) (aAi,j + bYi,j)di,j

(31)

where the flow rate F is in tons per hour, the pipe cross-sectional area A is in square meters, the velocity V is in meters per second, a ) 3603.4, and b ) 124.46. The distance matrix between the two ends is also given in Table 10. It should be noted that the velocity of the flows in each pipe was assumed to be 1 m/s, and the flow rate of each existing stream was assumed to be above a lower bound, LB ) 0.5 t/h. The mass balance equations, which will hold irrespective of the Boolean variables, are as follows: (2) Mass balance eqations FR(i)cout(i, c) ) ML(i, c), FR(j) )

∑

PTi,j.F +

∑

∑

TU1

TU2

TU3

discharge

0 0 0 30 150 50 170

0 0 0 40 125 80 90

0 0 0 35 120 40 100

30 40 35 0 80 30 35

150 125 120 80 0 65 55

50 80 40 30 65 0 40

170 90 100 35 55 40 0

Process Streams process unit

flow rate (t/h)

contaminant

mass load (g/h)

PU1

13.1

PU2

32.7

PU3

56.5

H2S oil SS H2S oil SS H2S oil SS

5109 131 327.5 548706 3597 1308 1412.5 5650 1977.5

Removal Ratio of the Treatment Units treatment unit

RR(H2S)

RR(oil)

RR(SS)

TU1 TU2 TU3

99.9% 90% 0

0 90% 95%

0 97% 20%

j, ja ∈ TU, j * ja (PTi,j.Fcout(i, c)) +

PTi,j.Y)true

(TTja,j.Fcout(ja, c)),

j, ja ∈ TU, c ∈ C, j * ja

TTja,j.Y)true

cout(j, c) ) cin(j, c)(1 - RR(j, c))

∑

PU3

TTja,j.F,

TTja,j.Y)true

∑

PU1 PU2 PU3 TU1 TU2 TU3 discharge

PU2

i ∈ PU, c ∈ C

PTi,j.Y)true

FR(j)cin(j, c) )

Distance Matrix (m) PU1

(UDu.Fcout(u, c)) e (

UDu.Y)true

∑

no. of sub-NLP global average average major average objects solver model runs NLP count iteration time (s)

UDu.F)EL(c),

21

UDu.Y)true

u ∈ U, c ∈ C

CONOPT MINOS

(32)

Figure 11. Superstructure of case 4.

Table 11. Statistical Results of Case 4

GDP ODP GDP ODP

56/100 87/100 4/100 89/100

1298 1386 1348 1385

6/13 7/14 7/14 7/14

155 146 153 144

1790



[

(3) Disjunctive equations UDi.Y ) true

∑

FR(i) ) UDi.F +

[

PTi,j.F ∨

]

PTi,j.Y)true

UDi.Y ) false

[

]

FR(i) )

∑

UDj.Y ) true

FR(j) ) UDj.F +

[

i ∈ PU, j ∈ TU

PTi,j.F ,

PTi,j.Y)true

∑

FR(j) )

∑

]

TTj,ja.F ,

TTj,ja.Y)true

[ [ [ [

]

TTja,ja.F ∨

TTja,ja.Y)true

UDj.Y ) false

that the formulation in the ODP improves the convergence in the global optimality and that the nested method obtained an optimal solution of 482 688.97 for most runs. The structural configuration is shown in Figure 12. Although the authors have modified the problem, information about the best known solution remains unavailable. Nevertheless, we believe that the global solution is the best solution as it has been proven by a subsequent examination to strengthen the TS search ability. The results of 10 runs from different initials all reached the same solution, each with more than 7000 NLP counts and 800 s on the average. 6. Conclusions

j, ja ∈ TU, j * ja

(4) Object-oriented equations UDi.Y ) true UDi.F g LB

(33)

] ] ]

,

Wcost(i) ) (b + aw(i)/V)DM(i, d)/CE

i ∈ PU, d ∈ discharge

UDj.Y ) true UDj.F g LB

,

UDj.C ) (b + aUDj.F/V)DM(j, d)/CE


PTi,j.Y ) true PTi,j.F g LB

(34) Acknowledgment

,

PTi,j.C ) (b + aPTi,j.F/V)DM(i, j)/CE

i ∈ PU, d ∈ discharge TTja,j.Y ) true TTja,j.F g LB

This research was supported by 973 Program of China (No. 2009CB320603), National Science Foundation of China (No. 60704029), and 863 Program of China (No. 2007AA04Z192).

]

Literature Cited

,

TTja,j.C ) (b + aTTja,j.F/V)DM(ja, j)/CE

The object-oriented disjunctive programming (ODP) model has been proposed as an improvement of the traditional GDP model. It can generate concise and well-posed sub-NLPs after the Boolean variables have been specified by removing the equations and variables associating it to the nonexisting units. A nested method has also been proposed for solving ODPs as well as GDPs based on the master and slave decomposition strategies. In the outer-loop, a heuristic algorithm iterates for the structure optimization represented by the Boolean variables. In the inner-loop, a gradient-based NLP solver deals with a subNLP optimization generated after the Boolean variables have been specified. The proposed modeling and solving methods have been successfully applied to process synthesis of HENS and WNS.


(5) Variable bounds cout(i, c), FR(j), cin(j, c), cout(j, c), TUcost(j) g 0 UDi.Y, UDj.Y, PTi,j.Y, TTja,j.Y ∈ {true, false},

(35)

i ∈ PU, j, ja ∈ TU, c ∈ C, j * ja We applied the nested method to solve this problem. The initial structure was randomly generated. While case 4 only has one object more than case 3, the objective function is more complex. The parameters for TS were set almost the same as that in case 3, except that parameter K was increased from 5 to 7 to enlarge the neighborhood. The statistical results of a total of 100 runs are presented in Table 11. We can see

(1) Grossmann, I. E. Review of nonlinear mixed-integer and disjunctive programming techniques. Optim. Eng. 2002, 3, 227–252. (2) Tu¨rkay, M.; Grossmann, I. E. Logic-based MINLP algorithm for the optimal synthesis of process networks. Comput. Chem. Eng. 1996, 20, 959–978. (3) Lee, S.; Grossmann, I. E. New algorithms for nonlinear generalized disjunctive programming. Comput. Chem. Eng. 2000, 24, 2125–2141. (4) Vecchietti, A.; Grossmann, I. E. LOGMIP: a disjunctive 0-1 nonlinear optimizer for process system models. Comput. Chem. Eng. 1997, 21, 427–432. (5) Floudas, C. A.; Akrotirianakis, I. G.; Caratzoulas, S.; Meyer, C. A.; Kallrath, J. Global optimization in the 21st century: advances and challenges. Comput. Chem. Eng. 2005, 29, 1185–1202. (6) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by simulated annealing. Science 1983, 220, 671–680. (7) Glover, F. Taboo Search- Part I. INFORMS J. Comput. 1989, 1 (3), 190–206. (8) Glover, F. Taboo Search- Part II. INFORMS J. Comput. 1990, 2 (1), 4–32.

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010 (9) Goldberg, D. E. Genetic Algorithm in search, optimization, and machine learning; Addison-Wesley: Reading, MA, 1989. (10) Athier, G.; Floquet, P.; Pibouleau, L.; Domenech, S. Process optimization by simulated annealing and NLP procedures application to heat exchanger network synthesis. Comput. Chem. Eng. 1997, 21, 475– 480. (11) Tayal, M. C.; Fu, Y.; Diwekar, U. M. Optimal design of heat exchangers: a genetic algorithm framework. Ind. Eng. Chem. Res. 1999, 38, 456–467. (12) Lin, B.; Miller, D. C. Solving heat exchanger network synthesis problems with Taboo Search. Comput. Chem. Eng. 2004, 28, 1451–1464. (13) Chen, X.; Li, Z. H.; Yang, J.; Shao, Z. J. Nested Taboo Search (TS) and sequential quadratic programming (SQP) method, combined with adaptive model reformulation for heat exchanger network synthesis (HENS). Ind. Eng. Chem. Res. 2008, 47, 2320–2330. (14) Rosenthal, R. E. GAMSsA user’s guide; GAMS Development Corp: Washington, DC, 2008. (15) Yee, T. F.; Grossmann, I. E. Simultaneous optimization models for hear integration---I. Area and energy targeting and modeling of multistream exchangers. Comput. Chem. Eng. 1990, 14 (10), 1151–1164. (16) Yee, T. F.; Grossmann, I. E. Simultaneous optimization models for hear integration---II. Heat exchanger network synthesis. Comput. Chem. Eng. 199014, (10), 1165–1184. (17) Ciric, A. R.; Floudas, C. A. Heat exchanger network synthesis without decomposition. Comput. Chem. Eng. 1991, 15 (6), 385–396.

1791

(18) Lee, S.; Grossmann, I. E. A global optimization algorithm for nonconvex generalized disjunctive programming and applications to process systems. Comput. Chem. Eng. 2001, 25, 1675–1697. (19) 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. (20) Linnhoff, B.; Flower, J. R. Synthesis of heat exchanger networks: Systematic generation of energy optimal networks. AIChE J. 1978, 24 (4), 633–654. (21) Huang, C. H.; Chang, C. T.; Ling, H. C. A mathematical programming model for water usage and treatment network design. Ind. Eng. Chem. Res. 1999, 38, 2666–2679. (22) Lee, S.; Grossmann, I. E. Global optimization of nonlinear generalized disjunctive programming with bilinear constraints: applications to process networks. Comput. Chem. Eng. 2003, 27, 1557–1575. (23) Wang, Y. P.; Smith, R. Design of distributed effluent treatment systems. Chem. Eng. Sci. 1994, 49 (18), 3127–3145. (24) Gunaratnam, M.; Alva-Argez, A.; Kobossis, A.; Kim, J.-K.; Smith, R. Automated design of total water systems. Ind. Eng. Chem. Res. 2005, 44 (3), 588–599.

ReceiVed for reView June 23, 2009 ReVised manuscript receiVed September 21, 2009 Accepted October 15, 2009 IE901010A

Object-Oriented Disjunctive Programming with a Nested Heuristic and

Recommend Documents