MIT Open Access Articles Optimization of a dual mixed refrigerant process using a nonsmooth approach

This article uses a nonsmooth ﬂ owsheeting methodology to create simulation and optimization models for dual mixed refrigerant processes. New improved operating conditions are obtained using the primal-dual interior-point optimizer IPOPT, with sensitivity information calculated using new developments in nonsmooth analysis to obtain generalized derivative information using a nonsmooth generalization of the vector forward mode of automatic differentiation. Several optimization studies are performed with constraints on both the minimum temperature difference ( D T min ) and total heat exchanger conductance ( UA max ) used to represent the trade-offs between energy consumption and the required heat transfer area. In addition, comparison is made with the conventional process simulator Aspen HYSYS using particle swarm optimization. Results show that the nonsmooth model was able to reduce the required compression power by 14.4% compared to the initial feasible design for the dual mixed refrigerant process, and by 20.4 e 21.6% for the dual mixed refrigerant process with NGL extraction. Furthermore, the solutions obtained from the nonsmooth model were 1.9 e 8.1% better than the design obtained by particle swarm optimization. Multistart optimization also shows that IPOPTconverges to the best known solution when starting from an initial feasible design.


Introduction
Natural gas is projected to play a central role in the shift towards green energy sources. According to the BP Energy Outlook from 2017 [1], it is the fastest growing fuel at 1.6% p.a, and is expected to overtake coal in overall demand by 2035. Considered a cleaner alternative to oil and coal due to a low sulfur content, lower CO2 emissions and no particle emissions, the main challenge with natural gas as a fuel is related to transportation. For shorter distances, pipelines are usually preferred, though it requires investments in infrastructure and limits the potential customers to predetermined locations. However, in order for natural gas to be traded as a global commodity, long distance transport using liquified natural gas (LNG) is necessary. LNG has also received increased attention as the solution to production of natural gas in remote offshore fields where pipeline transport to an onshore process facility is uneconomical. Nevertheless, liquefaction of natural gas is a very energy intensive process. Investments in expensive and proprietary technology such as cryogenic multistream heat exchangers and turbomachinery are required, and together with high operating costs, the liquefaction process generally accounts for about 30e40% of the total cost in the LNG chain [2]. The type of liquefaction process used largely depends on the purpose of the LNG plant. In general, LNG production plants are categorized into three types: base-load, peak-shaving and small-scale plants. For small-scale production, the capital costs are of primary concern, thus promoting designs with relatively few pieces of equipment such as single-mixed refrigerant (SMR) processes [3]. Compactness of the design and the small equipment inventory also make SMR processes attractive for floating LNG (FLNG) systems. In base-load production, on the other hand, the operating costs are significant and should be reduced by using more efficient designs. Typical processes considered for base-load production are therefore cascade, propane precooling (C3MR), and dual mixed refrigerant (DMR) processes.
Large-scale production of LNG demands more efficient liquefaction processes. In particular, as the train capacity goes up, operating costs will increase and account for a higher fraction of the total costs of liquefaction. Efficient design and operation are therefore needed to limit the total costs of these plants. Dual mixed refrigerant processes represent a promising technology for largescale production of LNG. With an additional multicomponent refrigerant mixture, they achieve small temperature driving forces in the precooling and liquefaction heat exchangers, as well as added versatility to the feed gas composition. Their compact design compared to other large-scale processes, i.e. refrigerant cascade and C3MR processes, also means that DMR is suitable for largescale FLNG units. Efficient and versatile designs are priorities in such large-scale LNG facilities at the expense of higher capital costs (though mitigated by economies of scale) and a larger equipment inventory compared to SMR processes. Although DMR processes already have been installed on large-scale FLNG processes such as Shell Prelude and Coral FLNG [4], they have received relatively little attention from the optimization community [5,6].
The high costs associated with production of LNG have resulted in numerous optimization studies for different liquefaction processes that are summarized in extensive reviews by Austbø et al. [5] and He et al. [6] Particular attention has been given to the Poly-Refrigerated Integrated Cycle Operations (PRICO) process [7], which is an SMR process consisting of one multistream heat exchanger (MHEX) and a simple vapor-compression refrigeration cycle. Its relatively basic design has made it a desirable test case for different tailor made MHEX models for optimization of LNG processes.
Kamath et al. [8] developed an equation-oriented (EO) MHEX model using concepts from pinch analysis such as composite curves. The model treats the MHEX as a heat integration problem with hot and cold utility duties set to zero. The authors then performed energy targeting by using the Duran and Grossmann model for simultaneous process optimization and heat integration [9]. Furthermore, the authors included a disjunctive model for phase regime detection that is represented in terms of complementarity constraints. Phase regimes of the inlet and outlet streams may vary during simulation and optimization of LNG processes, representing a modeling challenge as the phase regimes traversed in the MHEX are not known a priori. This also affects other process equipment such as compressors and valves that handle streams leaving the MHEX in a priori unknown phase regimes. Pattison and Baldea [10] proposed a different MHEX model using a pseudo-transient EO approach. Assuming a fixed relative sequence of stream temperatures prior to optimization, they constructed a series of enthalpy intervals for the composite curves. Temperatures for each enthalpy interval are then calculated from thermophysical property models by introducing a pseudo transient temperature variable and reorganizing the model equations into a system of differential-algebraic equations (DAEs). No disjunctive representation is needed for handling phase changes in the model. Instead, the model perturbs the time variable across the phase boundaries while keeping the temperature constant. Properties are then resolved with the solution from the previous time step as an initial guess. Hasan et al. [11,12] proposed a MHEX model using a superstructure of twostream heat exchangers and solving a heat integration problem with no external utilites. However, the model can only handle phase changes as long as the phases traversed in the MHEX are known a priori. Another superstructure model by Rao and Karimi [13] addresses this issue by including nonlinear constraints to ensure phase boundaries occur at the endpoint of each heat exchanger such that it operates in a specific phase regime. A process simulator (e.g. Aspen Plus) is used for the property correlations and therefore no disjunctive constraints are needed for phase detection.
The tailor-made MHEX models have so far been tested for SMR processes with a particular focus on the PRICO process. In addition, most of them require solving nonconvex nonlinear programs (NLPs) [8,10,13] or mixed-integer nonlinear programs (MINLPs) [11,12]. The complexity of DMR processes, which includes additional refrigerant streams, MHEXs, and splitting or mixing of streams, make them significantly more difficult to optimize. A large temperature span and phase changes in the MHEXs also result in a highly nonlinear variation in the enthalpy as a function of temperature. To account for this nonlinearity, the process streams are frequently partitioned into smaller stream segments with constant heat capacity flowrate. Examples presented by Watson et al. [14] show that using less than 20 segments to represent the two-phase region for the PRICO process leads to significant inaccuracy in the simulations resulting in an under-prediction of the necessary compression power. Therefore, additional variables must be added for each MHEX making the optimization problem large and challenging to solve. As a result, most optimization studies on DMR processes use conventional process modeling tools such as Aspen HYSYS [15] or Aspen Plus [16] for simulation and optimization, or together with an external software for optimization [5].
Khan et al. [17] optimized the specific compression power and total UAusing multi-objective optimization for a DMR process. Aspen HYSYS was used for process modeling, whereas Matlab was used for finding the Pareto-optimal solutions. Results show that multi-objective optimization decreases the specific compressor power and total heat exchanger conductance (UA) by 13% and 3% compared to the base case. Another optimization study was done by Hwang et al. [18] on a DMR process using Aspen HYSYS for process simulation and a hybrid optimization method combining a genetic algorithm with sequential quadratic programming (SQP). Process optimization resulted in a 34% reduction in total compression power compared to the original design. Morin et al. [19] optimized an SMR and a DMR process using an evolutionary search method and sequential quadratic programming. SQP was shown to obtain better results on average, particularly for the SMR process. However, it struggled to converge successfully for the DMR process without significant tuning, and was found to be higly sensitive to the problem formulation. The evolutionary search method, on the other hand, required relatively little tuning and was faster than SQP for the DMR process. Furthermore, it obtained solutions within 3.12% of the SQP solution. The authors therefore suggest using a hybrid optimization strategy, where an evolutionary search method would be used for obtaining a good starting point for the SQP algorithm. Additional optimization studies of DMR processes are reviewed in Austbø et al. [5] and He et al. [6].
This article uses the MHEX model developed by Watson et al. [20] that employs recent advances in nonsmooth analysis for handling composite curves and heat transfer area calculations. Similar to the models by Kamath et al. [8] and Pattison and Baldea [10], the heat exchanger is treated as adiabatic with no external utilities. Based on this, the authors reformulated the Duran and Grossmann [9] model into a single nonsmooth equation to prevent temperature crossovers in the MHEX. Equations for the overall energy balance and area calculations are also included in the model, reducing the degrees of freedom by two compared to the MHEX modules in Aspen HYSYS and Aspen Plus. Phase regime detection for the streams is handled using the nonsmooth mathematical operators min, max and mid [21], where the mid operator is a function that maps to its median argument. The MHEX model employs a hybrid simulation framework where two-phase substreams are solved sequentially in flash calculations nested in an EO environment. Flash calculations are solved using a nonsmooth extention of the inside-out flash algorithm by Boston and Britt [22] that handles phase-changes without relying on post-processing methods for phase detection. This procedure is summarized in a 3-paper series by Watson et al. [14,23] and Watson and Barton [24], where a mid-function is used for automatic outlet phase regime detection within the flash calculation procedure. The algorithm was shown to handle instances close to the critical point for which the conventional inside-out algorithm implemented in Aspen Plus with post-processing methods for phase detection failed to identify the correct single phase [23]. The purpose of the model is to provide an alternative to preliminary design and analysis of LNG processes currently conducted by commercial process simulation tools, e.g. Aspen HYSYS. No detailed geometrical considerations are therefore studied here. Instead, the model should be employed as an initial step in the design process to obtain good operating points for the design, as well as a rough cost estimate through installed compression power and total heat exchanger conduction UA, and then use more detailed models in the detailed engineering stage.
Model size remains an important issue in the existing multistream heat exchanger models. Binary variables or complementarity constraints for phase detection, and a MINLP formulation for enforcing the second law have resulted in these models only being used for studying the relatively simple PRICO process. For larger and more complex processes, the state-of-the-art is still to use commercial simulation tools, e.g. Aspen HYSYS, sometimes with an external optimizer. The nonsmooth modeling framework, however, exhibits several advantageous properties. As the flash calculations are nested in the flowsheet calculation, the number of segments chosen to model the two-phase region does not impact the model size. Furthermore, the reformulated Duran and Grossmann model places a single constraint on the problem, regardless of the number of stream segments used in the model. Consequently, additional stream segments can be included for an accurate representation of the two-phase region without increasing the number of variables in the model. The result is a multistream heat exchanger model that scales well with the number of streams and stream segments used in the analysis. Furthermore, auxiliary equipment, i.e. valves, compressors and mixers employ the same nonsmooth inside-out algorithm, and thus do not impact the overall size of the flowsheet model. The model was shown to work well for simulating [25] and optimizing [26] complex SMR processes. In this article, these models are expanded to that of a DMR process and improved using IPOPT [27]. Optimization studies of these processes have primarily used Aspen HYSYS together with an external derivative free solver. Others have used SQP, although it suffers from a low success rate, particularly without any prior tuning of the model [19]. Optimization cases with constraints on the total heat exchanger conductance UA max and minimum approach temperature DT min are considered in the analysis, and the performance of the algorithm is compared with the state-of-the-art, using optimization results from Aspen HYSYS with particle swarm optimization (PSO) [28,29].

