Computer-Aided Solvent Design for Reactions: Maximizing Product

Brooke , A. ; Kendrick , D. ; Meeraus , A. ; Raman , R. GAMS A User's Guide; GAMS Development Corporation: 1998. There is no corresponding record for ...
0 downloads 0 Views 230KB Size
Download PDF
5190

Ind. Eng. Chem. Res. 2008, 47, 5190–5202

Computer-Aided Solvent Design for Reactions: Maximizing Product Formation Milica Folic´,† Claire S. Adjiman,* and Efstratios N. Pistikopoulos Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom

A hybrid experimental/computer-aided methodology for the design of solvents for reactions, recently proposed by the authors [Folic´ et al., AIChE J. 2007, 53, 1240–1256], is extended. The methodology is based on the use of a few reaction rate measurements to build a reaction model, followed by the formulation and solution of an optimal computer-aided molecular design (CAMD) problem. The treatment of complex reaction systems, such as competing or consecutive reactions, is considered through the incorporation of a simple reactor model in the problem formulation. This approach is applied to two model reaction schemes, and it is shown that, in principle, it is possible to identify solvents that maximize product formation by enhancing the main reaction and suppressing byproduct formation. Since very few measurements are used to build the reaction model, the effect of uncertainty is tackled explicitly in a stochastic formulation of the CAMD problem. An approach to sensitivity analysis for the identification of the key model parameters is discussed. Using this information to generate scenarios, a stochastic optimization problem (whose objective is to determine the solvents with the best expected performance) is then solved. The final output consists of a list of candidate solvents that can be targeted for experimentation. The methodology is demonstrated on a Menschutkin reaction, which is a representative SN2 reaction. This shows that the uncertainty in the reaction model has little impact on the types of solvent molecules that have the best performance. Dinitrates are found to be a promising class of solvents, with regard to maximizing the reaction rate constant. 1. Introduction Given the important role played by solvents in chemical processing, significant research effort has been devoted to the development of systematic methods for solvent selection. Computer-aided molecular design (CAMD) has emerged as an attractive route for solvent selection. Several CAMD approaches have been proposed in the last two decades.1 The two most significant classes of techniques are “generate-and-test” and optimization-based methods. These have been successfully applied to a variety of solvent-based separation problems, allowing a much larger number of solvent molecules to be considered during separation system design than is possible by experimentation alone. On the other hand, there has been little work on the design of solvents that can enhance reaction rates, despite the significant gains that can be achieved by optimizing the choice of reaction medium.2 Reaction rate constants and the solubility of the substrates, catalysts, and products can vary by several orders of magnitude from solvent to solvent. In complex systems, a judicious choice of solvent has the potential to accrue large benefits; for instance, it may be possible to influence selectivity and thereby minimize byproduct formation. Solvents have been reported to affect the regioselectivity of 1,3-dipolar cycloadditions, particularly between nitrones and substituted alkenes.3 It may also be possible to reduce the number of processing steps required by “telescoping” them (that is, by carrying several reactions in a single solvent). This leads to a reduction in costs and, often, emissions. An example of the benefits that can be achieved can be found in the work of Rasmy et al.4 However, systematic solvent design methods for these complex problems are lacking, and industrial practice is currently mostly based * Author to whom correspondence should be addressed. Tel: +44 (0)20 7594 6638. E-mail: [email protected]. † Current address: R&D-Environmental, Haldor Topsoe A/S, Nymollevej 55, 2800 Lyngby, Denmark.

on experience and intuition when it comes to solvent choice during the development of new reaction routes. Recently, Gani et al.5 proposed a method for solvent selection for the promotion of organic reactions that combines knowledge from industrial practice and physical insights. CAMD is used to generate a list of solvent candidates, ranked according to their score. This method has been proven to be very effective for some application studies, but it requires a significant amount of information on both the solvents and the reactions to build a table with solvent scores for each reaction. Although it takes a large number of issues that affect performance into consideration, it is not possible to predict the solvent effect on the reaction rate in a quantitative manner. Using an empirical model of solvent effect on reaction rate constants, in our previous work,6–8 we have developed a systematic approach to solvent design for reactions, with the objective of maximizing the rate constant for a single reaction. This method is based on targeted experiments, model development, and candidate generation by optimization-based CAMD. The optimization problem is formulated so that a range of design criteria can easily be incorporated within the problem and considered simultaneously during the design phase. This allows the identification of optimal tradeoffs. The model of solvent effects on the reaction kinetics takes on a simple algebraic form and therefore is quick to evaluate. It can identify solvent rankings reliably, but it generally does not provide a high degree of quantitative accuracy. Stanescu and Achenie9 have proposed a two-step method based on a CAMD step in which candidate solvents are generated based on constraints on their physical properties, followed by density functional theory (DFT) solvation calculations to estimate the reaction rate and product yield in the candidate solvents. They have applied this approach to the Kolbe-Schmitt reaction, and they found that the product yield is largest without solvent. In solvents, the yield is limited by the reversibility of the reaction, especially in solvents with a

10.1021/ie0714549 CCC: $40.75  2008 American Chemical Society Published on Web 06/14/2008

Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 5191

