A hierarchical method to integrated solvent and process design of physical CO2 absorption using the SAFT‐γ Mie approach

Molecular-level decisions are increasingly recognized as an integral part of process design. Finding the optimal process performance requires the integrated optimization of process and solvent chemical structure, leading to a challenging mixed-integer nonlinear programming (MINLP) problem. The formulation of such problems when using a group contribution version of the statistical associating fluid theory, SAFT-γ Mie, to predict the physical properties of the relevant mixtures reliably over process conditions is presented. To solve the challenging MINLP, a novel hierarchical methodology for integrated process and solvent design (hierarchical optimization) is presented. Reduced models of the process units are developed and used to generate a set of initial guesses for the MINLP solution. The methodology is applied to the design of a physical absorption process to separate carbon dioxide from methane, using a broad selection of ethers as the molecular design space. The solvents with best process performance are found to be poly(oxymethylene)dimethylethers. © 2015 American Institute of Chemical Engineers AIChE J, 61: 3249–3269, 2015


Introduction
The choice of solvents or other functional chemicals greatly influences the performance of processes. These decisions are often made prior to the design of the flow-sheet structure and operating conditions on the basis of experience, heuristic rules, or a few exploratory experiments. A more systematic approach in selecting processing materials is computer-aided molecular design (CAMD), which aims to find the molecular chemical structure(s) whose properties match some user-specified target properties (cf. the overviews in Refs. 1 and 2). Although this approach can lead to improved molecular decisions, it is generally difficult to map the economics of a process to a few "key properties" of the desired molecules. For example, Kossack et al. 3 used the selectivity of an entrainer as an objective to represent the cost of an extractive distillation but found that this objective leads to an unfavorable entrainer choice. Papadopoulos and co-workers 4,5 proposed the use of multiobjective optimization with several key properties of the solvent as objectives. A Pareto-optimal set (a set of solutions representing best compromises of the objectives) of solvents was obtained, which served as a reduced solvent space for a subsequent optimization including the process model, in order to arrive at a better solution from a process-wide perspective.
The extent to which the discrete decisions relating to the choice of solvents or other additives are linked to choices that relate to process operation and process design can be further understood by considering the case of a solvent which is liquid only over a limited range of temperature. The choice of such a solvent as a reaction medium constrains the operating temperature that can be achieved in the reactor, and therefore limits reaction rates and productivity. Furthermore, the "best" values of the physical properties are intimately linked to process design and economic targets, especially when conflicting requirements arise. For example, in an absorption process, the solvent must offer high capacity and selectivity for the species to be absorbed, yet it must be easily separated. In the context of reaction-separation processes, yet more complex trade-offs are often necessary. [6][7][8] The integration of solvent and process design is, therefore, highly desirable to uncover better solutions. 2,9,10 Computer-aided molecular and process design (CAMPD), that is, the simultaneous optimization of a molecule's chemical structure and of the process structure and/or conditions, usually leads to complex mixed-integer non-linear Correspondence concerning this article should be addressed to J. Burger at jakob.burger@mv.uni-kl.de or C. S. Adjiman at c.adjiman@imperial.ac.uk.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. programming (MINLP) problems. To solve these problems using standard MINLP solvers is prohibitive in most cases of practical interest due to the relatively large number of discrete variables and strong non-linearities in the process models, which give rise to multiple local solutions and cause significant difficulties in initializing the non-linear problems (NLPs) that must be solved in the course of the MINLP algorithm. The development of CAMPD-specific solution strategies for MINLP problems of this type is, therefore, a challenging research field that has attracted increasing attention.
Hamad and El-Halwagi 11 first introduced the CAMPD problem and tackled it by optimizing linearized process models together with the solvent structure. Other approaches based on the solution of a simplified CAMPD problem have focused on reducing the generally large space of flow-sheet options and solvent candidates. Hostrup et al. 12 formulated several heuristic scenarios to eliminate process flow-sheet options and solvents before solving the reduced CAMPD problem numerically. Karunanithi et al. 13 reduced the discrete solvent space before explicitly considering the process. Chemically infeasible solvents were excluded, and constraints imposed by process operating conditions were used to eliminate further solvents based on the evaluation of the pure-component or binary-mixture properties.
Several solution strategies based on decomposition of the complex optimization problem have also been developed. Buxton et al. 14 suggested a decomposition-based algorithm in order to enable the solution of the CAMPD problem for non-linear process models. Two subproblems were solved iteratively in this approach: a primal NLP with fixed solvent structure to optimize the process conditions and a master mixed-integer linear problem to optimize the solvent chemical structure using a linearized process model. The authors noted that if the primal problem is solved in the standard way, the NLP algorithm often fails to converge due to the non-linearity of the process model. 14 A primal problem solution strategy consisting of several preprocessing steps was proposed, in which the solvent candidate at the current iteration is assessed by means of physical property tests. If any of these tests fail, an alternative reduced master problem is solved to reinitialize the solvent molecule. Giovanoglou et al. 15 extended this strategy to dynamic processes.
An alternative approach to address the complexity of the problem is to use stochastic algorithms as proposed by Marcoulaki and Kokossis, 16 who used a strategy based on simulated annealing 17 to solve the integrated solvent and process design problem. Stochastic approaches are robust in that they easily recover from the failed evaluations of the model that arise from non-linearity and infeasibility, and they are well suited to the presence of discrete optimization variables. However, these types of algorithms require many process evaluations to reach the optimal solutions, and thus, a high computational effort.
A more recent decomposition strategy is exemplified by the work of Eden at al., 18 Eljack et al., 19 and Bommareddy et al. 20 who formulated an alternative optimization problem in which the aim is to find the optimal physical properties that minimize the process costs without considering the solvent chemical structure explicitly. Once the optimal values of the properties have been identified, the chemical structures (pure components and mixtures) that possess these properties must be found by solving a separate CAMD problem. Another twostage approach, continuous molecular targeting (CoMT), was proposed by Bardow et al. 9 In the first stage, they described the solvent in their absorption process as a virtual component characterized with parameters based on the perturbed-chain statistical associating fluid theory (PC-SAFT) equation of state (EoS). 21 The parameters of the physical property model were treated as continuous optimization variables during the optimization of the process, formulated as an NLP. In the second stage, chemical structures were matched to identified model parameters. The CoMT approach has now also been applied to find better working fluids for organic Rankine cycles. 22,23 The challenge in both two-stage approaches 9,[18][19][20] is the second stage, as it is not guaranteed that the optimal molecule for the CAMD subproblem is optimal from the perspective of the overall problem.
In another recent example of CAMPD work, the statistical associating fluid theory for potentials of variable range (SAFT-VR) 24,25 has been used to optimize simultaneously the operating conditions of a physical absorption process for the high-pressure separation of carbon dioxide (CO 2 ) from methane (CH 4 ) and the molecular structure (number of carbon atoms) and composition of an n-alkane blend used as the solvent. 26,27 The number of carbon atoms characterizing the n-alkane was treated as a continuous parameter in the CAMPD resulting in an NLP. By invoking the principle of congruence, non integer values for the optimal chain length can be interpreted as representing a mixture of two n-alkanes with chain lengths differing by one carbon such that the average chain length in the mixture matches the optimal value.
While there have been encouraging advances in CAMPD techniques, the frequent use of a restricted design space and/or decomposition approaches in order to reduce the problem complexity can lead to suboptimal solutions. At the level of problem formulation, a particular challenge in CAMD and CAMPD is to identify a physical property model that is able to represent all solvent candidates and all relevant properties well. 28 The required physical properties are often predicted by group contribution (GC) methods, in which the solvent molecule is divided into a small number of chemical functional groups, such as hydroxyl OH, methyl CH 3 , methanediyl CH 2 ; for a comprehensive review on GC methods for use in design see Papaioannou et al. 29 By necessity, many authors use a separate GC method for every physical property, for example, vapor pressure, activity coefficients, enthalpy of vaporization, but the description of these different properties may be inconsistent as a result. For example, the liquid heat capacity of a substance can be calculated from the vapor heat capacity and the enthalpy of vaporization, or it can alternatively be modeled independently to increase its accuracy; modeling these three properties independently may be thermodynamically inconsistent. Moreover, most property-specific GC methods are restricted to pure-component properties, imposing the rather crude assumption that all mixtures considered exhibit ideal mixing. Finally, thermodynamic consistency across multiple phases is particularly important when one or more components are at near-critical conditions, and this cannot be achieved using different models for the vapor and liquid phases.
EoSs based on the GC formalism can help to overcome this inconsistency by providing a thermodynamic model from which all thermophysical properties (with the exclusion of transport properties) of all fluid phases can be determined consistently, without any assumptions of ideality. Advanced molecular EoSs such as the statistical associating fluid theory (SAFT) 30,31 are well suited to treating highly non-ideal, asymmetric, and associating systems commonly encountered in separation processes, and can be used to provide an accurate representation of liquid and vapor phases. These characteristics give SAFT EoSs a wide range of applicability. 32-35 GC methods have been proposed for several versions of SAFT 29,36 and can provide a consistent physical property description for use in CAMD and CAMPD. Recently, Perdomo et al. 37 used the SAFT-c GC EoS based on the square-well potential [38][39][40] in a CAMD study to find biofuels with improved performance.
Building on the advances that have been made in CAMPD and GC techniques in recent years, three distinct aims are considered in our article. First, we demonstrate that a GC version of SAFT can be used to solve CAMPD problems successfully on a thermodynamically consistent basis. In our work, the recently developed SAFT-c Mie GC EoS 41 based on the versatile Mie potential form 42 is used. Parameters for groups not previously reported are developed to expand the solvent design space.
As a second key objective, a novel hierarchical optimization (HiOpt) approach is proposed in order to overcome the complexity of the integrated solvent and process optimization problem. A reduced process model that satisfies only mass balance and thermodynamic constraints is used in an initial multiobjective optimization (MOO) step to determine a Pareto set, that is, a set of non-inferior solvents and process conditions. The Pareto set provides initial guesses for a second optimization step in which the original, full, model is used. The initial guesses for the solvent structure and the process conditions enable the solution of the challenging CAMPD problem.
The third and final aim of our work is to demonstrate the integrated use of the SAFT-c Mie model and of the hierarchical CAMPD approach to identify novel solvents and the corresponding process conditions for the separation of CO 2 from natural gas (principally containing CH 4 as the other component) by physical absorption. Given the large range of pressures characteristic of this process, which includes the critical pressure of CO 2 (73 bar), the use of an EoS that can model these conditions reliably is essential. A novel solvent is proposed for the process, which is predicted to lead to higher overall economic performance than conventional solvents.
The article is organized as follows. The proposed CAMPD methodology is described first, followed by the application of interest, namely the integrated design of a CO 2 absorption process and solvent; the simplified and detailed process models used for hierarchical decomposition are presented alongside the problem statement. Next, a description of the thermodynamic properties obtained with the SAFT-c Mie platform is presented and new group interaction parameters are developed to represent the desired solvent design space. Finally, the results of the application of the hierarchical CAMPD approach to the CO 2 absorption process are discussed before conclusions and perspectives are provided.