The dual mixed refrigerant process
The process studied in this article is a version of the AP-DMR process with single-stage compression for both refrigerant cycles. The process flowsheet is presented in Fig. 1 along with the variables considered in the optimization. It features three MHEXs for cooling the natural gas, though in practice the two cold MHEXs are normally merged into a single spiral-wound heat exchanger (SWHX). The SWHX is a cryogenic heat exchanger, commonly used in Dual Mixed Refrigerant (DMR) processes, that consists of one or several stream bundles wound around a mandrel. Separate tubes are used for the hot refrigerant and natural gas streams, whereas the cold refrigerant is flowing countercurrently outside the bundles. Normally, different refrigerants are used for cooling different parts of the heat exchanger in order to reduce the refrigerant circulation at the cold end. As a consequence, the heat transfer area required to cool the same quantity of natural gas is comparatively smaller. A warm mixed refrigerant (WMR) cycle consisting of ethane, propane and n-butane is used for precooling the natural gas and the cold mixed refrigerant (CMR) to an intermediate temperature T OUT HP;1 . After precooling, a CMR that consists mainly of nitrogen, methane, ethane and propane, is used for cooling the natural gas down to 120.15 K. For a better distribution of refrigerant in the SWHX, the CMR is sent to an adiabatic separator, where the liquid product is subcooled, throttled and used for providing cooling in MHEX 2. The vapor product, on the other hand, proceeds to MHEX 3 where it is condensed, subcooled and throttled to provide cooling in the cold end of the process.
The following nomenclature is used for the parameters and variables of the refrigerant streams in the DMR model: High pressure level of the warm (W) and cold (C) mixed refrigerants: P HP;W=C . Low pressure level of the warm (W) and cold (C) mixed refrigerants: P LP;W=C . Outlet temperatures of the high pressure refrigerant and natural gas streams from MHEXs 1 and 2: T OUT HP;1=2 . Outlet temperatures of the low pressure refrigerant from MHEXs 1,2 and 3: T OUT LP;1=2 . Component molar flowrates: f W=C . Total molar flowrate: F W=C .