high dielectric constant. The use of DFT in this approach can lead to a better accuracy than that which is possible with simpler models,8 but it is very time-consuming. These different attempts at developing CAMD approaches that are applicable to reactive systems highlight the central role of predictive models that can relate kinetics to solvent structure. These range from the empirical solvatochromic equation of Abraham and co-workers10–12 to ab initio approaches, such as those discussed by Cramer.13 These classes of models pose distinct challenges, in terms of their incorporation within CAMD, because of their varying levels of complexity, computational cost, and reliability. The objectives of this paper are to extend our earlier methodology8 to more-complex reaction systems and to test this through several examples. In our earlier work, the approach was applied successfully to an SN1 reaction, where it was shown that the uncertainty in the empirical model has little effect on the solution. While the results obtained were promising, further experience is required to test whether the proposed framework is broadly applicable across reaction classes. The empirical model, in which the solvatochromic equation and group contribution techniques are combined, is briefly reviewed in section 2. The formulation is then extended to more-complex reaction networks, with special emphasis on competing (or concurrent) reactions, in which more than one transformation is available for a given reactant, and consecutive reactions in which the product of the first reaction is an intermediate in the formation of the final product, or in which the desired product undergoes further transformation in the second reaction. This is the focus of section 3. The approach is then applied to model systems of competing and consecutive reactions in section 4. A methodology for handling uncertainty within the CAMD formulation is presented in detail in section 5. In section 6, it is applied to a Menschutkin reaction, which is typical of the SN2 class. 2. Reaction Model At the core of the proposed solvent selection methodology8 is a two-step approach. The first step is concerned with the development of a model of solvent effects on the reaction, and the second step involves the formulation and solution of an optimization-based computer-aided solvent design problem, which can be either deterministic or stochastic. The first step of the methodology is briefly reviewed in this section. The key issue in the first step (model development) is to identify a relationship that links solvent structure to the reaction rate in a way that quantifies solvent effects on a particular reaction. To do this, we use the multiparameter solvatochromic equation,12 which correlates properties of the solvent such as the empirical solvatochromic parameters (A, B, and S), a polarizability correction term (δ), and the Hildebrand solubility parameter (δH2 ) with the logarithm of the reaction rate constant: hδH2 logk ) logk0 + sS + dδ + aA + bB + 100

(1)

where logk0, s, d, a, b and h are reaction-specific constants that quantify the dependence of the rate constant on each solvent property. Following the notation of Zissimos et al.,14 A is the hydrogen bond acidity, B is the hydrogen bond basicity, and S is the dipolarity/polarizability. Two types of solvatochromic parameters are reported in the literature for a given compound, depending on how the measurements were obtained: “solute” and “solvent” values. In this work, solute values are used,

because many more have been tabulated. The work of Li et al.15 has shown that solute values can be used successfully to model the solvent in the Smx solvation model. Our earlier work on the solvolysis reaction has also shown the suitability of this approach.8 In the remainder of this paper, S, A, B, δ, and δH2 are called solvent properties, and logk0, s, d, a, b, and h are called reaction parameters. An important characteristic of the solvatochromic equation is that the reaction parameters and solvent properties are independent of each other. Thus, for each reaction, there is a fixed set of reaction parameters. Similarly, each solvent has a fixed set of properties, regardless of the reaction considered. The equation has been used previously to quantify solvent effects on the solvolysis of t-butyl halides11,12,16 and on Diels–Alder reactions.17 Both reaction classes have long been a key reference for theories of solvent effects on organic reaction rates. To obtain the reaction parameters, the model-building step involves gathering experimental data for a small set of predetermined solvents and generating a reaction model through regression. First, a small set of initial solvents for the reaction studied must be assembled. Although the number of solvents can be chosen arbitrarily, we have found that data in eight solvents yields sufficient information to build a suitable model, provided that the solvents are diverse, in terms of their polarity. As a measure of this, we use the ETN solvent polarity scale18 and we choose solvents that have ENT values distributed over the entire physical range. In addition, it is preferable to choose solvents with different functional groups. Wherever possible, literature data should be used at this stage to minimize experimental costs. In the absence of reliable data, experimental reaction rate constants for the solvents chosen are measured. To complete the relationship between solvent structure and reaction rate constant, solvent property values predicted from a knowledge of the molecular structure are used in the solvatochromic equation. The polarizability correction term (δ) can be calculated exactly based on molecular structure. The other solvent properties are correlated to their molecular structure using group contribution (GC) methods and correlations that are widely used for property prediction within CAMD-based frameworks. These methods are based on the principles of transferability and additivity, with atom groups such as CH2 and OH, used as building blocks. For further details on the specific methods used for solvent property prediction, the reader is referred to Folic´ et al.8 Forty-three functional groups (42 of which are UNIFAC groups) are available to construct solvent molecules. The use of UNIFAC first-order groups to build solvents makes integration with other solvent design approaches easier. The kinetic data and predicted solvent properties are used to regress the reaction parameter values (logk0, s, d, a, b, and h) and thereby obtain a model of solvent effects on the chosen reaction. The resulting model can predict the reaction rate constant in solvents other than those tested experimentally, based solely on their molecular structure. In addition to the optimal values of the reaction parameters which are used in the model, the regression provides confidence intervals for the parameters. These make it possible to quantify the uncertainty in the model. 3. Deterministic CAMD Problem Formulation In the second step of the solvent design methodology, a computer-aided solvent design problem is formulated based on the predictive model of solvent effects. Two formulations can be used. Using the optimal, or nominal, values of the reaction parameters, a deterministic formulation is derived. Alternatively,