HiOpt Approach to CAMPD Methodology
We propose a novel hierarchical optimization approach to solve CAMPD problems, which is referred to as HiOpt. Systematic hierarchical approaches for the design and optimization of complex chemical processes have long been used to make these problems tractable. 43,44 The design is structured in several stages of decisions in which the complexity of the process is increased in a stepwise manner. Hierarchical decision systems have found a broad range of applications, as seen for instance in the work of Franke et al., 45 Han et al., 46 or Burger et al. 47 A particular challenge in hierarchical design is to avoid eliminating good options too early in the design procedure. In order to address this issue, Burger and Hasse 48 proposed a design method that allows one to keep some flexibility when moving to the next level. Instead of focusing on one specific option, a Pareto-optimal set which contains several initial process options provides starting points for the design at the next level.
In our work, the methodology of Burger and Hasse 48 is adapted to solve the MINLP problem arising in CAMPD, which is given by min x;y f ðx; yÞ where the w-dimensional vector x consists of the continuous variables, typically representing operating conditions, equipment sizes, and physical properties of the process. The q-dimensional vector y of integer variables describes the process topology and the molecular chemical structure of the solvent. X and Y are bounded subsets of R w and N q , respectively. The design space of all possible combinations of x and y is given by X3Y. f is the objective function, typically an economic performance metric. The process model consists of the equations involving h p and h t representing process unit models and thermodynamic property models, respectively. The vector g represents constraints resulting from operability limits. Molecular feasibility and complexity rules are represented by the constraints h m and g m . As the process model and the objective function f are generally nonlinear, problem (CAMPD) is often very challenging to solve. In general, even with a local solver, a solution is only found if good initial guesses for the optimization variables x and y are specified. The process unit models denoted by h p ðx; yÞ ¼ 0 in problem (CAMPD) are referred to as the "full" process unit models. In addition to the full process unit models, reduced process unit models are required to develop the methodology. They preferably have fewer specifications and fewer, less complex, equations than the full models. These reduced process unit models are combined with the thermodynamic model (subset h t ðx; yÞ ¼ 0) to form a reduced process model. Using the reduced process model, a second optimization problem is formulated: min z;y f 1 1 ðz; yÞ; f 1 2 ðz; yÞ; f 1 3 ðz; yÞ; ::: complexity of the reduced process unit models denoted by h 1 p ðz; yÞ ¼ 0, which have fewer state variables and fewer degrees of freedom. The thermodynamic model (h 0 t ðz; yÞ ¼ 0) and the constraints (g 0 ðz; yÞ 0) are adapted to the reduced number of continuous variables z, leading to a reduction in the number of equations in these submodels. The molecular complexity and feasibility constraints, h m ðyÞ ¼ 0 and g m ðyÞ 0, are unchanged as they do not involve continuous variables. The reduced process unit models typically do not provide enough information to calculate the original objective function f. For example, it is not possible to size the process units based on the reduced information, and hence an estimate of capital costs is not available. Therefore, the problem is formulated as a MOO problem. Several surrogate objectives f 1 i represent different contributions to the cost of the process. Problem (MOO) does not have a unique optimal solution, but multiple non-inferior solutions, defining the Pareto-optimal set.
The choice of the objectives f 1 i depends strongly on the process studied and should be made carefully not to omit crucial cost contributions inadvertently. If the objectives are chosen appropriately, the solution of MOO can provide valuable information on the interactions between the process and the solvent. Generally, a large part of the design space is inferior, that is, not included in the Pareto-optimal set.
The proposed two-step HiOpt methodology, shown in Figure 1, incorporates problems (CAMPD) and (MOO). It is assumed that the objective function, the design space and all of the constraints in both problems are specified. In the first step, the Pareto-optimal set of problem (MOO) is determined, based on an approximation by a finite number of Pareto-optimal solutions. In the second step, these solutions provide initial guesses for the solution of problem (CAMPD). As a heuristic, molecules and process topologies that are represented by inferior values of y can be excluded from consideration in problem (CAMPD) to simplify the solution of the MINLPs in the final step. In a similar manner, tighter bounds for the continuous variables x can be chosen based on the designs identified in the solution of the MOO. In general, there is, however, no guarantee that the Pareto-optimal set contains the molecular structure and the process configuration that correspond to the global solution of problem (CAMPD).
The HiOpt approach can be understood as a generalization of the approach presented by Papadopoulos and Linke. 4 Instead of using a reduced process model, they omitted the process model completely in their approach. The objectives of their MOO problem were based on the physical properties of the solvent. Papadopoulos and Linke showed, for some simple processes, that the set of Pareto-optimal solvents obtained with their formulation of (MOO) includes the solution of the corresponding problem (CAMPD). The introduction of a reduced process model, as proposed in our work, enables the designer to choose objectives that represent the process cost in a more direct way. Process conditions, such as the operating pressures and temperatures, can be used in the objective functions and constraints of problem (MOO). The risk of omitting crucial cost contributions is thus reduced. In addition, the Pareto-optimal set of solvents may be smaller (compared to the case when only properties are considered as objectives) due to the constraints imposed by the reduced process model. Furthermore, the final solution of problem (CAMPD) becomes numerically less challenging as initial guesses for the process conditions are available.