The nonsmooth multistream heat exchanger model
The nonsmooth MHEX model is a generalization of the twostream countercurrent heat exchanger model [20], defined by the equations for the energy balance, minimum temperature difference DT min and total heat exchanger conductance UA. The energy balance for multiple hot and cold streams is given by Equation (1) where mCp H=C are the heat capacity flowrates for the hot (H) and cold (C) streams. Area calculations can also be defined analogously for MHEXs by assuming vertical heat exchange between the hot and cold composite curves [20].
Here, UAis the heat transfer conductance, DT LM is the log-mean temperature difference, K is the total number of enthalpy intervals and DQ k is the enthalpy change over interval k.
Location of the minimum temperature difference is ostensibly more challenging to calculate for MHEXs. For single phase heat transfer in a two-stream countercurrent heat exchanger, the minimum approach temperature occurs at the end points of the heat exchanger. Multiple stream inlets entails possible pinch candidates at all inlet temperatures. Furthermore, MHEXs are normally included in cryogenic applications with phase transitions, where the minimum approach temperature can also occur at interior points in the heat exchanger. Instead, the minimum approach temperature can be calculated using pinch analysis and heat integration with a pinch location algorithm. Different pinch location algorithms have been developed. However, most entail solving either a nonconvex NLP or an MINLP to global optimality. A reformulation of the simultaneous optimization and heat integration algorithm by Duran and Grossmann [9] was proposed by Watson et al. that translates the optimization problem into a single nonsmooth equation [20]: where P is the (finite) set of candidate pinch points and EBP p H=C are the enthalpies of extended hot/cold composite curves for pinch candidate p as defined in Watson et al. [20] as nonsmooth functions. The equation can be solved using nonsmooth numerical equation solvers, and no optimization problem must be solved to ensure feasible heat exchange.
Phase changes in MHEXs represent a known modeling challenge, as phase boundaries vary dynamically during simulation and optimization. Consequently, the phases traversed in the heat exchanger are not known a priori and must be determined on every iteration of the solver. Instead of using binary variables or disjunctions for detecting the correct phase regime, the MHEX model partitions each process stream into superheated (sup), two-phase (2p) and subcooled (sub) substreams whose inlet and outlet temperatures are determined by the nonsmooth equations [21]: where T IN=OUT are the inlet and outlet temperatures of the process stream, T IN=OUT sub=2p=sup are the corresponding inlet and outlet temperatures of the substreams, and T DP and T BP are the dew and bubble point temperatures of the process streams. Additional stream segments are used to improve the accuracy of the calculations. In the simulation and optimization studies in this article, five stream segments are used for the two single-phase regimes, whereas 20 stream segments are used for the two-phase region.
Stream temperatures in the two-phase region are calculated from pressure-enthalpy (PQ-) flash operations. Similarly, auxiliary equipment such as compressors and valves are modeled using isentropic and adiabatic flash calculations, respectively. Because of the aforementioned issue with phase regime detection and dynamic phase boundaries in MHEXs, the PQ-flash algorithm must be capable of handling instances of single phase flow. A nonsmooth extension of the well recognized Boston and Britt algorithm [22] that handles cases of single-phase flow is presented in a 3-paper series by Watson et al. [14,23] and Watson and Barton [24]. The algorithm includes a mid-function for correctly identifying the phase regime at the outlet, avoiding the use of the suggested post-processing methods for handling single-phase behaviour. In particular, the flash algorithm was shown to handle cases near the critical point robustly, where the conventional inside-out algorithm was prone to failure [23]. In addition, a methodology for calculating generalized derivatives analytically for the flash calculations, means the flash calculations can be nested and solved as subroutines in the model. The result is a smaller model size, where auxiliary equipment such as compressors and valves, as well as stream segments for the two-phase region are solved sequentially at every iteration.