5192 Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008

the confidence intervals can be used to generate a stochastic formulation that takes the uncertainty into account to arrive at a design. In this section, the deterministic formulation is considered. In this case, a mixed-integer problem (MIP) is posed, with the objective to identify a solvent in which the reaction rate constant under given conditions or a function of the reaction rate constant, is maximized. This may be product formation, yield, or another performance-related objective. The problem may be linear (MILP) or nonlinear (MINLP), depending on the objective function, as will be shown in sections 4 and 6.2 that involve these case studies. The generic formulation of this problem is given below: max f(x_, logk) s.t. h_1(x_, p_, logk) ) 0 g_1(x_, p_, logk) e 0 _e(logk, p_) ) 0 h2(n_, p_, _y) ) 0 g_2(n_, p_, _y) e 0 h_3(n_, _y) ) 0 g_3(n_, _y) e 0

pendix A. Based on the molecular complexity constraints given and the set of atom groups available, 9759 molecules are feasible. These include many chemical classes, such as alcohols, diols, triols, nitrates, amines, carboxylic acids, ethers, esters, aldehydes and ketones, both aliphatic and aromatic. An even larger space can be explored if some of the restrictions in molecular complexity constraints are lifted. We also include integer cuts (see Appendix A) in the formulation to allow the generation of successive solutions, giving a ranked list of candidate solvents. The introduction of process constraints and a generic objective function allows the consideration of more-complex reaction schemes, as shown in sections 3.1 and 3.2. 3.1. Competing Reactions. For competing reactions, we consider a model system of a continuous stirred tank reactor (CSTR) with a feed stream that consists of reactant C and a solvent. The following first-order reactions occur: k1

C 98 D (P1)

_x ∈ R p_ ∈ Rm n_ ∈ Rq logk ∈ Rr t

u _y ∈ {0, 1} where f is the objective function; h1 is a set of process model equality constraints; g1 is a set of process inequality constraints; e is a set of solvatochromic equations (one equation for each rate constant); h2 a set of structure–property equality constraints; g2 is a set of structure–property inequality constraints; h3 is a set of chemical feasibility and complexity equality constraints; g3 is a set of chemical feasibility and complexity inequality constraints; x is a t-dimensional vector of process variables (e.g., concentration); logk is an r-dimensional vector of logarithms of reaction rate constants; p is an m-dimensional vector of continuous variables denoting physical properties; n is a q-dimensional vector of continuous variables denoting the number of groups of each type in the molecule; and y is a set of binary variables (e.g., used to constrain the n variables to integer values). Structure/property constraints include the property prediction methods mentioned in section 2. Chemical feasibility constraints include the octet rule19 and the modified bonding rule,20 to ensure there are no free attachments and no double bonds between atom groups. These rules ensure that the molecule designed is structurally feasible, but they do not guarantee that the molecule is chemically stable. Instead, stability is assessed as part of the post-processing of results. Molecular complexity constraints include standard constraints on the maximum and minimum number of groups in the molecule, constraints on the maximum number of main and functional groups, as well as constraints that forbid or limit the joint occurrence of some groups. The choice of these constraints is based on the recognition that solvents are usually medium-size molecules, because they must have a useful liquid range. In addition, the group contribution methods used here do not account for proximity effects beyond the makeup of atom groups and are, therefore, most reliable when used for medium-size molecules. Structure–property and chemical feasibility and complexity constraints for this problem formulation are given in Ap-

k2

C 98 E where D is the desired product whose concentration should be maximized and E is the side product whose production should be as small as possible, compared to the production of D. The reaction rate constants corresponding to the desired reaction and side reaction are denoted by k1 and k2, respectively. The following objective function can be used in formulation P1: maxCD - CE

(2)

where Ci denotes the outlet concentration of species i (in units of mol/m3). For this case study, the process model (h1) equations in formulation P1 are as follows: k1 ) 10logk1

(3)

k2 ) 10logk2

(4)

CC )

CC0 1 + k1 + k2

(5)

CD ) k1CC

(6)

CE ) k2CC

(7)

Here, the x (process) variables are given as follows: CC, the outlet concentration of C; CD, the concentration of the desired product D; and CE, the concentration of side-reaction product E. There are two logarithms of the rate constants: logk1 and logk2. The parameter CC0 denotes the inlet concentration of reactant C. The remainder of the formulation is as presented in Appendix A. The resulting problem is a MINLP program, which is linear in the binary variables. It can be solved locally in GAMS21 with the outer-approximation algorithm.22 3.2. Consecutive Reactions. Another reaction system commonly encountered in industry is the case of consecutive reactions, where the reaction that yields the desired product is followed by another reaction, which consumes the desired product and, therefore, must be demoted, such as

Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 5193 k1

k2

Table 2. List of Solvents Used in Regression for Case 1b and the Chosen Rate Constant Valuesa

C 98 D 98 E where C is the reactant, D the desired product, and E is the undesired byproduct; k1 and k2 are rate constants. One possible objective is the maximization of the difference between the concentrations of products D (CD) and E (CE). The following set of process model (h1) equations must be added to the problem formulation to consider this objective: k1 ) 10logk1

(8)

k2 ) 10

(9)

logk2

CC )

CC0

(10)

1 + k1

CD )

k1CC 1 + k2

(11)

CE )

k1k2CC 1 + k2

(12)