Reduced process unit models
In our current work, reduced process unit models are used to calculate the values of variables characterizing the streams leaving the process units with respect only to mass balances and thermodynamic limits. Hence, the streams leaving the reduced process units are physically meaningful, even though the idealized treatment may incur high investment and operating cost. For distillation systems, these types of models are known as 1/1-models. [49][50][51][52] In such a model, a distillation column is assumed to be of infinite height and to operate under total reflux. Under these assumptions, the distillation products can be recovered to arbitrary purity as long as azeotropic and distillation boundaries, the thermodynamic limits, are respected. The 1/1-models of distillation columns can be represented with linear equations which can be solved efficiently and robustly. 52 In contrast to the infinitely large dimensions and energy requirements of the distillation columns represented by such models, the resulting product streams are often very close to those in realistic process models. 47,[53][54][55][56] Quantitative performance metrics can, therefore, be derived on the basis of the stream variables and can serve as the multiple objectives f 1 i in problem (MOO). The reduced process models require a much smaller number of physical property calculations than the full model but result in a feasible, or nearly feasible, and physically meaningful process stream table.

Mixed-integer MOO
There exist efficient deterministic algorithms to approximate the Pareto set of a MOO problem that involves only continuous variables, even if the Pareto set is non-convex. 57 By contrast, mixed-integer MOO (also called combinatorial MOO) is a rather new field of research in which a mainly heuristic treatment is used to approximate the Pareto-optimal set. 58 Fattahi and Turkay 59 recently proposed a deterministic algorithm based on the -constraint approach to calculate Pareto points in mixed-integer problems. The approach is, however, computationally expensive if the complete Pareto set has to be sampled.
In our work, a modification of the deterministic sandwiching algorithm as proposed by Bortz et al. 57 is used as a heuristic to achieve a compromise between performance and sampling quality. This algorithm was originally developed for continuous convex problems and is used to approximate the Pareto-optimal set with only a small number of Pareto points. Inner and outer approximations for the Pareto set are obtained iteratively, with increasing accuracy. The distance between the two approximations is used implicitly, via a ratio of distances, to test for convergence within a userdefined tolerance d S (a scalar that is greater than 1, with 1 corresponding to an exact match of the Pareto set). To calculate Pareto-optimal points, the algorithm makes use of the concept of weighted-sum scalarization, that is, the multiple objectives in problem (MOO) are replaced by their weighted sum f R , leading to a single-objective optimization problem that yields one Pareto-optimal point where N f denotes the number of objectives in the problem and w i is a positive number denoting the weight of objective f 1 i . The progress of the algorithm, including the selection of the weights w i is sketched in Figure 2 for the case of two objectives and is summarized as follows: 1. Find the extreme Pareto points using the unit vectors ð1; 0Þ T ; ð0; 1Þ T as weights and optimize the weighted sum. Construct the nadir point* N, cf. Figure 2a.
2. An inner approximation of the Pareto set is given by the line segment joining the two extreme Pareto points as shown in Figure 2a. The outer approximation is given by two line segments, each passing through one Pareto point and perpendicular to the weight vector associated with that Pareto point. The two segments intersect at point O. Point I is the intersection of the inner approximation and the line connecting point O and the nadir point N. If the Pareto set is convex, all Pareto-optimal solutions are located between the inner and outer approximations.
3. Two neighboring Pareto points form a facet. The quality of the approximation r S is determined for every facet. To do so, the inner and outer distances are calculated as the distances between point I and the nadir point N and point O and the nadir point N, respectively, cf. Figure 2a. The quality of approximation r S is defined as the ratio of the outer distance to the inner distance: r S ¼ ON IN : If the two approximations are the same, then r S is equal to 1. Otherwise, it is greater than 1. By construction, r S ! 1.
4. For the facet with maximal r S , refine the approximation: use the vector normal to the inner approximation as the weight vector w and optimize f R ðz; y; wÞ. If this yields a new Pareto point, replace the facet with two new facets around the new point, cf. panels b) and c) in Figure 2. If the optimization yields a previously known Pareto point, the outer approximation equals the inner approximation for this facet: 5. If the maximum of all r S is larger than the convergence tolerance d S , then go to step 3, otherwise stop.
The algorithm can also be applied to problems with more than two objectives; the reader is referred to Ref. 57 for a more general description. As weighted sum scalarization is used in the algorithm, only the convex hull of the Paretooptimal set is determined. Non-convex parts of the Pareto set, for example, point M in Figure 2d, are not found. In the absence of an algorithm to determine efficiently non-convex Pareto points, it is assumed that the convex approximation of the Pareto set contains a sufficient number of initial guesses to find a good solution of problem (CAMPD).