Simulation of the DMR process
Simulation Case I from Vikse et al. [30] was used as a starting point for the optimizer. The MHEX model [20] provides either two or three equations, depending on whether area calculations are included. Using a two-equation MHEX model is particularly useful for precooling and intermediate MHEXs during simulation as specifying either the UA-value or DT min in these heat exchangers is difficult. The process parameters in these MHEXs depend on the process conditions in the low temperature MHEX, which makes it hard for the user to select a pre-defined area or minimum approach temperature. Instead, during simulation of the DMR model, an efficient approach is to use the two-equation MHEX model to find the approach temperature for a set of initial conditions, and calculate the UA-value subsequently. Then, during process optimization, a three-equation MHEX model can be used to iterate around the initial area values as desired. No tear equations are required in this model as the pressure levels, material flows and compositions are set for the refrigerant streams similar to an equation oriented approach. In addition, the temperature is fixed after the cooler. The results of the simulations are discussed for each of the two variable sets below. Consequently, the DMR simulation model has seven equations, which were used to solve for the following unknown variables: T OUT HP;1 , T OUT HP;2 , T OUT LP;3 , P LP;C , F C , f W; ethane and DT min;2 [30].
The model converged to a feasible solution with 17.34 MW in total compression power after three iterations and a total simulation time of 68.7 s [30]. The results of the base case is presented in Table 1 and will be used as a future reference for the optimization studies. The variables solved for in the model are presented in italic. Furthermore, the total work, the individual compressor duties (W W=C ) and the UA 1=2 -values are calculated during postprocessing. The compressor duties are calculated using isentropic flash calculations, as well as a specified isentropic efficiency (h). The UA-values are calculated for MHEX 1 and 2 using Equation (2). The warm mixed refrigerant composition changes in the simulation as a result of varying the component molar flowrate of ethane (f W;ethane ). The corresponding natural gas composition, flowrate and pressure level are specified in Table 2.