Here, the x (process) variables are given as follows: CC, the outlet concentration of C; CD, the concentration of the desired product D; and CE, the concentration of side-reaction product E. There are two logarithms of rate constants: logk1 and logk2. Parameter CC0 is the inlet concentration of reactant C. 4. Application of Deterministic CAMD to Model Reaction Schemes The proposed deterministic formulation is applied to the reaction schemes discussed in sections 3.1 and 3.2. 4.1. Competing Reactions. Two model reaction systems are considered to demonstrate the application of the formulation from section 3.1. First, we consider a system of competing reactions where the desired reaction is favored in protic solvents while the side reaction is favored in aprotic solvents (case 1a). We then consider a case where the desired reaction rate is enhanced in polar solvents, while the side reaction is faster in apolar solvents (case 1b). The eight solvents used for regression were chosen based on their hydrogen-bond donor acidity value (A), for case 1a, or based on their polarity (ETN scale) for case 1b, so that half of the solvents are protic (or polar) and half of the solvents are aprotic (or apolar), for case 1a (or case 1b). The lists of solvents used in cases 1a and 1b are given in Tables 1 and 2, respectively. Logarithms of reaction rate constants were attributed in such a way that solvents in which a reaction is favored have high values of the rate constants, and solvents that do not affect the reaction favorably are assigned lower values. To be more realistic, the rate constant data for case 1a, shown in Table 1, were not ranked exactly in the same order as Table 1. List of Solvents Used in Regression for Case 1a and the Chosen Rate Constant Valuesa

a

solvent

A

logk

glycerol propane-1,2-diol aniline acetic acid acetone 1,2-dichloroethane benzene pentane

0.96 0.64 0.27 0.61 0.00 0.00 0.00 0.00

-0.5 -1.0 -1.5 -1.7 -4.0 -5.0 -6.0 -7.0

Protic solvents are assigned high experimental rate constants, because they are favored for the desired reaction.

solvent

EN T

log(k)

phenol glycerol propane-1,3-diol N-methylformamide dioxane benzene diethyl ether pentane

0.95 0.81 0.75 0.72 0.16 0.11 0.12 0.01

-0.5 -1.0 -1.5 -1.7 -4.0 -5.0 -6.0 -7.0

a Polar solvents are assigned high experimental rate constants, because they are favored for the desired reaction.

the A values. For instance, aniline, which has a lower value of A than acetic acid, has a slightly higher rate constant. The initial concentration of reactant C (CC0) is assigned a value of 1.5 mol/ m3. The reaction coefficients for the side reactions in both cases were obtained by assigning the reaction rate constants to the same solvents as those used for the main reaction, but in reverse order. The following equations were obtained for the two competing reactions in case 1a, through two linear regressions: logk1 ) -2.26 + 11.00S - 2.13δ + 24.19A - 1.23B -

9.76δH2 100 (13)

6.05δH2 logk2 ) -3.19 - 4.82S - 0.09δ - 16.52A - 1.85B + 100 (14) The regression and prediction statistics obtained (based on the eight solvents that were used for regression) are, for eq 13, R2 ) 0.97, SE ) 0.77, and AAPE ) 17.50%, and, for eq 14, R2 ) 0.99, SE ) 0.30, and AAPE ) 6.40%. (Note: AAPE stands for average absolute percentage error, which is defined N as AAPE (%) ) (1/N)∑i)1 (|Xpred,i - Xexp,i|/Xexp,i) × 100, where N is the number of compounds in the dataset, Xexp,i is the experimental value of the property (log(k)) for compound i and Xpred,i is the predicted value of the property for compound i.) Comparing the coefficients in eqs 13 and 14, one can see that the largest absolute difference is in reaction parameter a, which is multiplied by the hydrogen bond donor acidity A; this is the property that discriminates between protic and aprotic solvents. Equations 15 and 16 were obtained for case 1b: logk1 ) -4.85 + 3.05S + 1.17δ + 12.57A + 2.18B -

4.26δH2 100 (15)

logk2 ) -4.58 - 0.70S - 1.39δ - 21.76A - 2.85B +

7.49δH2 100 (16)

The regression and prediction statistics obtained (based on the eight solvents used for regression) are, for eq 15, R2 ) 0.99, SE ) 0.40, and AAPE ) 11.06%, and for eq 16, R2 ) 0.96, SE ) 0.83, and AAPE ) 26.19%. Equations 13 and 14 were included in the problem formulation (these are the e equations in formulation P1). Solving the resulting MINLP 100 times, with the addition of an integer cut after each solution, resulted in a ranked list of 100 solvents. The highest objective values were obtained for protic solvents

5194 Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008