The SAFT-c Mie EoS
A physical absorption process to separate CO 2 from CH 4 is studied. To obtain meaningful results, the thermodynamic property model must represent, accurately, the vapor-liquid equilibrium (VLE) of CO 2 , CH 4 , and the solvent of choice at high pressures and moderate temperatures. The SAFT family of equations has been used successfully to model such mixtures. Pereira et al. 27 studied the same process and successfully used a SAFT EoS for potentials of variable ranges (SAFT-VR) based on the square-well potential model 24,25 to represent the properties of both phases in the VLE region over a wide range of pressures, temperatures, and composition, when using n-alkanes as solvents. Nannan et al. 60 have proposed a model for use with the PC-SAFT EoS 21 for poly(oxyethylene)dimethylethers as solvents for the description of high-pressure CO 2 absorption processes. Qadir et al. 61 have also used a similar version of PC-SAFT for the description of another absorption process to separate CO 2 from hydrogen.
For the CAMPD problem considered here, the SAFT-c Mie GC EoS based on the Mie potential form 41 is used. In contrast to the SAFT versions considered so far in studies of such processes, SAFT-c Mie is a GC formulation and is thus ideally suited for the representation of a large set of solvents in the process. It can be used to describe the properties of the pure component as well as the mixtures for both the vapor and liquid phases on an equal footing. For the VLE of mixtures of CO 2 , CH 4 , and n-alkanes, it has been shown to give good agreement with experimental data. 62 The extension of the group parameter base to represent a large set of solvents is discussed later. Some basic characteristics of the SAFT-c Mie EoS are presented first for completeness.
In the SAFT-c Mie EoS, the molecules are described in terms of interaction potential parameters of the segments characterizing the various functional groups. Each molecule is represented by N k groups of type k and each group k consists of one or more identical segments as defined by the parameter m Ã k . Each segment is further characterized by a shape factor, S k ; the shape factor characterizes the extent to which the group's segments contribute to the overall thermodynamic properties of the molecule. The interactions between two segments are modeled using the Mie (generalized Lennard-Jones) potential form. The Mie pair interaction energy U kl of two segments of types k and l is given as a function of the distance r between the centers of the segments by 41,63 where kl is the depth of the potential well, r kl is the segment diameter, k r kl and k a kl are the repulsive and attractive exponents of the potential. The prefactor C kl is given by so that the minimum of the interaction energy is 2 kl . Six parameters, therefore, describe a (non-associating) functional group k: the number of segments m Ã k , the potential parameters r kk , k a kk ; k r kk , and kk , and the shape factor S k . This is augmented by unlike parameters to describe the interactions between two different types of segment k and l. The unlike parameters can be calculated using the following combining rules 41 In some cases, the parameters kl , k a kl , and k r kl are adjusted to best describe the target experimental data. In our work, the unlike dispersion energy parameter kl is either obtained from the relevant combining rule or estimated from experimental data, as described later. Other unlike interactions are obtained from the combining rules given in Eq. 4. In the SAFT-c Mie approach, it is also possible for the groups to have association sites that mediate strong directional interactions such as hydrogen bonds, 64,65 but this feature is not used in our current work. For further details of the SAFT-c Mie EoS, the reader is referred to Ref. 41. Description of the Case Study CO 2 absorption from natural gas A prototypical process for the removal of CO 2 from natural gas is studied. Natural gas is usually produced at high pressure and often has a high concentration of CO 2 (30 to 50 mol %); in extreme cases, it may have a CO 2 concentration of 70 mol %. 27 The CO 2 may be of natural origin or may stem from the injection of CO 2 into near-depleted reservoirs as part of enhanced oil recovery. 27 In order to introduce natural gas into the distribution system, for example, pipelines, the fraction of CO 2 must be reduced to below 3 mol %. 66 As the CO 2 concentration in the natural gas feed can be appreciable, removing it by chemical absorption with heat regeneration of the solvent is expensive. Physical absorption processes in which CO 2 is stripped from the loaded non-reactive solvent merely by reducing the pressure, without the application of heat, becomes a favorable option. 67 In such a process, the solvent should exhibit a high capacity for CO 2 , a high selectivity of CO 2 compared to natural gas and a low vapor pressure (to minimize solvent losses). Achieving these targets depends not only on the solvent's molecular structure but also on the process conditions. The optimal solvent can, therefore, only be identified with confidence by considering the integrated solvent and process design problem.
In our current work, we adopt the problem statement, the process flow sheet, and the process unit models of Ref. 27. The process flow sheet is given in Figure 3. The natural gas feed, Stream (1), is considered to be a binary mixture of CH 4 and CO 2 at high pressure. It is expanded and fed to a counter-current absorption column. The solvent absorbs CO 2 preferentially and a cleaned gaseous stream leaves the column at the top of the absorber as Stream (3). The pressure of the CO 2 -rich solvent leaving the bottom of the absorber, Stream (4), is then reduced to regenerate the solvent. The resulting CO 2 -rich vapor phase is separated from the liquid solvent in a flash drum and purged from the process via Stream (6). After regeneration, the lean solvent is pumped back to the absorber. To compensate for solvent losses in the gaseous product streams, a small feed, Stream (8), of fresh solvent is also supplied.
The specifications of the process are mainly adopted from Ref. 27; the full process model divided into fixed variables and optimization variables is specified in Table 1. The conditions of the natural gas feed and the fresh solvent feed are considered to be given and are not optimized. In the work of Ref. 27, it was found that low pressures in the flash drum lead to lower process costs. As the CO 2 outlet gas (Stream (6)) should have a minimum pressure of 0.1 MPa (ambient pressure), this low pressure is used for the flash drum. The number of theoretical stages in the absorption column is set to 10. The problem studied here can thus be understood as a retrofit design problem in which the number of stages is given by an existing absorption column. The absorber pressure and the solvent recycle flow rate are the continuous optimization variables, belonging to x in problem (CAMPD). Loose bounds are chosen for the optimization variables. They represent, however, values that are assumed to be achievable with standard process equipment. Additional constraints are imposed to increase the likelihood of obtaining a realistic process 27 : the cross-sectional area A ABS of the absorption column must be less than 30 m 2 .
the dynamic viscosity l of the solvent at every point in the process must be less than 100 cP. the temperature at every point in the process must be at least 10 K above the melting temperature T melt of the solvent. the mole fraction x ð3Þ CH4 of CH 4 in the purified gas must be at least 0.97 mol/mol.
The absorption column is modeled as isobaric using the equilibrium stage model. The flash drum is modeled as an adiabatic equilibrium stage. The pump is assumed to operate isentropically, and the expansion valves are modeled as adiabatic. The process unit models and detailed costing functions to calculate the net present value (NPV) of the process are given in Ref. 27.
The process model adopted from Ref. 27 is referred to as the full model in the remainder of our article. The intention of our work is to demonstrate the proposed hierarchical approach and to identify promising solvent candidates. Therefore, the flow-sheet structure is not optimized and is considered fixed in this study. It is assumed that the process flow sheet presented in Figure 3, despite its simplicity, is suitable to discriminate clearly between different solvent structures according to their potential for CO 2 removal.
As different solvents may lead to optimal performances under different process conditions, the pressure in the absorp-tion column, and the solvent recycle rate are treated as optimization variables. All other process variables are fully determined once these degrees of freedom are fixed. The problem studied can be concisely stated as: Given a feed stream of CO 2 and CH 4 (composition and state), find the solvent chemical structure, absorber pressure, and solvent recycle rate that maximise the net present value of the process, while meeting all constraints.

Reduced process model
In the full model, the absorption column is represented by one set of MESH equations for every stage. These MESH equations include for every stage: conditions for the phase equilibrium, mass balances, and an energy balance. In addition to the complexity inherent in this large number of nonlinear equations, the presence of the solvent recycle further hampers the convergence of the process simulation. The simulation is unlikely to converge without good initial guesses for the variables of the defined system of equations. The optimization problem (CAMPD), which invokes the process model numerous times, is even harder to solve.
To overcome this complexity, reduced process unit models for the absorption column and the flash drum are proposed. They are used in the multiobjective problem (MOO) in the first step of the HiOpt approach. The reduced models are based on the following assumptions: a. the absorption column is of infinite height; b. the solvent loss at the top of the absorption column is negligible compared to the solvent loss in the flash drum; c. the concentration of solutes (CO 2 , CH 4 ) in the lean solvent after the regeneration in the flash drum is negligible; d. the vapor phase in the flash drum is ideal; e. the temperature is constant throughout the process. As a consequence of Assumption (a), the maximum separation efficiency is achieved in the absorption column, i.e., the number of stages is infinite.
Assumptions (b) and (c) are reasonable if the pressure in the flash drum is very low compared to that in the absorption column. As seen in Table 1, the pressure in the flash drum is likely to be at least one order or magnitude smaller than that in the absorption column. If the temperatures and flow rates of both gaseous streams leaving the process are of the same order, then the solvent losses in the CH 4 product stream is expected to be small compared to the losses in the flash drum. Furthermore, the gas solubility in the flash drum is drastically reduced by the large pressure drop, leading to a very small concentration of solutes in the solvent leaving the flash drum. Assumption (d) is reasonable for low pressures. The pressure in the flash drum of 0.1 MPa is considered to be sufficiently low.
Assumption (e) is introduced in order to eliminate energy balances and other caloric calculations from the process in this step. In simulations with the full model, it is observed that the temperature is roughly constant throughout the process. This can be explained by the lack of significant heat input/output in the process, without units such as evaporators or condensers, in the process. The impact of the absorption enthalpy on the energy balances is partly compensated by pressure changes. For instance, the pressure drop between the absorber and the flash drum leads to a release of energy that is partly used to desorb the gas in the flash drum. This approximation is thus reasonable.
Reduced Model of the Absorption Column. A qualitative ternary diagram of the VLE for a typical CH 4 1CO 2 1solvent mixture at temperature T and pressure P ABS is shown in The calculation of the compositions and flow rates of all streams entering and leaving the absorption column requires the solution of two independent, and rather small, systems of equations: one VLE flash calculation to find the composition and flow rate of Stream (4), and a system of linear mass balance equations. Therefore, the calculations are highly robust and fast. Neither caloric calculations nor energy balances are required.
Reduced Model of the Flash Drum. Because of Assumption (c), the liquid phase in the VLE of the flash drum consists of pure solvent. As the vapor phase is ideal (Assumption (d)), the mole fraction x vap solv of the solvent in the vapor phase is given by where P FD denotes the pressure in the flash drum and P s solv ðTÞ the vapor pressure of the solvent, which is calculated using the SAFT-c Mie EoS. Once x vap solv is known, the flow rates and compositions of the liquid and vapor streams leaving the flash drum can be determined via mass balances.
All the steps required to characterize the process streams using the reduced model are listed in the Appendix; stream properties calculated with both the reduced and the full model are also compared. The differences are small enough to make the reduced model suitable for the generation of good initial guesses for the full model optimization.
Specifications for the Reduced Model. The specifications of the reduced model are given in Table 2. As the reduced model has fewer degrees of freedom than the full model, the number of specifications is decreased. The temperature T, assumed to be constant throughout the process, is an optimization variable. As the feed streams enter the process at about 300 K (cf. Table 1) and the adiabatic process does not contain any units for cooling, 300 K is used as a lower bound on the temperature. Although the pump provides a modest energy input to the process, no large temperature increase is expected and the upper bound is set to 350 K. The impact of these somewhat arbitrary bounds is discussed along the results of the case study.