Formulation of the optimization problem
Appropriate formulation of the optimization problem is important when analyzing LNG liquefaction processes. Design of heat exchanger networks frequently use a lower bound of the minimum approach temperature DT min as a measure of the level of heat integration in the process. A smaller DT min means more heat integration at the expense of a larger total heat transfer area and thus higher capital costs. Therefore, by placing constraints on the minimum approach temperature in the process one can manipulate the trade-offs between energy consumption and the required heat transfer area. For this reason, DT min has frequently been used as a specification in optimization of LNG processes. However, thermodynamic irreversibilities increase with both increasing driving forces and decreasing operating temperature [31,32], and as a result, the distribution of driving forces in the MHEXs is important. Jensen and Skogestad [33] showed that using a minimum temperature approach when optimizing LNG processes led to suboptimal utilization of the heat exchanger conductance. Later, Austbø and Gundersen showed that a problem formulation constraining the heat exchanger conductance rather than DT min , resulted in a better driving force distribution in the heat exchanger [32]. Kim and Gundersen [34] later expanded the results to optimization of dual mixed refrigerant processes, though here the total UA-value should be used rather than the heat exchanger conductance values for the individual MHEXs.
Another important consideration in optimization of LNG processes is the degree of superheating for the low pressure refrigerant. A certain degree of superheating is required to avoid the formation of liquid droplets in the compressor and excessive compressor wear. However, a large safety margin may degrade the optimal solution. Kim and Gundersen [34] plotted the specific compression power as a function of the minimum superheating value, and found that the solution was more sensitive to the minimum degree of superheating when using a DT min approach.
Furthermore, they discovered that the optimal superheating values for the WMR and CMR when using area constraints for the DMR process were 10 K and 12 K, respectively. A minimum superheating of 5 K is used for all optimization cases presented in this paper. The overall optimization formulation is similar to that of SMR processes [26]: DT sup ! DT sup;min ; x LB x x UB ; where DT sup is the degree of superheating for the two compressors and his a set of equations describing the DMR model. This set includes Equations (1)e(3), as well as energy balances for the individual stream segments. As the two-phase stream variables, and auxiliary equipment such as valves, mixers and compressors, are solved using nested flash calculations, their functions are not included in h, but are instead resolved at every iteration. The functions in hare described using the nonsmooth operators min, max and mid for modeling phase changes, and the minimum approach temperature in the MHEXs. The UA functions are given by Equation (2). Sensitivity based methods using classical derivatives or finite difference approximation are therefore prone to failure when encountering points of nondifferentiability. A traditional approach has been to either approximate the function around the kink points by using smooth approximations, or by reformulating the nonsmooth operators using disjunctions. Alternatively, generalized derivative elements can be used, which are extensions of the concept of derivatives to certain classes of nondifferentiable functions. One such generalized derivative is the Clarke Jacobian [35] that is defined for functions that are locally Lipschitz continuous on their domain. Another is the lexicographic (L-)derivative, which is defined for functions that satisfy conditions for lexicographic (L)smoothness [36]. The model presented in this article uses an automatic differentiation framework for calculating lexicographic directional (LD-)derivatives developed by Khan and Barton [37]. The lexicographic directional (LD)-derivative is a generalization of the classic directional derivative that is computed sequentially along the directions provided by a directions matrix. More importantly, the L-derivative can readily be obtained from the LDderivative. For an extensive review on evaluating LD-derivatives and their applications, the reader is directed to Barton et al. [38]. Nonsmooth optimization algorithms exist in the literature [39]. In particular, solvers can be divided into two main categories; subgradient methods and bundle solvers. Subgradient-based methods work similar to smooth methods (e.g. the steepest descent method) but where the gradient is replaced by an arbitrary subgradient of the function. Subgradient methods are advantageous in that they exhibit low storage requirements and are relatively easy to implement. However, convergence of the algorithms can be slow, there exists no rigorous stopping condition, and the selection of step size is challenging, predominantly due to the possibility of having nondecending step directions [40]. Rather than using an arbitrary subgradient at each point, bundle methods approximate the subdifferential of the function. This is achieved by gathering the subgradients obtained at previous iterations into a bundle. Different bundle methods have been proposed [40]; the proximal bundle solver [41], bundle Newton method [42] and a limited memory bundle method [43,44]. Implementations of the bundle Newton method for linearly constrained optimization and the limited memory bundle method for bound constraints exist. The proximal bundle method [41] is suitable for nonsmooth constrained optimization problems. Tests conducted on different convex and nonconvex nonsmooth optimization problems favored the proximal bundle method both in regards to efficiency and reliability [40]. However, attempts at using the MPBNGC v2.0 [41] solver for the nonsmooth single mixed refrigerant processes were unsuccessful [26]. Surprisingly, the interior-point optimizer, IPOPT [27], proved to be more suitable for optimizing the nonsmooth flowsheet models, despite its assumption that the objective function and constraints are twice continuously differentiable. Nevertheless, there are issues regarding the use of IPOPT for optimizing the nonsmooth models, particularly with regards to the termination criteria [26]. Specifically, no implementation currently exists of the nonsmooth analog for the dual feasibility calculations in IPOPT, such that termination at points of nondifferentiability is problematic [26]. However, Watson et al. [26] showed that the termination issues are avoided when the dual feasibility tolerance is suitably high. The authors also provided a list of changes to the default settings, which empirically was shown to provide better performance for the nonsmooth models. The full list of non-default settings is provided in Table 3.