(alcohols and amines). Aliphatic hydrocarbons and esters had the lowest objective function values. For example, the reaction giving product E was completely suppressed in 1,2-dihydroxy2-propene, whereas in ethane, it was more favorable than the reaction that gave product D (CE was higher than CD by 0.22 mol/m3). Similarly, for case 1b, eqs 15 and 16 are included in the problem formulation. The highest values for CD (e.g., for 1,2dihydroxy-2-propene, a value of 1.5) are obtained for alcohols (diols and triols), which are polar solvents, and the lowest (close to zero, and lower than the values of CE) are obtained for aromatic hydrocarbons and halosubstituted hydrocarbons, which are apolar solvents. 4.2. Consecutive Reactions. For the consecutive reactions, two cases are considered; when the desired reaction is favored in protic solvents while the consecutive reaction is favored in aprotic solvents (case 2a) and a case when the desired reaction is promoted and the consecutive reaction demoted in polar solvents (case 2b). The reaction models given by eqs 13–16 were used to predict the logarithms of reaction rate constants for the two cases. Solving the deterministic optimization problem for case 2a 100 times, with the addition of an integer cut after each solution, yielded a list of 100 molecules. By ranking those molecules according to the performance measure, the difference in performance between protic and aprotic solvents is highly visible. The highest objective function values belong to alcohols (diols and triols), followed by the N-monosubstituted amides, whereas the hydrocarbons and thiols (aprotic solvents) are located at the bottom of the list, with objective function values very close to zero. In some solvents (for example, 2-methylene1,8-octanediol), the reaction that yields product C is entirely suppressed. Similar results are obtained for case 2b; the highest objective function values (e.g., 1.47 mol/m3 for 1,1,2-butanetriol) were obtained for polar solvents, mainly alcohols (diols and triols), and the lowest objective function values were observed for halosubstituted aromatic and aliphatic hydrocarbons (apolar solvents). For those solvents, the objective function value is close to zero and less than the value of CE, which means that the reaction that yields E is more favorable than the reaction that yields product D. 5. CAMD under Uncertainty One issue that arises in the formulations presented is the fact that the reaction model is based on kinetic data in a few solvents only, so significant uncertainty is likely to be associated with the reaction parameters. This is confirmed by the 95% confidence interval that is calculated for each of the reaction parameters, which are usually very wide. Therefore, it is instructive to investigate the impact of uncertainty on the deterministic design. A two-step strategy has been proposed8 to quantify the effect of uncertainty on model reliability and to determine the optimal solvent candidate, given this uncertainty. The first step is determine the key reaction parameters and the “optimal” solvents most frequently identified across the uncertain parameter space. To ensure that any correlation between the different reaction parameters is captured, a global sensitivity analysis23 is performed. A stochastic (scenario-based) optimization problem then is formulated and solved, using representative scenarios distributed throughout the uncertainty range/space. The choice of these scenarios is based on the results of sensitivity analysis, as discussed in the remainder of this section.

5.1. Sensitivity Analysis. Global sensitivity analysis is performed first by sampling the uncertain parameter space and solving the deterministic design problem for each combination of the parameters sampled. A uniform distribution is used to sample the space if only one parameter is varied. When varying more than one parameter simultaneously, we use an approach proposed by Sobol’.24,25 The Sobol’ approach is based on a quasi-random generation of numbers, which allows a progressive and efficient sampling of multidimensional spaces. The main sensitivity metric considered is the number of different optimal solvents found by solving the deterministic design problem repeatedly for different reaction parameter (Ns) values. Given the large size of the design space (nearly 10 000 solvents), if uncertainty has a significant impact, it is expected that a large number of different solvents will be found. The occurrence frequency of each optimal solvent, which is defined as the number of times a given solvent is found divided by the total number of samples, also gives an indication of the importance of uncertainty. If the solution space is fragmented, with all optimal solvents recurring with low frequency, then uncertainty affects the solution greatly. If, on the other hand, a few solvents occur more frequently than others, the impact of uncertainty is not as significant. For a problem including a single reaction, the full uncertainty space for the design problem is five-dimensional and consists of parameters s, d, a, b, and h. The value of logk0 can be kept constant, because it does not have an impact on the optimal solvent structure. When there are several reactions, all six reaction parameters must be included in the uncertainty analysis. Therefore, a large number of sample points may be necessary to obtain a good estimate of the sensitivity metrics. It is possible to reduce the resulting computational cost by first computing the impact of each parameter individually and using this information to reduce the dimensionality of the problem. This is done through a one-dimensional global sensitivity analysis, where the value of one parameter is varied uniformly over its range, while keeping all other parameters constant at their nominal values. In this case, the relative importance of parameter θ is computed as Sθ )

Ns,θ U θ95 - θL95

(17)

where Ns,θ is the number of different optimal solvents designed U L over the set of sampled values of θ and θ95 and θ95 are the upper and lower limits of the 95% confidence intervals, respectively (for example, see Table 7, presented later in this paper). The parameters that have a larger Sθ value are key parameters, and the higher-dimensional sensitivity analysis can be focused on the reduced space of key parameters. While this analysis allows a reduction on complexity by focusing on key parameters, it ignores any correlation between parameters. This can only be captured by performing a multidimensional analysis, as described later in this paper. Strictly speaking, the evaluation of the sensitivity metric Ns is completed after all of the new designs have been identified. This is an infinite-dimensional problem, and an alternative convergence criterion is considered which is based on setting a limit on the rate at which new designs are being identified. A rate of 1 design per 100 sample points is used in this work, because, in our experience, this has proven to be a good indication that sufficient solvent diversity has been achieved. A byproduct of the sensitivity analysis is a virtual “solution map”, that is, a record of which optimal solutions occur over different regions of the multidimensional parameter space. In

Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 5195

