- Email: [email protected]
On the use of a mixed integer non-linear programming model for refrigerant design
Pergamon
Int. Trans. Opl Res. Vol. 4, No. 1, pp.45-54, 1997 © 1997 IFORS. Published by Elsevier Science Ltd All rights reserved. Printed in Great Britain PII: S0969--6016(97)00027-5 0969-6016/97 $17.00 + 0.00
On the Use of a Mixed Integer Non-linear Programming Model for Refrigerant Design NACHIKET CHURl and LUKE E.K. ACHENIE University of Connecticut, USA In this paper, a novel mixed integer non-linear programming model for single component refrigerant design is presented. At the heart of the approach is a new formulation for structural feasibility that allows multiple bonds, connectivity and isomers. The strategy defines a set of structural groups (consisting of atoms), subsets of which are combined to form refrigerant molecules. Molecules formed this way must obey structural and stability constraints. The design objective is to build a refrigerant molecule that has desired physical properties and performance characteristics. These attributes are formulated as mathematical programming constraints and performance objectives which involve both continuous and integer variables. With the current renewed interest in the environment, the suggested approach is applied to refrigerant design with an environmental constraint. The results indicate the viability and the flexibility of the approach. © 1997. Published by IFORS/Elsevier Science Ltd.
Key words: Refrigerant, computer aided molecular design, structural feasibility, mixed integer non-linear programming.
1. INTRODUCTION The search for new compounds within the chemical industry is in general a very prolonged and expensive process due to the large number of costly experimental studies needed to identify compounds with desired physical properties and performance attributes. In these cases computer aided design (Gani and Brignole, 1983; Brignole et al., 1986; Joback and Stephanopoulos, 1989; Gani et al., 1991; Odele and Machietto, 1993) has proven to be a viable supplement to the experimental approach. In this paper, a systematic and flexible mixed integer non-linear programming based approach to computer aided molecular design is presented. The paper starts with a description of a general strategy for refrigerant design, followed by the mathematical programming model, physical property constraints, a case study and finally discussions.
2. GENERAL STRATEGY FOR SINGLE C O M P O N E N T REFRIGERANT DESIGN A general strategy for forming a single component refrigerant is as follows: 1. Choose an initial set of structural groups that can potentially form single component refrigerants. 2. Select physical properties which must be set at prespecified levels, such as compatibility with polymers and lubricants, volatility, flammability or toxicity level. 3. Select a performance index, p, upon which refrigerants are evaluated (e.g. ozone depletion potential, coefficient of performance, corrosive tendency, economics, safety, pressure ratio, discharge temperature, capacity). 4. Do the refrigerant design. We achieve this by formulating and solving a mathematical program. 5. Design for best single component refrigerant, second best refrigerant, etc., which may be used in an experimental evaluation. 6. Experimentally synthesize single component refrigerant (if not known to exist) and verify the predicted performance characteristics. In this paper only a proof-of-concept study via computer simulation is presented.
Correspondence: Luke E. K. Achenie, Department of Chemical Entlineeriny, U-222, University of Connecticut, 191 Auditorium Road, Storrs, CT 06269, USA
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programmin9 Model
46
3. M A T H E M A T I C A L P R O G R A M M I N G M O D E L A general mathematical programming model for refrigerant design with a performance objective p(x,p) can be posed as max p(x,~')
(1)
x,f
subject to CS' ~< d
(2)
h(~) ~< 0
(3)
g(x) ~< 0
(4)
where ~' is a vector of all the structurally induced binary variables and x is continuous (process and design variables). C is a matrix of appropriate dimensions while d, h and g are vectors of appropriate dimensions. Equations of type equation (4) are not used in the case study. The compressor displacement CD is a good measure of the performance of a refrigerant in the cycle, and the ozone depletion potential ODP accounts for the environmental damage caused by certain compounds by depleting the Earth's ozone layer. Therefore, for a single component refrigerant design problem, an appropriate performance measure p is a function of CD, ODP and selected physical properties P as follows:
p = p(CD(~),ODP(~),P(~)).
(5)
3.1. Structural feasibility constraints For this design, a basis set of structural groups is chosen. Structural groups consist of connected atoms formed by a combination of functional groups, and whose net valence is at least 1. Examples of structural groups by this definition include C H 3 - - , CH2 F - and C H C I = . The selection of the basis set depends on the intended application and availability of accurate group contribution techniques for predicting the properties of interest. Among the attributes that the refrigerant should have is that it should not polymerize over a long period of time. This requirement is met by excluding double or triple bonds from the molecule. In addition, most refrigerants used for automotive air conditioning are aliphatic compounds. Based on these and other considerations, a representative set of structural groups can be identified: [-CH3--, C H 2 = , - - C H = , = C = , X], where X is a halogen atom. This set is by no means exhaustive. 3.1.1. Variables. To help with the exposition, we first define the following parameters: m vk Smax n nmax
= = = = =
number of structural groups in the basis set valence of the kth group in the basis set maximum valence of all the groups in the basis set number of groups in a molecule maximum number of structural groups in a molecule
We then specify a basis set of m structural groups having valencies {Vk}, with maximal valence SmaX. The maximum number of groups allowed in a molecule is nm,x. The actual number of groups, n, is obtained from the solution of the mathematical programming model. In addition, let us define the following binary variables:
Ulk =
Zijp =
{~ if the ith group in the molecule is of the kth kind otherwise
{~ if the ith group's jth site is attached to the pth group otherwise.
Finally, we define w to be a binary vector of length nm~x such that its first n terms are zero, while the
International Transactions in Operational Research Vol. 4, No. I Cl 4
F3
m3if~ 2 ",--2 1
tCll
Basis Set
231!H 7
C15
47
F6
21
CH2 CH3
543
FcCH
6
C1
Non-zero binary variables z w 112 317 713 8 722 9 217 412 736 : 224 512 nmax 231 245 617
u 16 24 35 46 56 65 73
Fig. 1. Example molecule.
remaining terms are one. This vector is introduced in the formulation to ensure that a connected molecule is formed from the structural groups. Figure 1 explains how these variables and parameters are to be interpreted for 1,1,l-trichloro-2,2difluoroethane. The bold numbers next to the groups are the g r o u p numbers in the molecule, and the numbers next to the bonds are the site numbers for the individual groups. A basis set of six groups (m = 6), namely CH3, CH2, CH, C, F and CI is chosen for illustration and only the non-zero terms of the variables are given. F o r example, the non-zero terms of u are u16, u24, u 35; a non-zero u 16 implies that the 1st g r o u p in the molecule is of type 6. F r o m the basis set we can see that the 6th g r o u p is C1. Similarly, u24 implies that the 2nd g r o u p is a C. We do not distinguish between the various ways in which C can be connected. Distinction between the various types of bonding ( > C < , - - C ~ = , = C = and > C = ) is d o n e within the property prediction procedures since bonding information is available from the z variables. The z variables are interpreted in a m a n n e r similar to u. A non-zero z 713 implies that the 7th group's 1st site is attached to the 3rd group. Since 1273and u35 are non-zero we know that the 7th g r o u p is of type 3 (CH) and the 3rd g r o u p is of type 5 (F). Thus a non-zero z713 indicates a b o n d between C H and F. Since there are a total of seven groups in the molecule, n = 7, and the first seven terms of w are zero and the remaining terms (w8... w.... ) are non-zero. In this case, n . . . . which is specified a priori, is at least 7. In summary, the parameters for this example are: m = 6; v = [1 2 3 4 1 1]; Smax = 4 = max{vk}; and n = 7. The purpose of the structural feasibility constraints is to generate molecules that do not violate basic feasibility criteria such as the Octet Rule. The molecules should not have any unattached sites or multiple bonds attached to the same site. O u r basic philosophy here is to develop a model that is linear in the integer variables. Using the variables defined above, these criteria are expressed as constraints of type equation (2): nmax
Z ~ /-'/ik~ nmax i=lk=l nmax Sma~
(6)
~.
i = I . . . nmax p=lj=l
(7)
k=l
i = 1 Smax
Z
Z
>1 - w,
i = 2 . . . nmax
(8)
p=lj=l nmax nl
nma x
2 Zu,,+ 2w,=nm.x
i=lk=l
W1 = 0 w i <~ w i + 1
(9)
i=1
i = 1 ...(nma x -- 1)
(lO) (11)
48
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programming Model Smax nmax
nmax
~ zijp- ~ z,v,p+M(U,k-- 1)~<0 j = vk p = I Smax
....
k= 1...m
(12)
Srnax
2 Zijp = 2 Zpji j=l"
i= l...n
p= I
1 . . . ( n m a x - - 1),p = (i + 1 ) . . . r i g a x
i:
(13)
j=l
~ x Zijp ~ 1
i = 1 ... nmaxJ : 1 ... Smax
(14)
i=
(15)
p=l
Uik - ~ Ui-1, k 4 0 k=l
2 . . . n m a x.
k=l
Equation (6) puts a limit on the number of structural groups that can constitute a molecule. The double-summation on the left-hand side is equal to n, the actual number of groups present in the molecule. Since a minimum of two groups are required to form a molecule and an upper limit of nmax is specified, 2 ~< n ~< nmax. Equation (7) is an implementation of the Octet rule. The expression on the left gives the number of bonds attached to the ith group, while that on the right states that if the ith group is of the kth type, its valence is vk. This ensures that the number of bonds attached to a group equals the valence of the group. Equations (8), (9), (11) and (10) ensure that only one molecule is formed. Due to the presence of equation (9), equation (6) becomes redundant. This is realized by constraining the ith group to be attached to one of the groups before it, that is, groups 1 to (i - 1) [equation (8)]. Thus the second group has to be attached to the first group, the third group has to be attached to either of the first two groups, and so on. The first group is always present [equation (10)], and the (i + 1)th group is present only if the ith group is present [equation (11)]. Equation (12) has to be introduced to account for different valencies of the groups. This equation is a linear analog of the non-linear equation Srnax
nmax
Uik ~ ~, zijp = 0 J=vk+lp=l
i = 1...n .... k = 1...m
(16)
which states that if the ith group is of the kth kind, then sites (vk + 1) to Smax of that group should not have any connections since they are non-existent. Equation (12) [which is equivalent to Z~=,tlZ~=,,~Zijp + M(U~k -- 1) ~< 0] is written in a form that is convenient from a computer programming point of view. In order to tighten the relaxation, M should be chosen to be the smallest positive number that is at least as large as the maximum attainable value of the first two terms in the equation. In our simulations we simply chose M large enough. Equation (13) is the symmetry constraint. Hence if a group is attached to a second group, the second group is automatically attached to the first one. Equation (14) ensures that a group's site can be attached at most once to some other group. Equation (15) is applied to force existence of the ith group if the (i - 1)th group is present. The structural feasibility constraints [equations (6)-(15)] are linear; hence they form a convex hull separating feasible molecular structures from infeasible ones. One property constraint that has to be incorporated is thermodynamic feasibility of type equation (3). A molecule is thermodynamically feasible at a particular temperature T if its Gibbs' Free Energy of formation Gf a t that temperature is negative, that is Gf(T) < 0.
(17)
For simplicity, Gf is calculated at standard temperature and pressure. While this constraint ensures stability for the generated molecule, it can also eliminate some potentially useful molecules. This may happen if a thermodynamically unstable molecule's dissociation kinetics are such that it can exist for time periods that are long enough for commercial use.
3.2. Characteristics of structural feasibility model The variables used to describe the molecule give nearly complete information about the molecule's composition and connectivity. One of the advantages of this new formulation becomes apparent when one looks at the various group contribution methods. A number of these methods make use of a
International Transactions in Operational Research Vol. 4, N o . I
Cyclic Molecule (one l,oop)
CHIN, / CH 2 CH ~
49
F
z!H
C~ Cyclic Molecule (two independent loops)
Unconnected Molecules (two sub-molecules)
.CH2
H/
CI
CH 2 F + CH 3 F
Fig. 2. Loops and sub-molecules.
certain amount of structural information in giving accurate predictions. For example, the group contribution method for liquid specific heat (Cheuh and Swanson, 1973) requires the addition of 4.5 to the value of Cp~ for any carbon group which fulfils the following criterion: 'A carbon group which is joined by a single bond to a carbon group which is connected by a double or triple bond with a third carbon atom.' It is obvious that this correction cannot be added without any connectivity information. A search of the open literature gave no previous C A M D formulation that was able to give the level of detailed molecular information required for the more accurate and more complex group contribution techniques. Another way in which structural information proves useful is in incorporating bonding constraints. One of the desired properties of a refrigerant is its stability, that is, the compound should not spontaneously decompose or polymerize. In many cases, polymerization can be eliminated by ensuring that the refrigerant does not have double or triple bonds. There is no systematic way of introducing this constraint mathematically in the absence of detailed structural information. With the new formulation, the constraint that has to be added to eliminate multiple bonds is: $
x
zij p <~ 1
i = 1 ...(nma ~ - 1),p = (i + 1)...nma x.
(18)
j=l
From graph theory (Busacker and Saaty, 1965) 'a connected planar graph with n vertices, a edges, and l regions (including the outside or unbounded region)' is defined by the Euler formula l = a - n + 2. We can use a variation of this relation to suppress the formation of a cyclic molecule if we so desire using the equation: l = a - n + b.
(19)
In this equation, l is the total number of independent loops, b is the number of 'sub-molecules' being designed, a is the total number of attachments, and n is the total number of groups. Obviously, b = 1 since we want a single connected molecule. A loop is a series of distinct groups that are attached in such a manner that one can go from one group to others, and return to the first group via a different route. A loop is independent if it cannot be described by a combination of other loops. Two groups are said to be attached if there exists at least one bond between them. These terms are explained with examples in Fig. 2. The total number of attachments (a) and the number of groups (n) are given by a =
Cip i= 1
p=i+ 1
i=Ik=l
(20)
50
N. Churi and L.E.K. A c h e n i e - - M i x e d Integer Non-linear Programmin9 Model OH
OH
OH
OH
N o n - z e r o t e r m s o f 'z' 112 122 211 221
112 122 211 221
Fig. 3. I n d i s t i n g u i s h a b l e molecules.
B a s i s Set 1
CH 3
2
CH 2
3
F
• N o n - z e r o t e r m s o f ' u ' a n d 'z' H
i1 22 33
Z
U
112 211 223 312
12 21 33
Z
113 122 211 311
Fig. 4. Multiplicity.
where {~ Ci p =
i f ESm'xz.. j = 1 13P
= 0
i = 1...(nm, x - 1 ) , p = ( i + 1)...nma x otherwise.
By using these relations in equation (19) and setting l = 0, the resulting molecule will not have any loops and will thus be acyclic. In summary, the formulation gives us a great deal of control over the features of the target molecules. We can distinguish between isomers, allow multiple bonds, and also specify whether or not we desire a cyclic molecule• One small limitation to the connectivity information is that one cannot distinguish between the individual bonds in the case of multiple bonds. For example, the two structures in Fig. 3 cannot be distinguished under the formulation. However, this is not a serious drawback since such cases are not common. In any case, no existing group contribution method is able to distinguish between such structures in a systematic manner. The formulation also does not consider any constraints based on bond angles. While it is possible to incorporate bond angles, we have not found any property prediction methods that can take advantage of this additional information. Figure 4 shows two possible ways in which the variables can describe a molecule. In general, there are several such multiplicities in specifying a molecule. This can lead to increased computational expense. The degree of multiplicity depends on the number of groups present in the molecule and also their valencies. Sine these different solutions describe the same molecule, they have identical objective function values. This fact can be used to eliminate multiplicity by applying a constraint to the objective function value of the type Pne,,, <<-Pola -- ~" However, such a constraint can cause problems with the termination criterion which depends on an increase in the objective value to detect the minima. Also, since other structures may have the same objective function value, there is a risk of
International Transactions in Operational Research Vol. 4, No. I
51
eliminating potentially good molecules. A better approach is to eliminate multiplicity by introducing an integer cut every time a new solution has the same objective function value as that obtained in the previous iteration. The integer cut given below eliminates $ from the search space: ~ - ~ 9~ ~< I BI - 1 iGB
(22)
iEN
where B = {il .~ = I/, N = {i[ ~i = 0 / and ] B ] is the cardinality of B. For the kth iteration, there will be k - 1 such constraints, each eliminating one of the previous solutions from the search space. A more rigorous approach is to introduce constraints that force the sites to be ordered. This will eliminate multiplicity by ensuring that there is only one way in which the variables can describe a molecule.
4. PHYSICAL P R O P E R T Y CONSTRAINTS For refrigerant design, the important refrigerant physical properties and process variables include (i) boiling point, (ii) flammability, (iii) materials and lubricants compatibility, (iv) stability, (v) toxicity, (vi) miscibility, (vii) the saturated liquid (or vapor) mixture pressure in the evaporator (or condenser), (viii) the saturated vapor temperature in the condenser, (ix) the compressor displacement, and (x) the ozone depletion potential. However, we have not used all of these in the case study. Constraints on physical property values are expressed as bounds and are of type equation (3). Thus, suppose one wants the refrigerant properties P1, P2--Pm to lie within a given interval. For example, Pt, P2 and P3 may correspond to the flammability, ODP and boiling point, respectively. This is expressed as: pmin <~Pi <~ Pmax,~/i = 1...m.
(23)
Note that Pi is implicitly a function of ~ which defines the molecular structure. The above equations mean that a property, Pi should lie between the lower and upper bounds p~in and pmax respectively. This allows us to accept refrigerants which have a boiling point within a given range, for instance. In the case study considered, constraints on the evaporator and the condenser pressures (Joback and Stephanopoulos, 1989) have been included. The lowest pressure in the cycle should be greater than the atmospheric pressure; this reduces the possibility of air and moisture leaking into the system. A high system pressure increases the size, weight and cost of equipment (Dossat, 1981). A pressure ratio of 10 is considered to be the maximum for a refrigeration cycle (Perry et al., 1984). Consistent with this, we choose the system pressure Pvp(Te) at the evaporating temperature and the system pressure Pvp(Tcd) at the condensing temperature to satisfy Pvp(Te)/> 1.4 bar
(24)
Pvp(Tcd) <~ 14 bar
(25)
where T~ and Ted are functions of the molecular structure 9A larger enthalpy of vaporization reduces the amount of refrigerant required. The heat of vaporization of Freon 12 is 18.40kJ/mol (Ashrac, 1972). By giving a lower limit of 15.0kJ/mol, the generated molecule will not be significantly worse than Freon 12 in terms of H~. This relaxation is not done to account for errors in the group contribution method, but to prevent elimination of potentially good molecules whose Hv is only slightly lower than that of Freon 2: Hv(T~) >1 15.0 kJ/mol.
(26)
A low liquid specific heat reduces the amount of refrigerant that flashes upon passage through the expansion valve. This is evaluated at T~,a,. The value of the specific heat for Freon 12 obtained from a group contribution method at this temperature is 27.12 cal/mol. As above, the limit is a relaxation of this value: C p ~ ( T ~ ) ~ 30.0 cal/mol.
(27)
By the Clean Air Act of 1995, refrigerants should have ODPs below 0.2 by the year 2000. Hence this value is used as the upper limit for ODP. The ODP of Freon 12 is 1.0.
52
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programming Model ODP <~0.2.
(28)
For a more practical design, materials and lubricants compatibility, toxicity and flammability ought to be incorporated into the mathematical programming model. The constraints discussed so far are by no means exhaustive, and several different ones can be added to achieve a specific refrigerant design objective.
5. CASE STUDY
5.1. Process conditions The following conditions are used for the refrigeration cycle (Duvedi and Achenie, 1996): (i) evaporating temperature Te = - 1.1 °C, (ii) condensing temperature Ted = 43.3 °C, (iii) average temperature Tmea, = 21.1 °C, and (iv) saturated conditions are assumed so that superheat temperature Ts equals T~. The relevant physical properties are hv(T~), Cp,(Tm~,,), Pvp(T¢) and Pvp(Tca). Selection of parameter values depends on the application to which the designed compounds are to be applied. For the case study, the aim is to identify hydrochlorofluorocarbons (HCFCs) that have properties comparable to, or better than, Freon 12, but with significantly lower ODP. The basis set for this case study consists ofCH3, CH2, CH, F and CI. This gives m = 5, Smax = 3, and v = [1 2 3 1 1]. The maximum number of groups in the molecule (nmax)is specified as eight, which corresponds to the highest substitution possible for propane from the groups in the basis set. This limit comes from the observation that higher molecular weight molecules do not have vapor pressure in the range desired for refrigerants. A particular form of the performance objective in equation (5) used in this case study is
f = Cp,(T,ne~,)/nv(Te)
(29)
which is to be minimized. A larger enthalpy of vaporization reduces the amount of refrigerant required. A low liquid specific heat reduces the amount of refrigerant that flashes upon passage through the expansion valve. Pressure ratio, discharge temperature, capacity and COP sensitivity to cycle variability also need to be reflected in the performance criteria for a more realistic cycle. However, in this proof-of-concept study, we will use a simpler model. Notice that equation (3) is non-linear in the integer variables. In the case study considered, t,he physical property constraints are of this form. By simply adding constraints of type ~ = ~ where ~ is continuous, equation (3) becomes non-linear in ~ instead of ~. Also, the objective function becomes a non-linear function of ~. Therefore the sub-problem is a non-linear program for this case study.
5.2. Solution methodology The locally optimal augmented penalty outer approximation (AP/OA) algorithm (Viswanathan and Grossmann, 1990) is employed to solve the MINLP. Initially, a relaxation of the M I N L P is solved. Subsequently, a finite sequence of non-linear programming (NLP) and mixed integer linear programming (MILP) problems are solved. In the N L P sub-problems, derivative information is obtained by perturbation. The flowchart for the algorithm is given in Fig. 5.
5.3. Results and discussions Table 1 shows the results obtained for the problem described above. The values of the property constraints for the molecules formed are also given. The first row gives the locally optimal solution to the relaxed NLP. However, it has no physical significance since the solution variables are not integers. The first feasible molecule formed is CH3--CI with an objective function value of 9.8513. The next iteration gives F - - H C = C H - - F with an objective function value of0.8911. Since there is an increase in the objective value, iterations stop at this point with CH3--C1 as the locally optimal molecule.
International Transactions in Operational Research Vol. 4, No. 1 Initial Guess
53
]
I ] RelaxedNLP I Yes
I
I 'Su ro lem I
Fig. 5. AP/OA algorithm: minimization problem.
Once a locally optimal solution is found, the next best solution is found by introducing an integer cut in the formulation. In this manner, it is possible to get a series of solutions in decreasing order of optimality as given in Table 2. Since a globally optimal solution is not obtained for the NLP sub-problem, it is conceivable that the sequence of local minima may not be monotonic. We did not however encounter this problem in the case study. As can be seen from Table 2, not all molecules may be suitable for practical purposes. Molecules like F - - H C = C H - - F may polymerize over the long run and hence cannot be used as refrigerants. It is possible to eliminate such molecules by introducing bonding constraints [equation (18)] into the formulation. In any case, the formulation can give a series of molecules that conform to all applied constraints and which can be transferred to the experimental stage for further evaluation. Due to the non-convexity of the formulation, the global optimum cannot be guaranteed when using the AP/OA algorithm for solving the MINLP model. However, current advances in global optimization [see, for example (Floudas and Visweswaran, 1993)] can be employed to obtain globally optimal solutions to the MINLP formulation. As the problem size increases, there will be a need to consider decomposition and parallel branch and bound algorithms (Pekny and Miller, 1992). The formulation for the single component refrigerant is somewhat computationally intensive. For the case study considered, there are 240 binary variables, 124 linear constraints and 6 non-linear constraints. On an IBM RISC 320H workstation the computation time was of the order of 10-30s for the initial relaxed problem, and 10-15min on the subsequent MILP/NLP sequence. This is mainly due to the bond multiplicity alluded to earlier. To alleviate this problem, an integer cut was Table 1. Single component refrigerant design Objective
h,(T,)
P,p(T~)
Itn.
Molecule
Cp,/h,
(kJ/mol)
Cp,(T,,,,,) (cal/mol)
(bar)
P,p(Tcd) (bar)
OOP
0 1 2
-CH3--CI F--HC=CH--F
0.5726 0.8513 0.8911
19.23 20.44 20.20
11.01 17.40 18.00
3.77 2.44 2.37
14.00 9.17 9.92
0.01 0.00
Table 2. Single component refrigerant design: series of feasible molecules Objective No.
Molecule
Cp,/h,
1 2 3 4
CH3--CI F--HC=CH--F CH3--CH2--F F~H(CH2) 2
0.8513 0.8911 1.1282 1.2458
h,( T,)
C p,(Tmean)
(kJ/mol)
(cal/mol)
Pvo(Te) (bar)
Pvp(Tcd) (bar)
ODP
20.44 20.20 17.78 18.88
17.40 18.00 20.06 23.52
2.44 2.37 4.19 3.51
9.17 9.92 13.07 13.72
0.01 0.00 0.00 0.00
54
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programming Model
introduced every time bond multiplicity was detected. While it is possible to eliminate bond multiplicity, the additional constraints are either non-linear, or if linear, result in an increase in the problem size. At any rate, we are formulating integer constraints to deal with this problem in a rigorous way.
6. CONCLUSIONS The model proposed for single component refrigerant design (with environmental considerations being incorporated) has been shown to be effective. The proposed structural feasibility constraints for the single component refrigerant design give nearly complete molecular information; thus we can easily accommodate any group contribution based property prediction techniques that require detailed molecular bonding information. Sensitivity of the analysis to uncertainties in the physical property estimation is very important but falls outside the scope of this paper. These will need to be addressed in the future through sensitivity analysis.
REFERENCES American Society of Heating, Refrigeration and Air-Conditioning Engineers, Inc. (1972). ASHRAE Handbook of Fundamentals. Brignole, E.A., Bottini, S. & Gani, R. (1986). A strategy for the design of solvents for separation processes. Fluid Phase Equilibria 29 125-132. Busacker, R.G. & Saaty, T. L. (1965). Finite Graphs and Networks: An Introduction with Applications. McGraw-Hill, New York. Chueh, C.F. & Swanson, A.C. (1973). Estimation of liquid heat capacity. Canadian Journal of Chemical Engineering 5, 596-600. Dossat, R.J. (1981). Principles of Refrigeration. Wiley, New York. Duvedi, A. & Achenie, L. E. K. (1996). Design environmentally safe refrigerants using mathematical programming. Chemical Engineering Science 51, 3727-3739. Floudas, C.A. & Visweswaran, V. (1993). A primal-relaxed dual global optimization approach. J. Optim. Theory and Appl. 78(2), 187. Gani, R. & Brignole, E.A. Molecular design of solvents for liquid extraction based on UNIFAC. Fluid Phase Equilibria 13, 331-340. Gani, R., Nielsen, B. & Fredenslund, Aa. (1991). A group contribution approach to computer-aided molecular design. AICHE J. 37, 1318-1331. Joback, K. G. & Stephanopoulos, G. (1989). Designing molecules possessing desired physical property values. In Proceedings FOCAPD 1989, Snowmass Village, Colorado, FOCAPD. Odele, O. & Macchietto, S. (1993). Computer aided molecular design: a novel method for optimal solvent selection. Fluid Phase Equilibria 47-54. Pekny, J. F. & Miller, D. L. (1992). A parallel branch and bound algorithm for solving large asymmetric travelling salesman problem. Mathematical Programming 55, 17-33. Perry, R.H., Green, D.W. & Maloney, J. O. (1984). Perry's Chemical Engineers' Handbook, sixth edn. McGraw-Hill, New York. Viswanathan, J. & Grossmann, I.E. (1990). A combined penalty function and outer-approximation method for MINLP optimization. Comp. & ChE. 14, 769-782.
Int. Trans. Opl Res. Vol. 4, No. 1, pp.45-54, 1997 © 1997 IFORS. Published by Elsevier Science Ltd All rights reserved. Printed in Great Britain PII: S0969--6016(97)00027-5 0969-6016/97 $17.00 + 0.00
On the Use of a Mixed Integer Non-linear Programming Model for Refrigerant Design NACHIKET CHURl and LUKE E.K. ACHENIE University of Connecticut, USA In this paper, a novel mixed integer non-linear programming model for single component refrigerant design is presented. At the heart of the approach is a new formulation for structural feasibility that allows multiple bonds, connectivity and isomers. The strategy defines a set of structural groups (consisting of atoms), subsets of which are combined to form refrigerant molecules. Molecules formed this way must obey structural and stability constraints. The design objective is to build a refrigerant molecule that has desired physical properties and performance characteristics. These attributes are formulated as mathematical programming constraints and performance objectives which involve both continuous and integer variables. With the current renewed interest in the environment, the suggested approach is applied to refrigerant design with an environmental constraint. The results indicate the viability and the flexibility of the approach. © 1997. Published by IFORS/Elsevier Science Ltd.
Key words: Refrigerant, computer aided molecular design, structural feasibility, mixed integer non-linear programming.
1. INTRODUCTION The search for new compounds within the chemical industry is in general a very prolonged and expensive process due to the large number of costly experimental studies needed to identify compounds with desired physical properties and performance attributes. In these cases computer aided design (Gani and Brignole, 1983; Brignole et al., 1986; Joback and Stephanopoulos, 1989; Gani et al., 1991; Odele and Machietto, 1993) has proven to be a viable supplement to the experimental approach. In this paper, a systematic and flexible mixed integer non-linear programming based approach to computer aided molecular design is presented. The paper starts with a description of a general strategy for refrigerant design, followed by the mathematical programming model, physical property constraints, a case study and finally discussions.
2. GENERAL STRATEGY FOR SINGLE C O M P O N E N T REFRIGERANT DESIGN A general strategy for forming a single component refrigerant is as follows: 1. Choose an initial set of structural groups that can potentially form single component refrigerants. 2. Select physical properties which must be set at prespecified levels, such as compatibility with polymers and lubricants, volatility, flammability or toxicity level. 3. Select a performance index, p, upon which refrigerants are evaluated (e.g. ozone depletion potential, coefficient of performance, corrosive tendency, economics, safety, pressure ratio, discharge temperature, capacity). 4. Do the refrigerant design. We achieve this by formulating and solving a mathematical program. 5. Design for best single component refrigerant, second best refrigerant, etc., which may be used in an experimental evaluation. 6. Experimentally synthesize single component refrigerant (if not known to exist) and verify the predicted performance characteristics. In this paper only a proof-of-concept study via computer simulation is presented.
Correspondence: Luke E. K. Achenie, Department of Chemical Entlineeriny, U-222, University of Connecticut, 191 Auditorium Road, Storrs, CT 06269, USA
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programmin9 Model
46
3. M A T H E M A T I C A L P R O G R A M M I N G M O D E L A general mathematical programming model for refrigerant design with a performance objective p(x,p) can be posed as max p(x,~')
(1)
x,f
subject to CS' ~< d
(2)
h(~) ~< 0
(3)
g(x) ~< 0
(4)
where ~' is a vector of all the structurally induced binary variables and x is continuous (process and design variables). C is a matrix of appropriate dimensions while d, h and g are vectors of appropriate dimensions. Equations of type equation (4) are not used in the case study. The compressor displacement CD is a good measure of the performance of a refrigerant in the cycle, and the ozone depletion potential ODP accounts for the environmental damage caused by certain compounds by depleting the Earth's ozone layer. Therefore, for a single component refrigerant design problem, an appropriate performance measure p is a function of CD, ODP and selected physical properties P as follows:
p = p(CD(~),ODP(~),P(~)).
(5)
3.1. Structural feasibility constraints For this design, a basis set of structural groups is chosen. Structural groups consist of connected atoms formed by a combination of functional groups, and whose net valence is at least 1. Examples of structural groups by this definition include C H 3 - - , CH2 F - and C H C I = . The selection of the basis set depends on the intended application and availability of accurate group contribution techniques for predicting the properties of interest. Among the attributes that the refrigerant should have is that it should not polymerize over a long period of time. This requirement is met by excluding double or triple bonds from the molecule. In addition, most refrigerants used for automotive air conditioning are aliphatic compounds. Based on these and other considerations, a representative set of structural groups can be identified: [-CH3--, C H 2 = , - - C H = , = C = , X], where X is a halogen atom. This set is by no means exhaustive. 3.1.1. Variables. To help with the exposition, we first define the following parameters: m vk Smax n nmax
= = = = =
number of structural groups in the basis set valence of the kth group in the basis set maximum valence of all the groups in the basis set number of groups in a molecule maximum number of structural groups in a molecule
We then specify a basis set of m structural groups having valencies {Vk}, with maximal valence SmaX. The maximum number of groups allowed in a molecule is nm,x. The actual number of groups, n, is obtained from the solution of the mathematical programming model. In addition, let us define the following binary variables:
Ulk =
Zijp =
{~ if the ith group in the molecule is of the kth kind otherwise
{~ if the ith group's jth site is attached to the pth group otherwise.
Finally, we define w to be a binary vector of length nm~x such that its first n terms are zero, while the
International Transactions in Operational Research Vol. 4, No. I Cl 4
F3
m3if~ 2 ",--2 1
tCll
Basis Set
231!H 7
C15
47
F6
21
CH2 CH3
543
FcCH
6
C1
Non-zero binary variables z w 112 317 713 8 722 9 217 412 736 : 224 512 nmax 231 245 617
u 16 24 35 46 56 65 73
Fig. 1. Example molecule.
remaining terms are one. This vector is introduced in the formulation to ensure that a connected molecule is formed from the structural groups. Figure 1 explains how these variables and parameters are to be interpreted for 1,1,l-trichloro-2,2difluoroethane. The bold numbers next to the groups are the g r o u p numbers in the molecule, and the numbers next to the bonds are the site numbers for the individual groups. A basis set of six groups (m = 6), namely CH3, CH2, CH, C, F and CI is chosen for illustration and only the non-zero terms of the variables are given. F o r example, the non-zero terms of u are u16, u24, u 35; a non-zero u 16 implies that the 1st g r o u p in the molecule is of type 6. F r o m the basis set we can see that the 6th g r o u p is C1. Similarly, u24 implies that the 2nd g r o u p is a C. We do not distinguish between the various ways in which C can be connected. Distinction between the various types of bonding ( > C < , - - C ~ = , = C = and > C = ) is d o n e within the property prediction procedures since bonding information is available from the z variables. The z variables are interpreted in a m a n n e r similar to u. A non-zero z 713 implies that the 7th group's 1st site is attached to the 3rd group. Since 1273and u35 are non-zero we know that the 7th g r o u p is of type 3 (CH) and the 3rd g r o u p is of type 5 (F). Thus a non-zero z713 indicates a b o n d between C H and F. Since there are a total of seven groups in the molecule, n = 7, and the first seven terms of w are zero and the remaining terms (w8... w.... ) are non-zero. In this case, n . . . . which is specified a priori, is at least 7. In summary, the parameters for this example are: m = 6; v = [1 2 3 4 1 1]; Smax = 4 = max{vk}; and n = 7. The purpose of the structural feasibility constraints is to generate molecules that do not violate basic feasibility criteria such as the Octet Rule. The molecules should not have any unattached sites or multiple bonds attached to the same site. O u r basic philosophy here is to develop a model that is linear in the integer variables. Using the variables defined above, these criteria are expressed as constraints of type equation (2): nmax
Z ~ /-'/ik~ nmax i=lk=l nmax Sma~
(6)
~.
i = I . . . nmax p=lj=l
(7)
k=l
i = 1 Smax
Z
Z
>1 - w,
i = 2 . . . nmax
(8)
p=lj=l nmax nl
nma x
2 Zu,,+ 2w,=nm.x
i=lk=l
W1 = 0 w i <~ w i + 1
(9)
i=1
i = 1 ...(nma x -- 1)
(lO) (11)
48
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programming Model Smax nmax
nmax
~ zijp- ~ z,v,p+M(U,k-- 1)~<0 j = vk p = I Smax
....
k= 1...m
(12)
Srnax
2 Zijp = 2 Zpji j=l"
i= l...n
p= I
1 . . . ( n m a x - - 1),p = (i + 1 ) . . . r i g a x
i:
(13)
j=l
~ x Zijp ~ 1
i = 1 ... nmaxJ : 1 ... Smax
(14)
i=
(15)
p=l
Uik - ~ Ui-1, k 4 0 k=l
2 . . . n m a x.
k=l
Equation (6) puts a limit on the number of structural groups that can constitute a molecule. The double-summation on the left-hand side is equal to n, the actual number of groups present in the molecule. Since a minimum of two groups are required to form a molecule and an upper limit of nmax is specified, 2 ~< n ~< nmax. Equation (7) is an implementation of the Octet rule. The expression on the left gives the number of bonds attached to the ith group, while that on the right states that if the ith group is of the kth type, its valence is vk. This ensures that the number of bonds attached to a group equals the valence of the group. Equations (8), (9), (11) and (10) ensure that only one molecule is formed. Due to the presence of equation (9), equation (6) becomes redundant. This is realized by constraining the ith group to be attached to one of the groups before it, that is, groups 1 to (i - 1) [equation (8)]. Thus the second group has to be attached to the first group, the third group has to be attached to either of the first two groups, and so on. The first group is always present [equation (10)], and the (i + 1)th group is present only if the ith group is present [equation (11)]. Equation (12) has to be introduced to account for different valencies of the groups. This equation is a linear analog of the non-linear equation Srnax
nmax
Uik ~ ~, zijp = 0 J=vk+lp=l
i = 1...n .... k = 1...m
(16)
which states that if the ith group is of the kth kind, then sites (vk + 1) to Smax of that group should not have any connections since they are non-existent. Equation (12) [which is equivalent to Z~=,tlZ~=,,~Zijp + M(U~k -- 1) ~< 0] is written in a form that is convenient from a computer programming point of view. In order to tighten the relaxation, M should be chosen to be the smallest positive number that is at least as large as the maximum attainable value of the first two terms in the equation. In our simulations we simply chose M large enough. Equation (13) is the symmetry constraint. Hence if a group is attached to a second group, the second group is automatically attached to the first one. Equation (14) ensures that a group's site can be attached at most once to some other group. Equation (15) is applied to force existence of the ith group if the (i - 1)th group is present. The structural feasibility constraints [equations (6)-(15)] are linear; hence they form a convex hull separating feasible molecular structures from infeasible ones. One property constraint that has to be incorporated is thermodynamic feasibility of type equation (3). A molecule is thermodynamically feasible at a particular temperature T if its Gibbs' Free Energy of formation Gf a t that temperature is negative, that is Gf(T) < 0.
(17)
For simplicity, Gf is calculated at standard temperature and pressure. While this constraint ensures stability for the generated molecule, it can also eliminate some potentially useful molecules. This may happen if a thermodynamically unstable molecule's dissociation kinetics are such that it can exist for time periods that are long enough for commercial use.
3.2. Characteristics of structural feasibility model The variables used to describe the molecule give nearly complete information about the molecule's composition and connectivity. One of the advantages of this new formulation becomes apparent when one looks at the various group contribution methods. A number of these methods make use of a
International Transactions in Operational Research Vol. 4, N o . I
Cyclic Molecule (one l,oop)
CHIN, / CH 2 CH ~
49
F
z!H
C~ Cyclic Molecule (two independent loops)
Unconnected Molecules (two sub-molecules)
.CH2
H/
CI
CH 2 F + CH 3 F
Fig. 2. Loops and sub-molecules.
certain amount of structural information in giving accurate predictions. For example, the group contribution method for liquid specific heat (Cheuh and Swanson, 1973) requires the addition of 4.5 to the value of Cp~ for any carbon group which fulfils the following criterion: 'A carbon group which is joined by a single bond to a carbon group which is connected by a double or triple bond with a third carbon atom.' It is obvious that this correction cannot be added without any connectivity information. A search of the open literature gave no previous C A M D formulation that was able to give the level of detailed molecular information required for the more accurate and more complex group contribution techniques. Another way in which structural information proves useful is in incorporating bonding constraints. One of the desired properties of a refrigerant is its stability, that is, the compound should not spontaneously decompose or polymerize. In many cases, polymerization can be eliminated by ensuring that the refrigerant does not have double or triple bonds. There is no systematic way of introducing this constraint mathematically in the absence of detailed structural information. With the new formulation, the constraint that has to be added to eliminate multiple bonds is: $
x
zij p <~ 1
i = 1 ...(nma ~ - 1),p = (i + 1)...nma x.
(18)
j=l
From graph theory (Busacker and Saaty, 1965) 'a connected planar graph with n vertices, a edges, and l regions (including the outside or unbounded region)' is defined by the Euler formula l = a - n + 2. We can use a variation of this relation to suppress the formation of a cyclic molecule if we so desire using the equation: l = a - n + b.
(19)
In this equation, l is the total number of independent loops, b is the number of 'sub-molecules' being designed, a is the total number of attachments, and n is the total number of groups. Obviously, b = 1 since we want a single connected molecule. A loop is a series of distinct groups that are attached in such a manner that one can go from one group to others, and return to the first group via a different route. A loop is independent if it cannot be described by a combination of other loops. Two groups are said to be attached if there exists at least one bond between them. These terms are explained with examples in Fig. 2. The total number of attachments (a) and the number of groups (n) are given by a =
Cip i= 1
p=i+ 1
i=Ik=l
(20)
50
N. Churi and L.E.K. A c h e n i e - - M i x e d Integer Non-linear Programmin9 Model OH
OH
OH
OH
N o n - z e r o t e r m s o f 'z' 112 122 211 221
112 122 211 221
Fig. 3. I n d i s t i n g u i s h a b l e molecules.
B a s i s Set 1
CH 3
2
CH 2
3
F
• N o n - z e r o t e r m s o f ' u ' a n d 'z' H
i1 22 33
Z
U
112 211 223 312
12 21 33
Z
113 122 211 311
Fig. 4. Multiplicity.
where {~ Ci p =
i f ESm'xz.. j = 1 13P
= 0
i = 1...(nm, x - 1 ) , p = ( i + 1)...nma x otherwise.
By using these relations in equation (19) and setting l = 0, the resulting molecule will not have any loops and will thus be acyclic. In summary, the formulation gives us a great deal of control over the features of the target molecules. We can distinguish between isomers, allow multiple bonds, and also specify whether or not we desire a cyclic molecule• One small limitation to the connectivity information is that one cannot distinguish between the individual bonds in the case of multiple bonds. For example, the two structures in Fig. 3 cannot be distinguished under the formulation. However, this is not a serious drawback since such cases are not common. In any case, no existing group contribution method is able to distinguish between such structures in a systematic manner. The formulation also does not consider any constraints based on bond angles. While it is possible to incorporate bond angles, we have not found any property prediction methods that can take advantage of this additional information. Figure 4 shows two possible ways in which the variables can describe a molecule. In general, there are several such multiplicities in specifying a molecule. This can lead to increased computational expense. The degree of multiplicity depends on the number of groups present in the molecule and also their valencies. Sine these different solutions describe the same molecule, they have identical objective function values. This fact can be used to eliminate multiplicity by applying a constraint to the objective function value of the type Pne,,, <<-Pola -- ~" However, such a constraint can cause problems with the termination criterion which depends on an increase in the objective value to detect the minima. Also, since other structures may have the same objective function value, there is a risk of
International Transactions in Operational Research Vol. 4, No. I
51
eliminating potentially good molecules. A better approach is to eliminate multiplicity by introducing an integer cut every time a new solution has the same objective function value as that obtained in the previous iteration. The integer cut given below eliminates $ from the search space: ~ - ~ 9~ ~< I BI - 1 iGB
(22)
iEN
where B = {il .~ = I/, N = {i[ ~i = 0 / and ] B ] is the cardinality of B. For the kth iteration, there will be k - 1 such constraints, each eliminating one of the previous solutions from the search space. A more rigorous approach is to introduce constraints that force the sites to be ordered. This will eliminate multiplicity by ensuring that there is only one way in which the variables can describe a molecule.
4. PHYSICAL P R O P E R T Y CONSTRAINTS For refrigerant design, the important refrigerant physical properties and process variables include (i) boiling point, (ii) flammability, (iii) materials and lubricants compatibility, (iv) stability, (v) toxicity, (vi) miscibility, (vii) the saturated liquid (or vapor) mixture pressure in the evaporator (or condenser), (viii) the saturated vapor temperature in the condenser, (ix) the compressor displacement, and (x) the ozone depletion potential. However, we have not used all of these in the case study. Constraints on physical property values are expressed as bounds and are of type equation (3). Thus, suppose one wants the refrigerant properties P1, P2--Pm to lie within a given interval. For example, Pt, P2 and P3 may correspond to the flammability, ODP and boiling point, respectively. This is expressed as: pmin <~Pi <~ Pmax,~/i = 1...m.
(23)
Note that Pi is implicitly a function of ~ which defines the molecular structure. The above equations mean that a property, Pi should lie between the lower and upper bounds p~in and pmax respectively. This allows us to accept refrigerants which have a boiling point within a given range, for instance. In the case study considered, constraints on the evaporator and the condenser pressures (Joback and Stephanopoulos, 1989) have been included. The lowest pressure in the cycle should be greater than the atmospheric pressure; this reduces the possibility of air and moisture leaking into the system. A high system pressure increases the size, weight and cost of equipment (Dossat, 1981). A pressure ratio of 10 is considered to be the maximum for a refrigeration cycle (Perry et al., 1984). Consistent with this, we choose the system pressure Pvp(Te) at the evaporating temperature and the system pressure Pvp(Tcd) at the condensing temperature to satisfy Pvp(Te)/> 1.4 bar
(24)
Pvp(Tcd) <~ 14 bar
(25)
where T~ and Ted are functions of the molecular structure 9A larger enthalpy of vaporization reduces the amount of refrigerant required. The heat of vaporization of Freon 12 is 18.40kJ/mol (Ashrac, 1972). By giving a lower limit of 15.0kJ/mol, the generated molecule will not be significantly worse than Freon 12 in terms of H~. This relaxation is not done to account for errors in the group contribution method, but to prevent elimination of potentially good molecules whose Hv is only slightly lower than that of Freon 2: Hv(T~) >1 15.0 kJ/mol.
(26)
A low liquid specific heat reduces the amount of refrigerant that flashes upon passage through the expansion valve. This is evaluated at T~,a,. The value of the specific heat for Freon 12 obtained from a group contribution method at this temperature is 27.12 cal/mol. As above, the limit is a relaxation of this value: C p ~ ( T ~ ) ~ 30.0 cal/mol.
(27)
By the Clean Air Act of 1995, refrigerants should have ODPs below 0.2 by the year 2000. Hence this value is used as the upper limit for ODP. The ODP of Freon 12 is 1.0.
52
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programming Model ODP <~0.2.
(28)
For a more practical design, materials and lubricants compatibility, toxicity and flammability ought to be incorporated into the mathematical programming model. The constraints discussed so far are by no means exhaustive, and several different ones can be added to achieve a specific refrigerant design objective.
5. CASE STUDY
5.1. Process conditions The following conditions are used for the refrigeration cycle (Duvedi and Achenie, 1996): (i) evaporating temperature Te = - 1.1 °C, (ii) condensing temperature Ted = 43.3 °C, (iii) average temperature Tmea, = 21.1 °C, and (iv) saturated conditions are assumed so that superheat temperature Ts equals T~. The relevant physical properties are hv(T~), Cp,(Tm~,,), Pvp(T¢) and Pvp(Tca). Selection of parameter values depends on the application to which the designed compounds are to be applied. For the case study, the aim is to identify hydrochlorofluorocarbons (HCFCs) that have properties comparable to, or better than, Freon 12, but with significantly lower ODP. The basis set for this case study consists ofCH3, CH2, CH, F and CI. This gives m = 5, Smax = 3, and v = [1 2 3 1 1]. The maximum number of groups in the molecule (nmax)is specified as eight, which corresponds to the highest substitution possible for propane from the groups in the basis set. This limit comes from the observation that higher molecular weight molecules do not have vapor pressure in the range desired for refrigerants. A particular form of the performance objective in equation (5) used in this case study is
f = Cp,(T,ne~,)/nv(Te)
(29)
which is to be minimized. A larger enthalpy of vaporization reduces the amount of refrigerant required. A low liquid specific heat reduces the amount of refrigerant that flashes upon passage through the expansion valve. Pressure ratio, discharge temperature, capacity and COP sensitivity to cycle variability also need to be reflected in the performance criteria for a more realistic cycle. However, in this proof-of-concept study, we will use a simpler model. Notice that equation (3) is non-linear in the integer variables. In the case study considered, t,he physical property constraints are of this form. By simply adding constraints of type ~ = ~ where ~ is continuous, equation (3) becomes non-linear in ~ instead of ~. Also, the objective function becomes a non-linear function of ~. Therefore the sub-problem is a non-linear program for this case study.
5.2. Solution methodology The locally optimal augmented penalty outer approximation (AP/OA) algorithm (Viswanathan and Grossmann, 1990) is employed to solve the MINLP. Initially, a relaxation of the M I N L P is solved. Subsequently, a finite sequence of non-linear programming (NLP) and mixed integer linear programming (MILP) problems are solved. In the N L P sub-problems, derivative information is obtained by perturbation. The flowchart for the algorithm is given in Fig. 5.
5.3. Results and discussions Table 1 shows the results obtained for the problem described above. The values of the property constraints for the molecules formed are also given. The first row gives the locally optimal solution to the relaxed NLP. However, it has no physical significance since the solution variables are not integers. The first feasible molecule formed is CH3--CI with an objective function value of 9.8513. The next iteration gives F - - H C = C H - - F with an objective function value of0.8911. Since there is an increase in the objective value, iterations stop at this point with CH3--C1 as the locally optimal molecule.
International Transactions in Operational Research Vol. 4, No. 1 Initial Guess
53
]
I ] RelaxedNLP I Yes
I
I 'Su ro lem I
Fig. 5. AP/OA algorithm: minimization problem.
Once a locally optimal solution is found, the next best solution is found by introducing an integer cut in the formulation. In this manner, it is possible to get a series of solutions in decreasing order of optimality as given in Table 2. Since a globally optimal solution is not obtained for the NLP sub-problem, it is conceivable that the sequence of local minima may not be monotonic. We did not however encounter this problem in the case study. As can be seen from Table 2, not all molecules may be suitable for practical purposes. Molecules like F - - H C = C H - - F may polymerize over the long run and hence cannot be used as refrigerants. It is possible to eliminate such molecules by introducing bonding constraints [equation (18)] into the formulation. In any case, the formulation can give a series of molecules that conform to all applied constraints and which can be transferred to the experimental stage for further evaluation. Due to the non-convexity of the formulation, the global optimum cannot be guaranteed when using the AP/OA algorithm for solving the MINLP model. However, current advances in global optimization [see, for example (Floudas and Visweswaran, 1993)] can be employed to obtain globally optimal solutions to the MINLP formulation. As the problem size increases, there will be a need to consider decomposition and parallel branch and bound algorithms (Pekny and Miller, 1992). The formulation for the single component refrigerant is somewhat computationally intensive. For the case study considered, there are 240 binary variables, 124 linear constraints and 6 non-linear constraints. On an IBM RISC 320H workstation the computation time was of the order of 10-30s for the initial relaxed problem, and 10-15min on the subsequent MILP/NLP sequence. This is mainly due to the bond multiplicity alluded to earlier. To alleviate this problem, an integer cut was Table 1. Single component refrigerant design Objective
h,(T,)
P,p(T~)
Itn.
Molecule
Cp,/h,
(kJ/mol)
Cp,(T,,,,,) (cal/mol)
(bar)
P,p(Tcd) (bar)
OOP
0 1 2
-CH3--CI F--HC=CH--F
0.5726 0.8513 0.8911
19.23 20.44 20.20
11.01 17.40 18.00
3.77 2.44 2.37
14.00 9.17 9.92
0.01 0.00
Table 2. Single component refrigerant design: series of feasible molecules Objective No.
Molecule
Cp,/h,
1 2 3 4
CH3--CI F--HC=CH--F CH3--CH2--F F~H(CH2) 2
0.8513 0.8911 1.1282 1.2458
h,( T,)
C p,(Tmean)
(kJ/mol)
(cal/mol)
Pvo(Te) (bar)
Pvp(Tcd) (bar)
ODP
20.44 20.20 17.78 18.88
17.40 18.00 20.06 23.52
2.44 2.37 4.19 3.51
9.17 9.92 13.07 13.72
0.01 0.00 0.00 0.00
54
N. Churi and L.E.K. Achenie--Mixed Integer Non-linear Programming Model
introduced every time bond multiplicity was detected. While it is possible to eliminate bond multiplicity, the additional constraints are either non-linear, or if linear, result in an increase in the problem size. At any rate, we are formulating integer constraints to deal with this problem in a rigorous way.
6. CONCLUSIONS The model proposed for single component refrigerant design (with environmental considerations being incorporated) has been shown to be effective. The proposed structural feasibility constraints for the single component refrigerant design give nearly complete molecular information; thus we can easily accommodate any group contribution based property prediction techniques that require detailed molecular bonding information. Sensitivity of the analysis to uncertainties in the physical property estimation is very important but falls outside the scope of this paper. These will need to be addressed in the future through sensitivity analysis.
REFERENCES American Society of Heating, Refrigeration and Air-Conditioning Engineers, Inc. (1972). ASHRAE Handbook of Fundamentals. Brignole, E.A., Bottini, S. & Gani, R. (1986). A strategy for the design of solvents for separation processes. Fluid Phase Equilibria 29 125-132. Busacker, R.G. & Saaty, T. L. (1965). Finite Graphs and Networks: An Introduction with Applications. McGraw-Hill, New York. Chueh, C.F. & Swanson, A.C. (1973). Estimation of liquid heat capacity. Canadian Journal of Chemical Engineering 5, 596-600. Dossat, R.J. (1981). Principles of Refrigeration. Wiley, New York. Duvedi, A. & Achenie, L. E. K. (1996). Design environmentally safe refrigerants using mathematical programming. Chemical Engineering Science 51, 3727-3739. Floudas, C.A. & Visweswaran, V. (1993). A primal-relaxed dual global optimization approach. J. Optim. Theory and Appl. 78(2), 187. Gani, R. & Brignole, E.A. Molecular design of solvents for liquid extraction based on UNIFAC. Fluid Phase Equilibria 13, 331-340. Gani, R., Nielsen, B. & Fredenslund, Aa. (1991). A group contribution approach to computer-aided molecular design. AICHE J. 37, 1318-1331. Joback, K. G. & Stephanopoulos, G. (1989). Designing molecules possessing desired physical property values. In Proceedings FOCAPD 1989, Snowmass Village, Colorado, FOCAPD. Odele, O. & Macchietto, S. (1993). Computer aided molecular design: a novel method for optimal solvent selection. Fluid Phase Equilibria 47-54. Pekny, J. F. & Miller, D. L. (1992). A parallel branch and bound algorithm for solving large asymmetric travelling salesman problem. Mathematical Programming 55, 17-33. Perry, R.H., Green, D.W. & Maloney, J. O. (1984). Perry's Chemical Engineers' Handbook, sixth edn. McGraw-Hill, New York. Viswanathan, J. & Grossmann, I.E. (1990). A combined penalty function and outer-approximation method for MINLP optimization. Comp. & ChE. 14, 769-782.