Optimization of the base case with a UA max formulation
The first example looks at improving the operation of the feasible design from Table 1 by varying the pressure levels, refrigerant compositions, and intermediate temperatures for the same maximum total heat exchanger conductance (UA max ). The models are programmed in the Julia v0.6.0 programming language and run on a Dell Latitude E5470 laptop in the Ubuntu v16.10 environment with an Intel Core i7-6820HQ CPU at 2.7 GHz and 8.2 GB RAM. New and improved operating conditions are obtained using both the nonsmooth approach and PSO together with the process simulator Aspen HYSYS. However, as pointed out by Vikse et al. [25], differences exist between Aspen HYSYS and Aspen Plus, where the latter uses the same property package as the nonsmooth model. Whereas simulations in Aspen Plus showed close correlation to results from the nonsmooth models, results in Aspen HYSYS deviated significantly at lower temperatures [25]. In particular, this was accredited to Aspen HYSYS using different ideal gas enthalpy correlations from Aspen Plus and the nonsmooth models. However, Aspen HYSYS also includes the possibility of loading property models directly from Aspen Plus. Therefore, simulations in Aspen HYSYS are done with underlying property models from Aspen Plus to better compare the results from PSO with the nonsmooth models. Termination criteria for the PSO algorithm were set to 500 maximum iterations, a function tolerance of 10 À8 and a maximum number of stall iterations of 40. The latter is the number of iterations with relative change in objective function value less than the overall function tolerance before termination.
Accurate representation of the process streams is ensured by adding an appropriate number of stream segments in the MHEXs. Five segments are used for representing the single-phase liquid and vapor substreams and 20 stream segments for the two-phase region in each MHEX. However, as the two-phase stream variables are solved sequentially in nested flash calculations, they do not impact the total model size. As a result, the nonsmooth optimization model consists of 108 variables, of which P 3 i¼1 s i ðn sup þ n sub À 2Þ ¼ 88(where s i are the number of process streams in MHEX i and n sup=sub are the number of segments for the superheated and subcooled regions) are the temperatures of the individual stream segments and the remaining 20 variables are other process variables such as compositions, pressure levels and intermediate stream temperatures. Variable bounds for the process variables are presented in Table 4. Upper and lower bounds on the temperature of the individual stream segments are set to 350 K and 100 K, respectively.
A solution to the DMR process was obtained with IPOPT after 135 iterations and a total run time (including the time needed for initializing and simulating the DMR model to obtain the initial feasible design) of 1173 s. The new operating conditions resulted in Table 3 Non-default settings for IPOPT used in this work [26].   Fig. 2. The composite curves in the new design are significantly closer in MHEX 3, though slightly less so in the intermediate MHEX.
The result is a design where the overall temperature difference is closer to being proportional to the operating temperature, thus reducing the exergetic losses due to irreversible heat transfer as discussed by Austbø and Gundersen [32]. Lower refrigerant flowrates and a higher degree of separation before the intermediate MHEX also result in less self refrigeration and a smaller overall heat exchanger duty. The DMR model was also optimized using PSO with the same model built in Aspen HYSYS. The same initial feasible design and variable bounds were used to compare the performance of the two optimization strategies. A total run time of 20.9 h was needed before the algorithm obtained a solution with a total compression work of 15.11 MW, 1.8% higher than the solution with the nonsmooth model (NM). PSO obtains a process with larger approach temperatures in all three MHEXs. The high pressure level of the cold mixed refrigerant is also significantly higher than what was obtained with the nonsmooth model. Morin et al. reported a similar difference in objective values between using the evolutionary search method and SQP for optimization of DMR processes [19]. Improved results from the two optimization methods are presented in Tables 5 and 6 along with the data from the initial feasible design. The UA max calculated from the simulation, is used as an upper bound on the total heat exchanger conductance. Consequently, any improvements must come from a redistribution of the available heat transfer area in the individual MHEXs.
Next, the DMR model is optimized for different maximum total heat exchanger conductance values, using both the nonsmooth approach and PSO for comparison. Variable bounds and IPOPT settings remain the same as in Tables 3 and 4, and the analysis is done with maximum total heat exchanger conductance values of 6, 9 and 12 MW/K, respectively. Solutions were obtained for the nonsmooth model in the three cases. However, as the initial point is the base case in Table 1 with a UA max ¼ 5:39MW/K, the solution will vary considerably from the initial design for larger heat exchanger conductance values. Consequently, the number of iterations and hence the total solution time is larger. The first case converged after 106 iterations and a total run time of 838 s. The second case required 221 iterations to solve, resulting in a total run time (including the time needed for initialization and simulation to obtain the initial feasible design) of 2067 s. For the last case, with a UA max ¼ 12MW/K, a solution was obtained in 137 iterations, corresponding to a total run time (including initialization) of 1158 s.  Results for the three cases are presented in Tables 7 and 8. The driving force plots for the three cases are shown in Fig. 3. As can be observed, the driving force profiles show a similar shape for the three cases, as well as for the base case. However, with increasing heat exchanger conductance values, the obtained designs become tighter with smaller driving force temperatures throughout the process.
The three cases were also optimized using Aspen HYSYS with PSO. For UA max ¼ 6MW/K, a solution was obtained after 22.3 h with an objective function value of 14.84 MW, which is 2.4% higher than what was obtained by the nonsmooth model. In the second case, PSO obtained a design after 25.3 h with a total compression power of 13.79 MW, or 1.7% higher than the corresponding nonsmooth solution. PSO obtained a solution with 13.26 MW for the third case, which is 1.1% higher than the solution obtained by the nonsmooth model. The optimization algorithm required 23.9 h to solve this case. It should be mentioned that PSO terminated after reaching the maximum number of iterations in all cases above. Consequently, no local optima were obtained, and the PSO algorithm terminated with the best known solution within 500 iterations. A function tolerance of 10 À6 was also tested, however, in that case the PSO algorithm converged to significantly suboptimal points compared to the nonsmooth models.