the two-dimensional case, this map can be visualized, as shown in Figure 1 where eight solvent molecules, represented by eight different symbols, are shown. The map shows some dominant designs, which cover relatively large regions of space (e.g., the triangle), and less-common designs (such as the cross and rhomboid). In the case where the design problem is linear, it can be shown that each design occurs over a convex region of the parameter space or a union of disjoint convex regions.26 The virtual solution map is used to identify a set of representative scenarios, as discussed in the next section. 5.2. Scenario Generation and Problem Formulation. After we have gathered and analyzed the results of the sensitivity analysis, we can proceed to formulate and solve a stochastic optimization problem over the uncertain parameter space. The objective function is the expected performance of the solvent, and it is approximated by taking the average performance over a finite set of scenarios. The problem is formulated for a small number of representative scenarios that are distributed throughout the uncertainty space. For this purpose, the solution map obtained by sensitivity analysis is used. A convex hull is identified for each optimal solvent and one or more scenarios are selected within each hull. A weight factor is assigned to each scenario based on the fraction of the volume of the overall space represented by the scenario. Further details can be found in Folic´ et al.8 When the representative scenarios have been determined, we can formulate and solve the stochastic optimization problem, for the purpose of maximizing the expected value of the performance objective over the entire uncertain parameter space. The problem is written as follows: max

_x,logk,n_,y_

s.t.

M

∑

1 w f(x_, logki) M i)1 i h_1,i(x_, p_i, logki, n_, _y) ) 0

(i ) 1,

. . . , M)

g_1,i(x_, p_i, logki, n_, _y) e 0 (i ) 1, . . . , M) _e(logki, p_i) ) 0 (i ) 1, . . . , M) h2,i(n_, p_i, _y) ) 0 (i ) 1, . . . , M) g_2,i(n_, p_i, _y) e 0 (i ) 1, . . . , M) h_3(n_, _y) ) 0 g_3(n_, _y) e 0 _x ∈ Rt p_i ∈ Rm×M _n ∈ Rq log(ki) ∈ Rr×M u _y ∈ {0, 1}

(P2) where wi is the weight factor for scenario i and M represents the number of scenarios. All other symbols are as previously defined, and all the constraints are as given in Appendix A. The subscript i indicates which scenario is being considered and the weight factors in the formulation correspond to the ratio of the area (volume) of the particular molecular design convex hull to the total sum of the area (volume) of convex hulls. The solution obtained via stochastic optimization is more reliable than the one obtained by deterministic optimization in the sense of average performance. We also include integer cuts in the formulation (see Appendix A.1) to allow the generation of successive solutions, giving a ranked list of candidate solvents.

Figure 1. Map of solutions resulting from sensitivity analysis in a twodimensional space (θ1 × θ2). Each symbol represents a different molecule.

6. Application of Stochastic CAMD to a Menschutkin Reaction Both the deterministic and the stochastic problem formulations are applied to a Menschutkin reaction case study. Menschutkin reactions are SN2 reactions between haloalkanes and tertiary amines, yielding quaternary ammonium salts, and they are known to be highly influenced by the change in reaction medium.27 Quaternary ammonium salts are ionic liquids, which means they contain only ionic species, not molecular species. The particular reaction studied here is between tripropylamine ((C3H7)3N) and methyliodide (CH3I): (C3H7)3N + CH3I f CH3(C3H7)3N+ + IReaction rate data were collected from the kinetic study published by Lassau and Jungers.28 In total, data for 59 solvents were gathered. Significant solvent effects were observed, as can be seen in Table 3. A theoretical study of solvent effects on Menschutkin reactions, using ab initio calculations and Monte Carlo calculations, has been performed by Castejon and Wiberg.29 The authors also conducted a study of the reaction of amines with trimethylsulfonium salts, using ab initio methods, which predicted the effect of two solvents (water and DMSO) very well.30 6.1. Reaction Model Building. The set of 59 solvents studied by Lassau and Jungers28 is very diverse. The solvents are presented in a rank-ordered list in Table 3, based on the experimentally measured rate constant values. Eight solvents were chosen to build the reaction model, based on their diversity in polarity and in functional groups. These are listed in Table 4. The remaining 51 solvents are used for verification of the model (that is, for obtaining the statistics upon extrapolation of the model to the entire set of 59 solvents). The following solvatochromic equation was obtained through regression performed with the eight solvents from Table 4: logk ) -5.73 + 3.80S - 0.03δ + 3.17A + 0.09B -

0.44δH2 100 (18)

The statistics obtained for this regression were satisfactory. R2, which is the square of the Pearson product moment correlation coefficient, is 0.86. The average absolute percentage error (AAPE) calculated after extrapolation of eq 18 to predict logkvalues for all 59 solvents in the data set is 18.77%. The rank of each solvent is also reported in Table 3, where 1 denotes the solvent with the largest rate constant, and 59

5196 Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 Table 3. Solvent Number, Solvent Name, Experimental logk Value, Predicted logk Value, and Associated Ranking for the 59 Solvents for the Menschutkin Reaction

Table 4. Eight Solvents Used to Build the Reaction Model and Their ETN Solvatochromic Polarity Values for the Menschutkin Reaction

Solvent number

name

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

benzyl cyanide 1,1,2,2-tetrachloroethane N,N-dimethylformamide nitroethane acetophenone 1,2-dichloroethane benzaldehyde 2,5-hexanedione phenylpropanone phenyl-4-butanone-2 1,1,2-trichloroethane nitropropane propionitrile cyclopentanone butyronitrile 1-methylnaphthalene 1,4-dichlorobutane cyclohexanone acetone iodobenzene 1,2-dibromoethane 2,4-pentanedione o-dichlorobenzene 2-butanone anisole bromobenzene 1,1-dichloroethane ethanethiol chlorobenzene benzylalcohol styrene ethoxybenzene 3-heptanone dioxane m-dichlorobenzene allylchloride tetrahydrofuran bromoethane ethyl acetate benzene 1-bromobutane 1-chlorobutane methanol toluene ethanol 1-chlorohexane propanol ethylbenzene p-xylene m-xylene n-butanol cumene 1,3,5-trimethylbenzene ether isoprene dibutyl ether cyclohexene cyclohexane hexane