Solvent design space
The set of all solvent candidates for the optimization is referred to as the solvent design space. In previous work, 27 the solvent space was the family of the binary mixtures of n-alkanes. This was in part due to the adoption of a modelling framework based on the representation of homologous series of compounds 68,69 within the SAFT-VR EoS, rather than a GC method, so that molecules from a different chemical family could not easily be modeled at the time. In addition, the direct solution of the CAMPD problem made it challenging to consider a broader design space. Although it was found that the CO 2 separation process could be operated with n-alkanes with an economically viable positive NPV, n-alkanes are rarely used 70 as agents for acid gas separation. One of the widely used commercial solvents for physical absorption of CO 2 is Selexol, which consists of poly(oxyethylene)dimethylethers. 67 Selexol has a better CO 2 /CH 4 selectivity compared to the n-alkanes and is thus likely to perform better in the process under current consideration. Driven by this motivation, the solvent space is extended to linear alkyl ethers. Molecules containing up to 22 alkyl groups and 21 ether groups are considered to ensure that all liquid linear alkanes and ethers are included. This is a broad set of structures within the families considered: it contains more than 200 candidate molecules with vapor pressures and melting points that make them suitable for use as a solvent in the process.
Although the solvent design space is limited to ether structures, it contains some interesting chemicals whose performance and suitability can be directly compared in the study: the family of n-alkanes (CH 3

Group definition
The final modeling component required before solving the CAMPD problem is a representation of the solvent design space. The functional groups used in this work are given with their SAFT-c Mie parameters in Table 3. They are chosen to represent CH 4 , CO 2 , and the solvent design space. A simple representation of the solvent design space could be obtained using CH 3 , CH 2 , and O groups. However, such a simple representation does not allow for the distinction of structural isomers, that is, structures with the same number of oxygen atoms but at different positions along the alkyl chain. Such isomers exhibit significant deviations in their physical properties. For instance, the vapor pressures of methylpropylether and diethylether, two structural isomers, can differ by up to 20%. 74 In order to allow a more accurate representation of the solvent properties, two O groups are defined: eO (end oxygen) represents an oxygen atom located next to a terminal CH 3 group in the molecule, and cO (center oxygen) represents an oxygen atom located between two CH 2 groups. An example of group partitioning is shown for dioxymethylenedimethylether in Figure 5.
The SAFT-c Mie parameters characterizing the like and unlike interactions of CH 2 and CH 3 were developed and presented in Refs. 62 and 73; they are adopted here without further modification. We also adopt the parameters of like and unlike interactions of CO 2 and CH 4 which were developed in other work and will be published elsewhere. 72 The parameters of the like interactions of cO and eO and the remaining parameters of the unlike interactions shown in Table 4 are developed here. We discuss the development of the cO and eO parameters due to the unusual treatment of the O group, but do not report in detail the development of the parameters for CO 2 and CH 4 .
To reduce the dimensionality of the parameter estimation problem, the eO and cO groups are assigned the same values of the parameters m Ã k , r kk , k a kk ; k r kk , and S k . The groups therefore only differ in the dispersion energy (potential depth) parameters eOeO and cOcO , respectively. The like interaction parameters are reported in Table 3. For groups k and l, the unlike interaction parameters r kl , k a kl , and most of the k r kl are calculated using appropriate combining rules (cf. Eq. 4). The unlike interaction parameters kl and k r kl that are adjusted to experimental fluid-phase data to improve the overall description are given in Tables 4 and 5.

Parameter estimation
The parameters of the new groups cO and eO are estimated from corresponding experimental data obtained from For stream labels, shown as superscripts, refer to Figure 3. the DETHERM database. 75 Details of the implementation, variance models, and error function are given in the Appendix. Estimation from Pure Component Data. In a first step, the vapor pressures and liquid densities (in either saturated or compressed states) of pure ethers are considered to obtain the like parameters for the eO and cO groups and their unlike interactions parameters with each other and with the CH 2 and CH 3 groups. The compounds used in the estimation and the temperature ranges of the experimental data are specified in Table 6. The original data sources can be found in Refs. 47, 74, 76-89.
The following seven parameters are estimated simultaneously from the data: S eO ; r eO;eO ; k r eO;eO ; eO;eO ; cO;cO ; eO;CH 3 , and eO;CH 2 . The remaining eO and cO group parameters are either fixed to constant values or coupled during the estimation according to the following constraints  (14) where the superscript CR denotes the value obtained by applying the relevant combining rule, cf. Eq. 4. Equations 12 and 13 are equivalent to imposing the same deviation from combining-rule behavior (which is often represented by the binary interaction parameter k ij for a geometric-mean rule) for the interaction energy of the eO and cO groups with the alkyl groups. The percentage absolute average deviation (%AAD) from the experimental data obtained with the optimal set of parameters is given in Table 6. The calculated vapor pressures and liquid densities are compared to experimental values in Figures 6-11. The calculated vapor pressures are found to be in good agreement with the data for the majority of the ether structures despite the small number of parameters, cf. Figures  6-8. Although the %AAD is larger than 10% in some cases, the model clearly provides good qualitative agreement with the data, which is crucial for solvent design. Structural isomers, such as methylpropylether and diethylether, are ranked correctly, cf. Figure 6. The deviations for the liquid density are generally below 1%, cf. Parameter Estimation from Gas Solubility Data. In a second step, high-pressure vapor-liquid equilibria and solubility data of CO 2 and CH 4 in the ether solvents are considered: the unlike interaction parameter eO;CO 2 is estimated from experimental solubility data of CO 2 in poly(oxyethylene)dimethylethers, 91,92 and the unlike interaction parameter eO;CH 4 is fitted to high-pressure VLE data of CH 4 and methylal (CH 3 OCH 2 OCH 3 ). 93 The following parameters were coupled during the estimation to provide values for the unlike interaction parameters cO;CO2 and cO;CH4 cO;CO2 CR Comparisons of the calculations with the SAFT-c Mie EoS and the experimental data 91-93 used in the estimation are given in Figures 12-14. A good description of the solubility of CO 2 is obtained with the SAFT-c Mie approach for the available datasets, cf. Figures 12 and 13. It can be seen in Figure 14 that the calculated VLE for the CH 4 1methylal mixture exhibits some deviations from the experimental values in the higher-pressure critical region of the mixtures. The pressure of the process in this study, however, never exceeds 10 MPa. The unlike repulsive exponent parameter k r eO;CH 4 was calculated to comply with the combining rule given in Eq. 4 during the estimation. The alternative of adding this parameter to the set of estimated parameters does not improve the quality of the description of the VLE significantly. Despite the general deviations in the VLE, the solubility of CH 4 in methylal is represented well (left part of Figure 14). The set of experimental solubility data of the mixtures of CO 2 and CH 4 with alkyl ethers is limited; additional measurements, particularly with longer ethers, would allow for a more detailed assessment of the proposed model.
Because of the deviations observed for the vapor pressure and the paucity of validation data for solubility, the proposed model is not suitable for the full detailed design of process equipment. Nevertheless, it is accurate enough to reflect the effects of the solvent chemical structure on the process   performance and trends in processing conditions. Due to the strong physical basis of the SAFT-c Mie approach and the comparably small number of parameters, the predictive capabilities of the thermodynamic model are expected to be wellsuited to identify novel solvent candidates. Once the optimal solvent candidate has been identified, the thermodynamic model should be adjusted if a more detailed process design is desired.

Representation of the solvent design space
Having identified the functional groups that are used to characterize the solvent molecules, we can now define the vector y of integer variables that represent how many times each functional group occurs in the solvent molecule. We first note that the number of CH 3 end groups, n CH 3 is fixed and equal to 2 for all molecules in the solvent design space as only linear molecules are considered. The integer decision variables are n CH 2 ; n eO ; n cO ½ where n k denotes the number of groups of type k. The bounds on these variables are listed in Table 1. The solvent molecule comprises a minimum of one CH 2 group based on the physical argument that a molecule smaller than propane would be too volatile and would not be suitable as a solvent. A maximum of two eO groups is allowed in the molecules, since there are only two CH 3 groups. Given that only linear molecules are considered in this case study, only one molecule feasibility constraint (g m ðyÞ 0)   needs to be imposed. We ensure that each cO group is associated with two CH 2 groups as follows n cO 2n CH 2 11 0 (17)

Hierarchical Solvent and Process Design for CO 2 Absorption
Optimization with the reduced model Problem Formulation. As the first step of the HiOpt approach, problem (MOO) is formulated and solved. The equations of the reduced model described above serve as the reduced process unit model (h 1 p ðz; yÞ ¼ 0Þ. The vector z of continuous variables is given by where T is the temperature of the isothermal reduced process model in K, and P ABS is the absorber pressure in MPa. The vector of integer variables is the same as for the full model, cf. Eq. 16. The bounds on z and y are given in Table 2. The reduced model provides neither internal flow rates in the absorption column nor energetic variables such as the    power consumption of the pump. The NPV, which is related to the sales of purified gas minus the capital and operating costs for the separation, 27 cannot be calculated due to this lack of information and the NPV cannot, therefore, be used as an objective in the optimization. To address this issue, three objective functions that represent different contributions to the NPV are used at this stage of the hierarchical optimization. The gas sales are directly proportional to the molar flow rate of Stream (3) (the "product flow rate"), which is, therefore, chosen as the first objective, to be maximized max z;y All molar flow rates F ðiÞ are specified in kmol/s. To account for the sizes of the process units, and thus, to penalize the capital cost of the absorber, pump, and flash drum, the recycle mass flow rate of solvent in kg/s (the "solvent flow rate") is used as a second objective to be minimized min z;y (20) where the molecular mass of the solvent in kg/kmol is denoted by MW solv . A final contribution to the overall process cost is the solvent loss. It is equal to the mass flow rate of solvent that must be supplied to the process and is penalized by minimizing a third objective min z;y f 1 3 ðz; yÞ ¼ F ð8Þ Á MW solv (21) The price of the solvent is not accounted for in this work. It is assumed that the differences in the production costs of the solvents are negligible. It is difficult to estimate production costs for chemicals without large errors, especially for components which are not being produced at present. The results, however, give an indication of the molecular structure of the optimal solvent and can be seen as a recommendation for further studies focused on the development of new cost efficient production processes for the optimal solvents.
The constraints g 0 ðz; yÞ 0 are defined according to the description of the case study where l is the solvent's viscosity (in cP) and x ð3Þ CH 4 is the mole fraction of methane in Stream (3). The viscosity is estimated by the method of Sastri and Rao. 94 Further, there is a constraint on the normal melting point T melt of the solvent (in K) for which two cases are considered: In the base Case 1, the process temperature is requested to be at least 10 K above the melting temperature of the solvent. It was found that this constraint is active during the optimizations. Thus, the accuracy of the estimated melting points is crucial for the validity of the results. In our current work, the GC method of Marrero and Gani 95 is used. In the Appendix, we show that the melting point is underestimated for some compounds in the solvent space. An additional safety margin of 20 K is added in Case 2 to reduce the risk of obtaining a process in which the solvent solidifies. One could also adopt the chance-constrained approach of Maranas 96 to treat the effect of this uncertainty more formally.
Computation of the Pareto Set. The reduced process model is implemented in gPROMS. 97 The SAFT-c Mie property model is implemented in Fortran 62 and accessed via a gPROMS Foreign Object. gPROMS provides an outer approximation/equality relaxation/augmented penalty (OAERAP) solver for MINLP problems, 98,99 which is used to solve all optimization problems in our work. The sandwiching algorithm for the MOOs is implemented in MATLAB. 100 The convergence tolerance for sandwiching is chosen to be d S ¼ 1:03. The algorithm identifies 10 Pareto-optimal points to approximate the Pareto set.
Results and Discussion. For Case 1 with melting point Constraint (24), the values of the three objectives for the 10 Pareto-optimal points are depicted in Figure 15. The product flow rate and the solvent flow rate are shown on the horizontal and vertical axes, respectively. Different ranges of values of the solvent loss are represented by the different symbols. Solutions that lie near the lower right-hand corner of the graph are desirable from the point of view of objectives f 1 1 and f 1 2 . The solvents corresponding to the Pareto-optimal points are listed in Table 7 using an identifier (ID); the IDs are also shown in Figure 15. The values of the continuous optimization variables in z at the Pareto-optima are also listed in Table 7.
Among the 10 Pareto-optimal points, there are six different solvent structures. These principally correspond to poly (oxymethylene)dimethylethers (IDs 1-5), the chemical fam-ily with the largest oxygen content among all solvents in the solvent space. This indicates the positive effect of oxygen on the process performance. All Pareto-optimal structures have an oxymethylene group (CH 3 OA) at each end. The solutions with ID 1-5 suggest a trade-off between the product flow rate and the solvent flow rate. In addition to the solutions with structures ID 1-5, the optimization based on the reduced model yielded only one other Pareto-optimal point (ID 6). Additional calculations confirm that the solutions with ID 1-5 are also obtained if the solvent loss is omitted as an objective, that is, if a two-dimensional Pareto set is calculated. Thus, the objective of minimizing the solvent loss is not strictly in conflict with the other two objectives. Indeed, if one were to start with solvent ID 1 and focus on improving the product flow rate (by moving to the right in Figure 15), then the solvent loss would be reduced at the same time, that is, the third objective would improve. When the solvent loss is minimized without considering the other objectives, solvent ID 6 is found. Note that the solvent losses in the two solutions with IDs 5 and 6 are almost the same.
The small solvent loss arising when molecule ID 6 is used can be explained by its large molecular weight and low volatility. Larger molecules are not feasible due to the constraint on the melting point. If high product flow rates are required, then molecules with a high oxygen content (ID 5) are favorable. A high oxygen content increases the selectivity of the solvent, so that methane loss is reduced. If a low solvent flow rate is desired, then the smaller methylal molecule (ID 1) is preferred. However, its rather small size leads to a high vapor pressure and the largest solvent loss of all Paretooptimal points.
In addition to the chemical structure of the solvent, the calculations yield optimal process conditions. The process temperature T hits the lower bound, or is very close to it, for all Pareto points. As expected for a physical absorption process of this type, this indicates that low temperatures are advantageous for all of the objectives. The solvent capacity is higher at lower temperatures, leading to a lower solvent flow rate. Lower temperatures also reduce the amount of solvent that is evaporated and lost in the regeneration step. Furthermore, lower temperatures lead to larger product flow rates because the selectivity of the Pareto-optimal solvents is higher. The choice of the lower bound on temperature, therefore, has a crucial effect on the results of the optimization. The numbers (IDs) next to the symbols correspond to the numbering of solvents in Table 7. The value of 300 K was estimated based on the fact that the process does not contain any cooling units and on the assumption that the temperature should not fall far below the temperature of the feed streams. The simulations of the Paretooptimal solvents/conditions using the full model, discussed in detail in the next subsection, reveal that the temperature in the process does not fall below 285 K in most cases. A reassessment of the MOOs with 285 K as lower bound for T is found not to change the final results of the study. The upper bound of 350 K has no influence on the optimization results. The pressure in the absorber at the Pareto points depends on the weighting of the product flow rate and the solvent flow rate. If the solvent flow rate is weighted heavily, solvent ID 1 is obtained and the pressure reaches its upper bound (7.5 MPa). A high-pressure difference between the absorption and desorption leads to a high cyclic capacity, that is, the concentration difference between loaded and lean solvents is high. By contrast, if the product flow rate is weighted heavily, solvent ID 5 is obtained and the pressure reaches its lower bound (1.0 MPa). This indicates that the selectivity of the Pareto-optimal solvents is better at low pressures.
For the highest product flow rates (ID 5) and the lowest solvent loss (ID 6), the optimization yields rather large molecular structures, cf. Table 7. The melting point Constraint (24) is active for these structures. The results of Case study 2 with Constraint (25) are given in the Appendix. The results are essentially identical to that of Case study 1. Constraint (22) on the viscosity is far from being active for every Pareto-optimal solution.
From the investigation based on the reduced model, it is clear that trade-offs between the different objectives are closely related to the chemical structure of the solvent and the operating conditions. The restriction to non-inferior Pareto points excludes a large set of solvents. Although obtained with a reduced model, the Pareto set provides valuable insights into solvent/process interrelationships. Furthermore, the Pareto-optimal set provides initial guesses for the next stage of the hierarchical design.
Optimization with the full model are a subset of the set x of continuous process variables. The bounds on these variables are given in Table 1. The NPV of the process serves as single objective f ðx; yÞ (for details on the cost calculation see Ref. 27). The constraints gðx; yÞ 0 contain the Constraints (22)- (24), which have already been used in the optimization with the reduced model. Similarly, the constraints h m ðyÞ ¼ 0 and g m ðyÞ 0 are unchanged. The full model includes hydrodynamic information to determine the size of the absorber. 27 Hence, an additional constraint on the column cross-section is included (cf. the description of the case study) where A ABS is the cross-section in m 2 .
Calculation of Initial Guesses. Any Pareto-optimal solution given by T as specifications would thus yield a solution that violates Constraint (23) on the product purity. To create feasible initial guesses for the optimization, a process simulation with the full model is performed for each Pareto-optimal point using the following specifications: The results of these simulations are given in Table 8, ranked in order of decreasing NPV. The solvent flow rate F ð10Þ is found to be higher compared to the solutions at the Pareto-optimal points (cf. Table 7): this increase is necessary to reach the required product purity with an absorption column of finite height. As the energy balances are taken into account, the temperature in the flash drum T FD is different from the 300 K specified in the Pareto-runs with the reduced model.
CAMPD Optimizations. The simulations of the 10 initial guesses yield feasible optimization variables y and z to start the solution of the MINLP based on the full model. Using more than one initial guess increases the likelihood of finding the global optimum or at least a high-quality solution. The computational effort of running the 10 independent optimizations identified in this case study is tractable, despite the complexity of each MINLP. For cases where a larger number of Pareto points are found, only the top ranked Pareto-optimal solutions should be taken as initial guesses to reduce the numerical effort.
The optimizations are performed in gPROMS. Six out of ten optimization runs converged including the top five in the ranking, cf. Table 8, and all six optima are identical. Additional information on the solver settings, the numerical performance and on the procedure to generate convergent simulations as initial guesses is provided in the Appendix.

Results and Discussion
The optimal solution is given in Table 9. The NPV of the process obtained with this solvent is clearly economically attractive with a NPV of 1721 3 10 6 US$. The optimal solvent, penta(oxymethylene)dimthylether (OME 5 , Pareto-ID 5), belongs to the set of Pareto-optimal solvents. The NPV is only slightly higher than the NPV obtained for the topranked initial guess given in Table 8. Note that solvent structure and the process conditions of these initial guesses stem from the Pareto optimization with the reduced model. This indicates that the reduced model and the three objective functions provide a good representation of the key cost contributions.
The optimal solvent OME 5 has a rather high oxygen content of 48 mass %. The optimization was repeated with a solvent space restricted either to the n-alkanes or poly(oxyethylene)dimethylethers (Selexol) for comparison. For the n-alkanes, the optimization yields n-decane as the best solvent, with an oxygen content of 0 mass % and a NPV of 871 3 10 6 US$. (In previous work, 27 using the same process model but a different property model, n-dodecane with a NPV in the same order had been found.) For the poly(oxyethylene)dimethylethers, the optimization yields tri(oxyethylene)dimethylether (CH 3 O(CH 2-CH 2 O) 3 CH 3 ) as the optimal solvent, with an oxygen content of 36 mass % and a NPV of 1577 3 10 6 US$. The trend that solvents of higher oxygen content perform better in the process is confirmed by these findings. It can be explained by the polar nature of ether oxygen, which results in a higher selectivity of the solvent toward CO 2 .
In this work, the production costs of all solvents are assumed to be identical to the production costs of n-alkanes, which are adopted from Ref. 27. The compound OME 5 is presently not produced in a large-scale industrial process. It can be produced from methylal and trioxane. 101 A conceptual design for the large-scale production process for OME 5 has been recently presented. 47 It would be interesting to estimate production costs on this basis to assess how this affects the CAMPD solution. We note, however, that the solvent losses in the process are small, and as a consequence the solvent cost may be neglected compared to the gas sales and the investment cost.
The solvent is assumed to be a pure substance. It was argued in Ref. 27 that it might be advantageous to use a mixture of several substances as solvent. The melting point of a solvent blend would be lower than that of the pure components allowing for a lower temperature bound in the process. Every component contributes specific properties to the solvent mixture which might lead to a synergistic effect on the overall performance. In addition, a mixture of similar substances might have lower production costs than a pure substance. Inte-grated solvent mixture and process design is a challenging task which is part of our future research activities.
The calculation of the pure component properties with the SAFT-c Mie EoS is excellent for OME 5 , cf. Table 6. The quality of the SAFT-c Mie predictions for the mixture, in particular the solubility for (OME 5 1 CH 4 1 CO 2 ) is not known at this point due to the lack of experimental data. The solubility of CO 2 in OMEs should be measured experimentally to assess the potential of the OMEs as solvents for CO 2 absorption.
The best solvent for the feed conditions and the flow-sheet topology considered in this work is OME 5 . If, for instance, the CO 2 concentration in the feed were changed or a second flash stage were to be added to the process to reduce the methane losses, then the optimization would need to be repeated. It is, however, reasonable to assume that other feed conditions or process topologies would also favor solvents with a high oxygen content.
The CAMPD analysis provides useful insights into the key performance drivers and is helpful to identify trends and promising solvent candidates. Studies of this type provide guidance in directing further experiments, modeling, and simulation. A realization of the optimal process of the CAMPD analysis would, however, require further experimental and design investigations.

Conclusions
The integrated, computer-aided, design of molecules, and processes, the problem (CAMPD), can often be formulated as a MINLP, albeit a challenging one. In our work, a novel hierarchical two-step methodology, HiOpt, is proposed to solve the CAMPD problem. In the first step, reduced models are derived for the unit operations of the process. If the original objective function can no longer be calculated from the reduced model, several surrogate objectives are defined which represent different contributions to the original objective. The Pareto-optimal set of best compromises between these objectives is determined using MOO. In the second step, the original MINLP problem is solved. The Paretooptimal solutions are then used as initial guesses for the solution of the full MINLP. This helps to make the problem more tractable and to identify favorable local solutions.
For the first time in CAMPD, a SAFT-based GC approach, the SAFT-c Mie EoS, 41 is used to model the physical properties of the substances involved in the process. The SAFT-c Mie EoS can be used to reliably represent the thermodynamic properties of the fluid phases for broad families of molecules with a relatively small number of parameters. The method eliminates the need for a large number of independent GC models for every property.
The proposed methodology is applied to the CAMPD of a physical absorption process to separate CO 2 from methane in natural gas. The solvent space contains the set of all linear alkyl ethers, comprising long-chain highly oxygenated compounds including commercial solvents and n-alkanes. Poly(oxymethylene)dimethylethers (OMEs), and in particular penta(oxymethylene)dimethylether, are identified as promising solvents for CO 2 absorption. They exhibit the highest oxygen content of all linear alkyl ethers. Further CAMPD studies for CO 2 absorption should include more oxygenated compounds, such as ketones, carbonates, methanol, ethanol, and other alcohols or their mixtures. The interaction parameters for the structural groups required to model these families with the SAFT-c Mie EoS are currently being developed. A rationale for the reduced model of the absorber is provided in this Appendix. The generic absorption column considered is shown schematically on the left panel of Figure 16. Stream F denotes the inlet gas, Stream S is the inlet solvent, Stream R is the outlet (clean) gas, and Stream L is the outlet (loaded) solvent. Streams L i and R i denote the liquid and vapor, respectively, that leave stage i of the column. The column consists of N stages. For every stage i, R i and L i are in VLE. The component mass balances around the overall column and around each stage are where F j denotes the molar flow rate of stream j, x j its composition vector. It follows from Eq. A1 that Thus, the "difference stream" P, defined as the difference between the liquid and vapor stream at a given height in the column, is constant. Note that P is a virtual stream. Its flow rate F P or some elements of its composition vector x P may be negative.
The compositions of the external and internal streams for a column with five stages are shown schematically in the ternary diagram on the right panel of Figure 16, in which the tie lines of VLE are represented. Both the pressure and the temperature are assumed to be constant. The virtual stream P is located outside the triangle, and is given as the intersection of the lines (FL) and (RS). Starting from the clean gas R, the compositions of the internal streams can be constructed based on the VLE between the gas R i and liquid L i and Eq. A2: R i and L i lie on the same tie line and all lines through R i11 and L i for any i intersect at P.
For the purpose of calculations with the reduced model in our current work, the compositions of Streams R, F, and S are known. The larger the number of stages present in the column, the greater the achievable CO 2 concentration of the loaded solvent L. There is, however, a thermodynamic limit for the CO 2 concentration of L: it can never exceed the CO 2 concentration of point L 1 that lies on the tie line pointing to the feed F (cf. Figure 16). Once the concentration L 1 is reached in the liquid phase, line (FP) coincides with the tie line through L 1 . The compositions of both phases no longer change below the stage at which this occurs. Such a stage is referred to as the pinch point. When a column is assumed to be of infinite height (N ! 1), as is the case for our reduced model, there must be a pinch point in the column at which L 5 L 1 .

Comparison of Process Simulation Results Obtained with Reduced and Full Model
An analysis of the accuracy of the reduced models used in our work with respect to the full models is provided here. Using the reduced model of the absorption column in process simulations invokes the following steps: 1. Specify input parameters, including the composition x ð2Þ and flow rate F ð2Þ of the feed, the absorber pressure P ABS , the process temperature T, and the desired purity of product x ð3Þ CH 4 2. Solve a VLE flash at P ABS , T under the constraint that the compositions of both phases x liq ; x vap lie on a line passing through x ð2Þ in the composition space to determine the flow rates of the streams leaving the column, F ð3Þ and F ð4Þ , and the flow rate of fresh solvent F ð10Þ .
A comparison of the concentrations calculated with the reduced and the full model for the solvent penta(oxymethylene)dimethylether, one of the Pareto-optimal solvents uncovered in our work, is shown in Figure 17, overlaid with the phase envelope at the corresponding temperature and pressure. The pressure in the absorption column and the purity of the product gas stream are set to be identical in both models. The temperature T in the reduced model is set equal to the temperature T ð4Þ at the bottom of the absorption column in the full model. The two models are in good agreement in terms of the stream concentrations. The main differences are that in the reduced model the lean solvent is purer and the CO 2 mole fraction in the purge Stream (6) is higher than in the full model.

Implementation of the Parameter Estimation
To estimate parameters of the chemical groups for the use in the SAFT-c Mie EoS and pure component data, the parameter optimization problem is implemented in the software gPROMS. 97 The SAFT-c Mie calculations are outsourced into a foreign object. gPROMS uses a maximum likelihood method to estimate the model parameters. Constant relative variance models are used for the experimental data with errors of 3% for the vapor pressure and 1% in the liquid density.
To estimate parameters of the SAFT-c Mie groups for the gas solubility data, a stand-alone code 102 for the calculation of the fluid-phase equilibrium is used. Therein the method of least squares with equal weighting of all experimental points is used.

Estimation of Melting Points
The method used to estimate the melting points (Marrero and Gani 95 ) is based on a different partitioning of the groups compared to the partitioning used in the SAFT-c Mie approach. The groups CH 2 , CH 3 , CH 3 O, and CH 2 O are used to represent the solvent molecules. We, therefore, define a vector n 0 k ¼ ½n 0 CH2 ; n 0 CH3O ; n 0 CH2O ; n 0 CH3 of numbers of groups that provides an equivalent representation, for use in the calculation of the melting point. The two representations are linked by a set of linear equations of type h m ðyÞ ¼ 0: During the solution of the optimization problems considered in our work, we find that the constraint on the melting point is active at the optimal solution. Thus, the accuracy of the estimation method is crucial. A comparison of the predictions and of the experimental values for the melting points of several compounds in the solvent design space is shown in Table 10. In some cases, the melting point is underestimated by up to 8.7 K. In other cases, however, an overestimation of the melting point by up to 10 or 15% is observed. Depending on whether a conservative solution or a more relaxed solution is desired, one may thus tighten or relax the melting point constraint. The tightening of the constraint is considered in our work as discussed in the main text and in the following section.

Optimization with a Tightened Melting Point Constraint
The Pareto-optimal points for Case 2, with the tightened Constraint (25) on the melting point, are calculated using the reduced model. The results in terms of the multiple objective functions are shown in Figure 18, and the solvent structures and process conditions are given in Table 11. The Pareto set exhibits  The value of T MG is estimated using the method of Marrero and Gani, 95 and T exp denotes the corresponding experimental values. 90 Figure 18. Pareto-optimal points of the MOO with the reduced model in Case 2 (tightened Constraint (25) for the melting point).
The numbers (IDs) next to the symbols correspond to the numbering of solvents in Table 11. similar characteristics to that obtained with the relaxed constraint of Case 1. With a high weight put on the solvent loss, a new solvent ID 7 replaces solvent ID 6 which violates Constraint (25) at the given process conditions. For a process temperature of 300 K, the solvent ID 5 also violates the Constraint (25). It is replaced in most cases by solvent ID 4 in the Paretooptimal set. For one combination of weights, solvent ID 5 remains the Pareto-optimal solution. In this case, the process temperature increases so that Constraint (25) becomes active. The positive effect of a larger molecular structure with a slightly greater oxygen content overcomes the negative effect of a slightly higher process temperature.

Implementation and Optimization Performance of the Full Process Model
The initialization strategy described here is implemented to converge the full process model within gPROMS based on the set of Pareto-optimal solvents generated and corresponding process information. The process model is initially set up at a refer-ence point with the specifications reported in Table 1 and with n-decane as a reference solvent n CH 2 ; n eO ; n cO ½ 0 ¼ 8; 0; 0 ½ To obtain solutions for other solvents/conditions, a dummy dynamic simulation is set up in which the specifications (A6) are changed continuously within some time horizon T H to the desired values, for example n eO ðtÞ ¼ n 0 eO 1t=T H n des eO 2n 0 where t is pseudo-time. We take advantage of the fact that the SAFT-c Mie implementation allows the specification of a noninteger number of occurrences of structural groups in a molecule. Although this has no direct physical meaning, it helps one to overcome numerical difficulties as the specification of a new solvent structure, which may lead to very different fluid-phase behavior from that found with n-decane as a solvent, can be achieved in a continuous, and thus, robust fashion. At the end of the dynamic simulation, the values of all of the process variables are saved for the initialization of the optimizations. The standard solution parameters of the OAERAP implementation within gPROMS are used. In all successful runs, the solver typically requires 2-3 MINLP iterations to identify a local solution. Optimizations with initial guesses with a solvent chemical structure and/or process conditions far from the optimum do not converge. In the course of the optimization in these cases, variables may take combinations of values that are infeasible inputs to the process model. Physically meaningful bounds on the variables aid in convergence of the simulations.