Optimization of the base case using a DT min formulation
In this section, the UA max constraint in Equation (7) is replaced by a DT min specification for all three MHEXs. Again, the feasible design in Table 1 is used as a starting point for the optimizer, with the bounds as specified in Table 4. Furthermore, the minimum approach temperature was specified to 3.50 K in all three MHEXs. IPOPT obtained a solution for the DMR model with a DT min formulation after 188 iterations and a total run time of 2363 s.
The solution requires a total compressor work of 14.51 MW and a total heat exchanger conductance of 11.54 MW/K. By comparison, a UA max formulation with UA max of 11.54 MW/K, resulted in a total compression power of 13.16 MW or 9.3% less than the solution predicted by the DT min formulation. Thus, energy is saved without extra investment, just by using a problem formulation inspired by thermodynamics that utilizes the heat exchanger area more efficiently. The composite curves and driving force plot for the new design are presented in Fig. 4. The new improved solutions for the nonsmooth model with the DT min and UA max formulations are summarized in Table 9.

Optimization of the dual mixed refrigerant process with NGL extraction
The last case considers another version of the DMR process where NGL extraction is included between MHEXs 1 and 2 (see Fig. 5). The process has been simulated and optimized for the feed gas compositions in Table 10. Better operating conditions are obtained using the formulation in Equation (7) with a UA max constraint. As the natural gas stream data is fixed in the optimization, the number of optimization variables remains unchanged from the previous model. However, an additional constraint is added to ensure a lean LNG product with more than 89% methane.
Flowsheet simulation is used to obtain feasible starting points for the optimizer. The following unknown variables were solved for in the simulation model: T OUT HP;1 , T OUT HP;2 , T OUT LP;3 , P LP;C , P LP;W , P HP;C , and UA 3 . A warm mixed refrigerant composition with 47.83% ethane, 34.17% propane and 18% n-butane was used for all three feed compositions. The molar flowrate was set to 1.55 kmol/s. The cold mixed refrigerant flowrate and composition also remained constant for the three cases with 7% nitrogen, 41.80% methane, 33.20% Table 7 Solutions of the nonsmooth model for different UAmaxvalues.    Tables 11 and 12. Optimization was done using the simulation results as a starting point and the calculated UA max as an upper bound on the total heat exchanger conductance. As the natural gas stream data are assumed fixed in this study, the decision variables in the model remain the same as in the base case. However, most of the variable bounds are changed from the base case, and are provided in Table 13 for clarification.
The first case considers a relatively light natural gas composition with 87.60% methane content. IPOPT obtained a solution after 300 iterations and a total run time of 2684 s including the time needed for initialization and simulation. The solution requires a total compression power of 15.57 MW, which is a reduction of 20.9%   compared to the initial feasible design. With an NGL extraction temperature of 220.00 K, the design also satisfies the requirement of a methane content larger than 89%. In the second case, a richer composition with 85.6% methane content is used. The optimization model converged after 600 s and 62 iterations to a solution with a total power requirement of 14.90 MW, which is a 21.1% reduction compared to the initial design. In the third case, a rich natural gas composition with 83.6% methane content is assumed where significant NGL separation is required to obtain a satisfactory LNG specification. A solution was obtained after 1794 s and 173 iterations. The results for the three cases are summarized in Table 14. Driving force plots of the three cases are also presented in Fig. 6.