logkexp logkpred predicted ranking -1.74 -1.84 -2.00 -2.12 -2.15 -2.20 -2.22 -2.23 -2.26 -2.26 -2.26 -2.32 -2.33 -2.44 -2.46 -2.50 -2.50 -2.57 -2.60 -2.66 -2.66 -2.66 -2.78 -2.79 -2.83 -2.83 -2.90 -2.91 -2.93 -3.01 -3.05 -3.08 -3.15 -3.21 -3.21 -3.23 -3.32 -3.40 -3.44 -3.52 -3.66 -3.66 -3.66 -3.80 -3.80 -3.89 -3.91 -3.99 -4.04 -4.06 -4.11 -4.12 -4.40 -4.70 -4.93 -5.18 -5.48 -5.93 -6.78

-0.94 -1.86 -2.05 -2.63 -1.96 -3.53 -1.94 -1.21 -1.84 -1.22 -2.66 -2.61 -2.52 -2.79 -2.49 -2.87 -3.53 -2.80 -3.45 -3.08 -3.50 -1.24 -3.48 -3.47 -3.43 -3.16 -3.73 -4.64 -3.74 -2.08 -3.97 -3.41 -3.60 -3.34 -3.48 -4.36 -4.07 -4.61 -3.90 -3.99 -4.64 -4.59 -3.91 -3.96 -3.78 -4.62 -3.70 -4.10 -4.03 -4.03 -3.72 -4.19 -3.89 -4.92 -5.18 -5.01 -4.86 -4.90 -5.74

1 6 9 14 8 29 7 2 5 3 15 13 12 16 11 18 30 17 24 19 28 4 26 25 23 20 34 52 35 10 41 22 31 21 27 48 45 50 38 42 53 49 39 40 36 51 32 46 43 44 33 47 37 56 58 57 54 55 59

denotes the solvent with the smallest rate constant. This can be used for the qualitative analysis of the results obtained from the regression. Benzyl cyanide is ranked first by prediction and by experiment. If we consider the first 10 experimentally tested solvents, 7 are among the top 10 solvents given by the model predictions. For the first 20 experimentally tested solvents, 17 are predicted to be within the top 20. Therefore, the model provides a ranking suitable for solvent design. Overall, given the statistics and rankings presented for the regression performed, the solvatochromic equation seems to give

solvent

EN T value

benzyl cyanide 1,1,2,2-tetrachloroethane nitropropane bromobenzene tetrahydrofuran benzene ethanol hexane

0.37 0.27 0.37 0.18 0.21 0.11 0.65 0.01

Table 5. Ranked List of Top 10 Solvents Obtained in the Deterministic Design Step for the Menschutkin Reaction Solvent rank

name

objective function

1 2 3 4 5 6 7 8 9 10

2-methylene-1,8-dinitrooctane 2-methylene-1,7-dinitroheptane 2-methylene-1,6-dinitrohexane 2-methylene-1,5-dinitropentane 2-methylene-1,4-dinitrobutane (3-nitro-2-butenyl)aniline 2-methylene-3-methyl-1,8-dinitrooctane 2-methylene-3-methyl-1,7-dinitroheptane 2-methylene-3-methyl-1,6-dinitrohexane 2-methylene-3-methyl-1,5-dinitropentane

0.94 0.93 0.92 0.91 0.88 0.77 0.76 0.75 0.75 0.74

good predictions of the effect of solvent on the logarithm of the reaction rate constant for the Menschutkin reaction, although a small set of solvents is used for regression. 6.2. Deterministic Design of Solvent Molecules. Having developed a model of solvent effects for the Menschutkin reaction, we now focus on generating new solvent candidates. Because there is only one reaction, the following simple objective function is used: max logk (19) where logk is given by eq 18. Solving the design problem, which involves the constraints given in Appendix A, we obtain a ranked list of candidate solvents that give high values for the reaction rate constant for the Menschutkin reaction. In Table 5, a list of the 10 best ranked solvents and their objective function values is shown. The three best solvents are presented in Figure 2. Nine out of the 10 candidate solvents in Table 5 are unsaturated aliphatic dinitrates, which points to this chemical class as a source of good solvents for the Menschutkin reaction. The unsaturated bond on C2 carbon in most of the structures reported in Table 5 may readily lead to isomerization into conjugation with the N1 nitrogen, forming a CdCH-NO2 group, which would react with nucleophiles. These compounds are thus not standardly used as solvents, and this explains their absence from the experimental list. As an alternative, saturated aliphatic dinitrates (1,n-dinitroalkanes, where n ranges from 2 to 9) are also identified as good candidates. Thus, 1,9dinitrononane has an objective function value of 0.64 and is predicted to be the best unsaturated molecule, followed closely by similar compounds. No dinitro compounds were present in the experimental lists given in Table 3; therefore, their performance must be tested in the laboratory for these results to be confirmed experimentally. Most of the other highly ranked solvent candidates contain at least one nitro group. They include, for instance, benzene derivatives, with a nitro group and another functional group such as the amine group. Compounds with a single NO2 group (e.g., nitroethane and nitropropane) are ranked high in the list of experimentally tested solvent shown in Table

Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008 5197

Figure 3. Number of sampled points versus the number of generated designs for the two-dimensional uncertainty space (s-h) for the Menschutkin reaction. Table 8. Results for the Two-Dimensional Sensitivity Analysis Runs for the Menschutkin Reaction

Figure 2. Structures of the three best molecules generated in the design step for the Menschutkin reaction. Table 6. Lower and Upper Bounds of the 95% Confidence Intervals for the Reaction Parameters for the Menschutkin Reaction in eq 18 95% Confidence Interval parameter

lower bound

upper bound

s d a b h

-3.46 -5.07 -13.54 -12.31 -11.34

11.07 5.02 19.87 12.48 10.47

Table 7. Sampled Parameter Range and Number of Distinct Designs from One-Dimensional Sensitivity Runs for All Five Parameters for the Menschutkin Reaction parameter varied

range

number of designs generated

sensitivity index, Sθ

s d a b h

14.53 10.09 33.41 24.79 21.81

9 2 3 5 11

0.62 0.20 0.09 0.20 0.50

3, which indicates that the properties of the nitro group affect this reaction favorably. The molecules designed should be considered with caution, although they give a good indication of what would be the best type of solvents to target for experimental verification. Because the reaction model is built based on kinetic data in eight solvents only, and then extrapolated in the design step to predict rate constants for nearly 10 000 compounds, high uncertainty in the rate constants can be expected. The bounds (upper and lower bound) of the 95% confidence intervals for the parameters in eq 18 are shown in Table 6. For example, the most uncertain parameter is a, for which the lower bound is -13.54 and the upper bound is 19.87. This wide range

parameters varied

size of parameter space

number of molecules generated

s and d a and b s and a s and b s and h

14.53 × 10.09 33.41 × 24.79 14.53 × 33.41 14.53 × 24.79 14.53 × 21.81

13 9 28 27 34

reflects a fairly high uncertainty. Robust molecular design is addressed in the next section through the strategy presented in section 5. 6.3. Design of Solvent Molecules under Uncertainty. In this section, the reliability of the reaction model developed for the Menschutkin reaction in section 6.1 is investigated with global sensitivity analysis to determine the key parameters and obtain a “‘map of solutions”. Based on these results, a small number of representative scenarios are identified. For those scenarios, we then solve a stochastic optimization problem and compare the resulting robust solutions with those obtained via deterministic optimization in section 6.2. Sensitivity analyses are performed in one, two and five dimensions. The results of a two-dimensional run and of the five-dimensional run are used to formulate and solve the problem of solvent design under uncertainty. 6.3.1. One-Dimensional Sensitivity Analysis. For the onedimensional sensitivity analysis, 1000 parameter combinations were sampled. No new designs were generated in the last 500 combinations. The numbers of designs generated from these runs, together with the width of the uncertainty range for each parameter, are shown in Table 7. It can be seen from the Sθ values that the design is almost insensitive to parameter a, more sensitive to parameters b and d, and notably more sensitive to parameters s and h. This is a similar result to that obtained for the solvolysis reaction.8 Here, we have a case of a reaction in which the product is an electrolyte, so the reaction involves the formation of ions and, most likely, the solvent dipolarity/ polarizability and cohesive energy density contribute the most to the value of logk. 6.3.2. Two-Dimensional Design under Uncertainty. 6.3.2.1. Sensitivity Analysis. To illustrate the proposed procedure graphically, two-dimensional sensitivity analyses are performed. The Sobol’ sequence24,25 is used to sample the uncertainty space, starting with 1024 (210) sampled points, and increasing the number in each run by 1024 points. A graph showing the change in the number of designs with the number of sampled points for the two-dimensional case where convergence was achieved with the slowest rate is presented in Figure

5198 Ind. Eng. Chem. Res., Vol. 47, No. 15, 2008

Figure 4. Solution map for s vs h, divided into clustered areas of different design solutions. All molecules generated are presented. The largest areas correspond to the following molecules: 1, glycerol, 2, 2-methylene-1,8-dinitrooctane, 3, 6,7-dimethyl-8-nitro-2-nitromethyl-1-octene; and 4, iso-butane.

3. This figure shows that convergence is achieved after 6144 points. The results of five two-dimensional runs are shown in Table 8. As expected, varying two parameters at a time results in an increase in the number of designs generated. If the two parameters varied describe different interactions in the system, more than twice the number of designs are obtained than when one combines parameters that describe similar types of interactions (s and d or a and b). The largest number of designs is obtained when varying simultaneously parameters s and h. The design is almost insensitive to parameter a, more sensitive to parameters b and d, and notably more sensitive to parameters s and h. This is consistent with the results obtained in the onedimensional analysis. 6.3.2.2. The Design Problem. Having identified s and h as the key parameters, we focus only on their combination to obtain a robust solution in the two-dimensional uncertainty space, while keeping the other parameters at their nominal values. The map of clustered solutions in the two-dimensional (s and h) space is shown in Figure 4. Ninety representative scenarios are obtained by dividing each cluster into three subareas. Out of 34 molecular design clusters, 4 molecules appear in too few runs (