Convergence characteristics
Optimization of DMR processes is significantly more challenging than for SMR processes. Derivative based methods in particular tend to struggle to converge, mainly due to nonlinear constraints and the size of the models. Morin et al. [19] reported this using an SQP algorithm for optimizing a DMR process in Aspen HYSYS. Tuning the algorithm for the particular model required attention from the designer and was therefore quite time consuming. However, once properly set up it achieved about a 3% decrease in objective function value compared to using a stochastic search method. The model in this paper is a standalone simulation and optimization tool tailor made for LNG processes. Consequently, it does not involve interfacing commercial software with a solver. Instead, the model uses a nonsmooth flowsheeting strategy and LD-derivatives to compute sensitivities at points of nondifferentiability. Improvements are done using the interior point algorithm IPOPT with sensitivity calculations provided by the LD-derivatives. However, IPOPT assumes twice continuously differentiable objective functions and constraints, which represents an issue with the termination criterion. The dual feasibility calculations are not valid at nonsmooth points [26], which can cause the algorithm not to converge and instead iterate in a negligibly small search space. As mentioned previously, this issue is resolved by increasing the dual feasibility tolerance to 0.1, thus a local optimum cannot be verified for these nonsmooth optimization models with a gradient based optimizer. However, in order to avoid this issue altogether, a nonsmooth optimization solver that can handle L-derivatives must be developed. Global optimality cannot be guaranteed for the solutions as a local solver is used for optimizing the nonsmooth models. Therefore, in order to evaluate the quality of solutions, multistart using 15 runs from randomly generated (and possibly infeasible) starting points was done for the cases with variable UA max values. The results show that none of the converged solutions were better than the one provided using the feasible starting point. For the case with a UA max value of 6 MW/K, multistart resulted in four runs converging to the same solution, two cases converging to a suboptimal design, and the remaining nine runs not converging at all. This corresponds to a success rate of 27%, which is slightly better than the 10e20% range reported by Austbø and Gundersen [32] with UA max constraints for the PRICO process. For the cases with heat exchanger conductances of 8 MW/K and 15 MW/K, six of the runs converged to the best known solution, whereas two runs converged to suboptimal solutions. This corresponds to a success rate of 40%. A slightly higher success rate was obtained with the DT min formulation, where 45% of the runs converged to the best known solution. Therefore, a possible strategy would be to relax the UA max constraint first, and use the solution as an initial guess before optimizing the driving force distribution in the MHEXs.
Setting bounds for the optimization variables is challenging for the DMR process due to the constraints and number of variables in the model. In particular, finding a good set of bounds on component molar flowrates, refrigerant pressure levels and minimum approach temperatures is challenging, and may influence the convergence. Simulation provides good initial points to the optimizer. In addition, it makes it easier for the user to set bounds that can increase the convergence ratio. In the case studies performed here, the pressure, DT min and temperature bounds are relaxed to include the solutions of different case studies. Significant improvements to the convergence can be made, however, through stricter bounds tailored for a specific case. The latter becomes easier by first considering the initial feasible design, and then determining the bounds accordingly. As no other tuning is necessary for IPOPT besides the settings provided in Table 3, finding bounds and optimizing the models is still less time-consuming and yield better results than running Aspen HYSYS with stochastic search methods such as PSO. The optimization studies here are only considering local optimization and no global solution can be guaranteed. However, analysis using multistart, as well as comparisons to PSO in Aspen HYSYS, show that the optimizer obtain good results despite lacking an analogous dual feasibility relation for nonsmooth functions.
In the optimization studies, some of the decision variables ended up on their respective lower bounds. As the bounds were relaxed further, the model experienced failures in the flash calculations, as the model approached regions for which no solution was obtained using the inside-out flash algorithm. Consequently, it imposes a limitation on the current implementation. Developing an alternative model formulation or heuristics for handling these extreme flash conditions to prevent failures in the overall model, is therefore a necessary next step.

Conclusions
Several stand-alone models for LNG processes have been reported in the literature, and these have been tested for single mixed refrigerant processes. However, these models often involve solving either a mixed integer nonlinear program or a nonlinear program with a large number of variables and constraints. Dual mixed refrigerant processes feature designs with small temperature differences, several multistream heat exchangers, intermediate throttling, and separation. The models are therefore significantly larger and more complex than single mixed refrigerant processes, and as a result, optimization studies in the literature tend to use commercial process simulators such as Aspen HYSYS along with a derivative based or stochastic search method. The nonsmooth flowsheet models presented here are stand-alone tools for both simulation and optimization of LNG processes, where functioning models already have been developed for different single mixed refrigerant processes. In this article, the nonsmooth flowsheeting strategy is used for developing a simulation and optimization model for dual mixed refrigerant processes. Different cases were studied using both DT min and UA max constraints, and improvements were done with IPOPT using sensitivity information provided by LD-derivatives. The new operating conditions resulted in a 14.4% decrease in total compression power for the DMR process and 20.9e21.6% for the DMR process with NGL extraction. The nonsmooth model also obtained 1.1e2.4% more energy efficient operating points than PSO in 25e80 times less solution time. Multistart for several of the cases showed that IPOPT obtained the best (known) solution when using a feasible starting point. Furthermore, the model achieved a higher success rate with a UA max constraint than what was reported by Austbø and Gundersen for the basic PRICO process [32].