Modular Framework for Simulation-Based Multi-objective Optimization of a Cryogenic Air Separation Unit.

A framework to obtain optimal operating conditions is proposed for a cryogenic air separation unit case study. The optimization problem is formulated considering three objective functions, 11 decision variables, and two constraint setups. Different optimization algorithms simultaneously evaluate the conflicting objective functions: the annualized cash flow, the efficiency at the compression stage, and capital expenditures. The framework follows a modular approach, in which the process simulator PRO/II and a Python environment are combined. The results permit us to assess the applicability of the tested algorithms and to determine optimal operational windows based on the resultant 3-D Pareto fronts.

Introduction

The commodities derived from air separation units (ASUs) are essential for many supply chains worldwide. For instance, He et al.[1] highlight the relevance of these units by reporting that they represent ∼5% the total power consumed in China. ASUs are an integral part of many processes and can be integrated into various systems such as natural gas systems[2−4] and storage power plants.[5,6] The commodities obtained from an ASU have important applications in the manufacture and health care sectors. Indeed, their production level has noticeably increased over the years, and it is projected to keep growing.[7,8]

Air is mainly composed of argon (Ar), nitrogen (N2), and oxygen (O2). Relevant uses of Ar can be found in welding applications such as tungsten inert gas (TIG) and gas tungsten arc welding (GTAW) because it generates an inert media, minimizing defects in final joints. Moreover, Ar is employed in fluorescent light bulbs,[9] the production of silicon semiconductors (as a protective gas),[10] and as a carrier in inductively coupled plasma spectroscopy.[11] On the other hand, the uses of N2 are mainly in the production of fertilizers, nitric acid, and nylon.[12] It is utilized as an inert gas to suppress the risk of explosions within a combustible material.[13] During the melting process of aluminum, N2 works as a substitute for chlorine and Freon to prevent oxide formation.[14] Finally, the main uses of O2 entail support systems in metallurgy, polymer synthesis, and medicine. Applications in metallurgy include laser-oxygen cutting of mild steel.[15] In addition, oxygen plasma treatment is applied to modify the surface of different polymers to improve their adhesive properties.[16] Regarding medical applications, hyperbaric oxygen therapy contributes in the healing of open wounds,[17] promoting an antimicrobial effect for suppressing bacterial growth on tissues.[18] More recently, the demand for medical O2 has increased during the COVID-19 outbreak in 2020–2021 because of its use in intermediate and intensive care units for mechanically ventilated patients.[19]

An ASU separates air into its principal components (Ar, N2, and O2) through a low-temperature liquefaction process. These products might vary in their composition, phase, and produced quantities. The selection of suitable operating conditions is crucial to guarantee the desired final product characteristics. Several contributions have performed simulation-based analyses of ASUs. For instance, Zhu et al.[20] implemented a low-order dynamic model for a cryogenic distillation column, primarily applied in N2 purification, using HYSYS as process simulator. Teague and Edgar[21] suggested a predictive model of a small pressure-swing absorption (PSA) system of air separation to address technical needs and operation in military aircraft applications. Jiang et al.[22] optimized a PSA separation system while tailoring optimization algorithms. Huang et al.[23] proposed a nonlinear model predictive control (NMPC) to adjust the operating conditions of a PSA system to overcome variations in the product demand.

Seeking to enhance the system’s design, Van Der Ham and Kjelstrup[24] improved the heat integration of the distillation columns by moving the low-pressure column down along the high-pressure column. Such changes decreased ∼23% in the total entropy production. Manenti et al.[25] studied the intensification of an ASU using the process simulator PRO/II. Main accomplishments included improving the O2 purity, recycle of the rich Ar stream, and a feasibility analysis of energy generation. Aneke and Wang[26] studied the incorporation of heat recovery cycles with different configurations at the compression stage to enhance the system’s energy efficiency. Fu et al.[27] conducted a process analysis and optimization based on a complete equation-oriented approach. Tesch et al.[28] evaluated the integration of LNG regasification into air separation processes to enhance the compression stage. An exergy and economic analysis was carried out using Aspen Plus. Finally, Young et al.[29] performed a detailed design and economic evaluation of a cryogenic ASU, identifying the relevance of each stage when computing the capital expenditures. Although open source state-of-the-art platforms, such as FOQUS,[30,31] allow the interaction with commonly used chemical engineering process modeling software (Aspen Plus, gPROMS, Thermoflex) for simulation-based optimization, the functionalities of such platforms are still under development. In this sense, a generic integration using multiobjective optimization (MOO) evolutionary algorithms and a high resolution ASU model in PRO/II has not being fully explored in the literature.

MOO algorithms permit to explore highly nonlinear systems with complex solution domains.[32−34] A comparative study between different evolutionary algorithms demonstrated that the nondominated sorting genetic algorithm III (NSGA-III) and the unified nondominated sorting genetic algorithm III (UNSGA-III) can achieve competitive results when compared with other approaches.[35] The performance of these algorithms was tested for an open-shop scheduling problem with resource constraints.[36] On the other hand, the multiobjective evolutionary algorithm based on decomposition (MOEAD) demonstrated superior performance in benchmark problems, such as the traveling salesman,[37] and in passive vehicle suspension optimization.[38] Similarly, Gaspar-Cunha et al.[39] applied a MOO methodology to minimize the cycle time, length of weld and warpage in an injection molding process using NSGA-III, UNSGA-III and MOEAD. A number of contributions have relied in these evolutionary optimization algorithms,[40,41] and the Python package pymoo(42) incorporates them for easy deployment.

In this work, three objective functions are optimized including the annualized cash flow (CF), the efficiency of the Rankine cycle (ef) at the compression stage, and the capital expenditures (CAPEX) of the facility. The first two objectives are maximized while the third is minimized. The achieved results provide guidance regarding the most adequate setup of operational variables. In a previous study, the reference vector-guided evolutionary algorithm (RVEA) was implemented to address the MOO of an ASU.[43] In this context, major contributions of the present work are summarized next:

A unified framework integrating the process simulator PRO/II and a Python environment for the MOO of an ASU, offering a variety of optimization algorithms, is proposed. It includes results handling and visualization.

The reliability of the high-resolution simulation is assured by verifying and replacing the main simulation at each iteration step.

A comparison under fair conditions between evolutionary algorithms for simulation-based MOO using quality indicators is detailed.

Operational windows for flexible ASU operation based on Pareto front results are provided.

The remainder of this contribution is organized as follows. The ASU process and its simulation strategy are described in section . The MOO problem formulation, its objective functions, MOO algorithms, quality indicators, and the proposed framework are presented in section . Thereafter, we present the results and findings in section . Finally, conclusions and future work are outlined in section .

Air Separation Unit and Compression System

Air separation can be performed with different technologies. Among them, the PSA, vacuum PSA, and cryogenic ASUs are the most relevant. In this contribution, a cryogenic ASU is considered. The proposed system incorporates air purification for the inlet stream, compression, and air separation. Cryogenic separation requires rigorous heat integration and an appropriate design of the separation columns to achieve an efficient setup. Figure illustrates the process flow diagram of the cryogenic ASU under study.

Figure 1

Process flow diagram of the air separation unit including the compression and purification system.

Process Description

Air impurities need to be removed from the feed stream before entering the compression section. Common impurities include carbon dioxide, water, and low molecular weight hydrocarbons (methane, ethane). The presence of impurities in the downstream stages can lead to maintenance problems and even equipment explosion.[44] A common method for air purification is molecular sieve adsorption, which mainly requires synthetic zeolites because they have the needed pore size for impurities retention. This front-end purification system consists of two vertical vessels (T-301 A/B) combined with a set of switch valves. The configuration of these valves depends on production requirements and maintenance scheduling. The affinity of a molecular sieve decreases over time, requiring regeneration. For this step, a stream of heated N2 gas is injected in counter-flow, at approximately 523 K. Depending on the available N2, regeneration could operate around a temperature of 407 K and pressure of 300 kPa. The regeneration gas is diffused uniformly over the surface of the bed for preventing accumulation of the residual components. Once the process concludes, a direct contact aftercooler (DCA) tower reduces the temperature of each molecular sieve.[45]

Thereafter, the pressure increases to 600 kPa before entering to the cooling section. This is achieved using three compressors (C-101/102/103) with an average compression ratio of 1.95:1 each. The heat generated is removed by inter coolers. Thus, the compressors operate close to isothermal conditions. The heat exchangers (E-101/102/103) reduce the air temperature to 313.15 K in the first two and to 303.15 K downstream of the heat exchanger E-103.

Heat is converted into electricity using a Rankine cycle (RC). The reduction of power consumption, obtained from the RC, is usually 10%. The efficiency of the compressors and the adiabatic expansion is assumed 80%.[26] As the temperature of the discharged air raises, the amount of energy generated increases as well.[46] Nevertheless, this configuration requires more power to attain the same pressure, limiting the reduction of energy consumption up to 0.2%. At the end of the compression stage, the air stream splits, and 10% of the flow rate goes to the turbine expander.

The separation stage consists of two distillation towers: high pressure (HP) and low pressure (LP) columns. The HP column (T-201) operates at 600 kPa. The flow rate of the bottoms product contains approximately 35% of O2, 1% Ar and processes about 60% of the main feed. It has a total condenser allowing a liquid reflux through the distillation column. A reboiler is not considered in the design as the air stream exiting the bottoms of the column exchanges heats with the upper stream before entering the Ar column.

The LP column (T-202) operates at 150 kPa to increase the relative volatility of O2 and N2. The composition of O2 in the stream is 99% as the side draw removes the Ar. The liquid N2 in T-201 supplies the reflux of T-202.[45] The Ar column feed comes from the bottoms of the LP column. The draw rate is about 20% of the air feed rate, and 4% is removed as Ar product. The pressure drop available for this distillation is limited by the temperature of the expanded rich liquid from the HP column.[44]

Process Simulation Strategy

The selection of the most suitable equation-of-state (EOS) allows us to avoid inaccuracies in the thermodynamic calculations of the system. Skogestad[47] suggested employing the Soave–Redlich–Kwong (SRK) EOS for processes with nonpolar components such as O2 and N2. On the other hand, Aneke and Wang[26] utilized the Peng–Robinson EOS available in Aspen Plus to simulate an ASU. At low-density gas regions, both equations fit experimental pressure–temperature–volume data accurately and offer a similar reliability on the whole range of temperature.[48] For the system under study, the SRK EOS handles the multicomponent vapor–liquid equilibrium. The simulation of the system starts at the end of the pretreatment section. The composition of the inlet air stream at the compression section is detailed in Table . The process parameters and variables considered in the simulation are detailed in the process description and are listed in Table .

Table 1

Feed Flow Rate Composition and Fraction (mol/mol)

component	molar fraction
N2	78.09
O2	20.95
Ar	0.93

Table 2

Process Parameters and Variables for the ASU Simulation in PRO/II

parameter	value	units
δT approach of subcoolers	1	K
δT approach of condenser reboilers	1	K
trays in HP column	44
trays in LP column	69
trays in argon column	55
pressure drops in distillation columns	10–15	kPa

The system is simulated using the commercial process simulator PRO/II in steady state. The high-resolution model incorporates quality constraints to ensure that products comply with the desired quality standards such as the degree of purity. In addition, two calculators are incorporated to simulate and include reference streams into the multistream heat exchanger E-202. In this sense, the path of each flow stream is traceable within the process, facilitating further analyses. Regarding the column trays efficiency, Biddulph[49] suggests an efficiency of 70% for all trays in ternary mixtures such as air. However, Zhu et al.[50] assumed a tray efficiency of 100% for an ASU with crude Ar separation. In this work, an efficiency of 100% for all trays is considered. Finally, pseudostreams are included on each distillation column to simulate energy integration in the whole heat exchanger network.

Simulation-Based MOO

Problem Formulation

The optimization problem considers three conflicting objectives; the annualized cash flow, efficiency of the Rankine cycle, and the capital expenditures. The first two objective functions are maximized, f1(x) and f2(x), while the third, f3(x), is minimized.

The optimization problem is written as followssubject towhere the equality constraints are denoted by h(x) and are essentially incorporated by the process simulation (mass and energy balances) and techno-economic evaluation. The inequality constraints, g(x), include the convergence verification of the simulation, design, and feasibility constraints and other constraints incorporated to the objective functions. The vector x ∈ R11 represents the decision variables vector, and it is constrained between a lower and upper bound xLB and xUB, respectively. The objective functions and the decision variables are described in more detail next.

The first objective function f1(x) is the cash flow, CF(x), which is calculated by subtracting the operational expenditures OPEX(x) of the process from the generated revenues:The revenues are calculated by summing the product of the ith produced air constituent’s price c and its flow rate F(x), with i = 1, ..., P. In this case, we consider three air constituents, O2, N2, and Ar, i.e., P = 3. The OPEX(x) is defined in eq .

The ef at the compression stage, ef(x), is the second objective function f2(x). The ef(x) is the ratio between the work generated by the expander Wexp(x) and the sum of the heat removed by the jth heat exchanger Qhex, (x) with j = 1, 2, ..., q, in this case q = 3, at the compression

Finally, the capital expenditures, CAPEX(x), correspond to the third objective function f3(x). It dennotes the capital investment required for the facility. The CAPEX(x) is calculated as follows; the sum of the eth cost bare module CBM, (x) of each rotatory or static equipment e = 1, ..., E. This total summation is multiplied by 1.18 times the ratio between the chemical engineering plant cost index (CEPCI) of 2019 and a base year (2001):

The operational expenditures OPEX(x), mentioned in eq , consist of 18% of the CAPEX(x), plus 1.23 times the raw material costs (RMC) and the utility costs (UC(x)), which is

Table enlists the prices of the air constituent products and the cost of the pretreatment for air conditioning. Table shows the utilities cost. The UC(x) and the parameters for computing the CAPEX(x) of each equipment were obtained from Turton et al.[51]

Table 3

Prices of Products and Raw Materials

process stream	price	units
O2	155.06	$/Ton
N2	86.51	$/Ton
Ar	195.00	$/Ton
air pretreatment	10.00	$/Ton

Table 4

Utilities from CAPCOST

utility	price	units
low-temperature refrigeration	8.49	$/GJ
electricity	18.72	$/GJ
refrigeration moderated to low	4.77	$/GJ

The vector x ∈ R11 gathers the 11 decision variables considered in this work, which also correspond to the degrees of freedom (DOF) of the problem (11 DOF). Table depict them with their respective lower xLB and upper bounds xUB, and the physical units, when it applies. Here, TE-101/102/103 corresponds to the outlet hot product temperature of the heat exchangers E-101, E-102, and E-103, respectively. LxE1/2 is the liquid fraction of the streams entering at stages 1 and 28, respectively, in the column T-202. ΔTdew is the hot product temperature increase above the dew point of the stream entering at stage 33 in the column T-202. Fcool is the flow rate of the coolant at the compression stage. R1/2 is the ratio of the coolant flows that split before entering the heat exchangers E-101 and E-102, respectively. Since R3 is defined by the equation R3 = 1 – R1 – R2, it is not considered as a decision variable. PC-104 is the discharge pressure of compressor C-104, and PP-101 is the discharge pressure of pump P-101.

Table 5

Upper and Lower Bounds of the Decision Variables of the ASU System

decision variables	xLB	xUB	units
TE-101	288	303	K
TE-102	288	303	K
TE-103	288	303	K
LxE1	0.8	1.0	-
LxE2	0.41	0.71	-
ΔTdew	3.0	5.0	K
Fcool	4.0	6.0	kmol/s
R1	0.1	0.4	-
R2	0.1	0.4	-
PC-104	150	250	kPa
PP-101	900	1100	kPa

Evolutionary MOO Algorithms

MOO addresses the balanced evaluation between conflicting objectives. In practice, it is not possible to attain a single optimal point capable of simultaneously optimizing all objective functions under study. In this sense, a set of optimal points offers a broader and flexible result. These optimal points could be selected on the basis of specific needs. In contrast, a common approach to solve multiple objective problems is to build a single function in the fashion of a weighted sum. However, the magnitude of these weights depends on an expert opinion who could provide subjective importance to each objective, introducing bias.

As an alternative, a set of optimal solutions permits us to overcome such limitations and to have a broader picture of the problem solution. This set of solutions is referred to as the Pareto-optima points, corresponding to nondominated solutions. Such a characteristic implies that no objective can improve without penalizing at least another one. The Pareto-optima points delimit the Pareto front (PF), which is the boundary between feasible and unfeasible solutions.[34]

The use of evolutionary multiobjective algorithms incorporates parallelization features such as vectorized evaluation, threaded-loop evaluation, and distributed evaluation. Such features contribute to debottlenecking the optimization routine. When considering a Python-scientific computing environment, the package pymoo shows competitiveness when interconnected to simulation-based systems. The framework offers updated single- and MOO algorithms and different features related to such optimizations following an object oriented programming strategy. Additionally, visualization and decision making are included in the package.[42]

In this work, the resultant PFs are illustrated as 3D surfaces because of the three evaluated objective functions. The following optimization algorithms are taken into consideration:

Nondominated Sorting Genetic Algorithm III

The Nondominated Sorting Genetic Algorithm III (NSGA-III) employs a mating-restriction scheme and establishes a multiple targeted search in advance. This feature ensures the diversity of the solutions. The NSGA-III characterizes by the combination of the parent and offspring populations. Thus, it allows us to conserve the fittest members of the previous generation. More details of this algorithm are provided by Deb and Jain.[52] The selected parameters for testing this algorithm are a population size of 78 individuals and 24 generations. To verify its accuracy, three iterations are carried out.

Unified Nondominated Sorting Genetic Algorithm III

The unified approach of NSGA-III (UNSGA-III) holds more flexibility in terms of the number of objective functions in the optimization problem. The normalization operators and niching-based selection procedure automatically defunct for mono-objective problems and become active for multi and multiobjective problems. The selected parameters for testing this algorithm are the same as those in in NSGA-III.

Multiobjective Evolutionary Algorithm Based on Decomposition

Finally, the multiobjective evolutionary algorithm based on decomposition (MOEAD) keeps the diversity of the trade-off solutions by setting weight vectors randomly. Each population member in the resulting subproblem is correlated with a weight vector. Zhang and Li[53] proposed two metrics referred as Penalised-Based Intersection (PBI) and Tchebycheff (TCH). The PBI metric is the measurement of the distance between a point and the ideal point produced by the weighted sum of perpendicular distance d2. Results are obtained by adding the distance d1 along the reference direction as depicted in the following equation

The TCH metric considers an ideal point denoted as z* and the weight vector w. The objective function of the jth subproblem is described as

The selected parameters for this algorithm are 20 neighbors, 0.7 probability of neighbor mating, and 24 generations per iteration. To verify its accuracy three iterations are carried out. These settings are selected for matching the total simulations/evaluations performed by the other tested algorithms and to compare them under fair conditions.

Quality Indicators

Four quality indicators (QIs) are selected to illustrate the performance of each optimization algorithm. These indicators include convergence, spread, uniformity, and cardinality.[54]

Two main approaches to quantify convergence consider the PF availability or not. If available, the distance to the PF (or some designated points) to a solution set are quantified. When the PF is not available, the convergence is estimated through the Pareto dominance between solution sets. The alternative strategy applied in this work is based on the Pareto dominance relation provided by the solutions, which is known as the C-indicator[55]where A and B correspond to the compared solution sets, and a ≤ b means that a weakly Pareto dominates b. In summary, the C-indicator gives a ratio of the number of elements in B that are weakly dominated by an element in A. As this quality indicator does not follow a commutative behavior, both C(A, B) and C(B, A) are calculated.

Uniformity is a QI that considers the even spacing between elements of a solution set. The uniformity of a solution set describes the distribution and its ability to represent appropriately a PF. The most common uniformity indicator of a solution set A = {a1, a2, ..., a} is known as spacing (SP),[56] and it is defined as followswhere the distance d(a, A/a) is the minimum 1-norm of the element a ∈ R to the rest of the elements in A, m is the number of objective functions as shownand d̅ is the mean of all d(a, A/a), with i = 1, ..., n.

The pread of a solution refers to the cover of a solution set. The maximum spread (MS) provides a measure of the range obtained in a solution setwith a, a′ ∈ R, and m as the number of objective functions. Thus, the spread is directly proportional to the maximum extent of each objective function.[57]

Finally, cardinality refers the number of elements in a solution set. The cardinality of solution is defined by the number of nondominated solutions in the set (N). The error ratio (ER) indicator proposed by Van Veldhuizen[58] establishes the proportion aswhere a is an element of the solution set A and the function e(a) is defined as followsIn this case, the expression a ∈ PF is assumed to be the nondominated solutions.

Framework Architecture

The proposed framework uses a first-principles high-resolution commercial simulator to obtain data from a process simulation. Thereafter, within a Python ecosystem further analyses take place. Building simulations with this kind of software are a common practice for process engineers. However, with the advent of more powerful computers and the availability of more and cutting edge scientific computing tools, it is possible to obtain better insights. Currently, it is possible to run simulators recurrently with Python programs using all kinds of optimization, statistical, or data science functionalities. This allows us to build more flexible and sophisticated designs and architectures considering more complex and realistic scenarios, which is an extraordinary basis for decision-making.

The interaction between PRO/II and the Python ecosystem occurs by using the Component Object Model (COM) implementation. The COM is a technology developed by Microsoft for building an interface in which different software components can interact. Certain industry standards use this technology, like CAPE-OPEN for process simulation technologies. The PRO/II simulator provides a Simulator COM Server that allows full read/write access to the PRO/II simulation database. Thus, any language or application that supports COM can use these functionalities provided by the PRO/II COM Server. The Python for Windows extensions package, or pywin32, allows accessing Window’s COM to control and interact with other applications. Once the communication is established, the PRO/II COM Server grants access to read and write in the objects and streams contained in the process simulation. The ability of manipulating a simulation model from Python permits either concurrently or iteratively evaluate the high-resolution model. The implemented framework in this work is illustrated in Figure . Similar frameworks for solving other PRO/II process simulation problems have been proposed by Jones et al.[59] and the authors.[34,43,60]

Figure 2

Framework with the main elements for the evolutionary MOO of an ASU simulation using the coupling of Python and PRO/II simulator.

Because the PRO/II simulations are manipulated from Python, it is important to guarantee a well-developed process simulation as described in section 2.2. PRO/II provides a graphical user interface to comply with this task. Once the base simulation is set, from a Python program the decision variables (operating conditions) are set accordingly to the MOO algorithm.

The Python code consists of two main functional blocks, the evaluation of the objective functions and the main program. In the main program, the hyper-cube domain in R11 is set. Following the documentation of pymoo,[42] a MyProblem class is defined from the base class Problem. In the MyProblem class, the init reserved method of certain dimension (i.e., 11) defines upper and lower bounds of the decision variables (i.e., between 0 and 1), number of objective functions, number of constraints, and other settings. In addition, it is required to define a method written as evaluate, which implements the mechanism of evaluating the cost function to return a vector with the cost values. In these methods, we use the function ASU that implements the multiobjective function, which is detailed later. Thereafter, the program creates several objects, problem of the class MyProblem, an algorithm object algorithm defining the optimization method (i.e., NSGA-III, UNSGA-III or MOEAD). Finally, the function minimize is called with problem and algorithm as part of the arguments. This runs the solver of the minimization problem and returns the resultant PF. A code snippet of the main block is depicted in Chart .

Chart 1

The interaction between PRO/II and Python occurs in the objective function ASU. A snippet of the objective function ASU is provided in Chart . The function receives a vector of decision variables x and returns a vector of values of the objective functions y ∈R3.

Chart 2

The program connects through the PRO/II Server to the simulation database. Then the 11 values given in the vector x are set into the simulation database, and the changes are committed. At this point, the database is closed and the simulation runs. Once the simulation finishes, the database is opened again and the resulting values, required to calculate the objective functions, are retrieved from the PRO/II database. Finally, the output vector of the functions y = [f1, f2, f3] is calculated and returned to the optimization algorithm. This function is iteratively executed by the function minimize until the stop criteria is met.

Sometimes an input can produce an unfeasible simulation setup, corrupting the simulation file. To avoid this problem and to evaluate the system under fair conditions a replacement mechanism places the original simulation before each iteration. Moreover, this technique ensures the independence of each simulation during the optimization process and prevents the original file from being unusable. At the end of the process, the results are stored in data frames using pandas, to be used later for postprocessing and visualization.

Results

The computational tools to test the performance of the proposed framework are the process simulator PRO/II Process Engineering 10.2 (64 bit) and Spyder (Python 3.7.6). The experiments were performed in a laptop PC Intel Core i7–8565U CPU @ 1.80 GHz with 16.00 GB of installed RAM. The problem considers two different constraint setups. The first evaluations are unconstrained. The second set of results include a constraint that guarantees only positive results of the CF objective function. Each evaluation establishes a 3D-figure and three 2D-projections to represent all the feasible solutions in the PF and their correlation. The three previously described MOO algorithms are evaluated under fair conditions.

Table shows the computational time spent by the algorithms per iteration, the number of feasible solutions verified by a conditional which evaluates its status (in each iteration a total of 1872 runs take place), and the ideal optimal points for the unconstrained and constrained results. From a computational time perspective there is not significant difference between algorithms. However, the MOEAD shows a slightly higher computational time. The addition of the constraint does not impact the evaluation time, yet it improves the ideal points specially for the NSGA-III and UNSGA-III results.

Table 6

Optimization Overview for Unconstrained and Constrained Evaluations

algorithm	NSGA-III	UNSGA-III	MOEAD	units
unconstrained
comp. time	3.45	3.49	3.54	hours per iter
feasible sim.	1586	1580	1608	-
CF	15.082	15.082	17.716	MM USD/year
eff	6.69	6.75	6.495	%
CAPEX	42.707	42.686	42.967	MM USD
constrained
comp. time	3.51	3.50	3.55	hours per iter
feasible sim.	1422	1465	1608	-
CF	16.390	16.208	17.716	MM USD/year
eff	6.70	6.65	6.49	%
CAPEX	42.804	42.780	42.967	MM USD

The best results for the optimal set for each evolutionary algorithm are selected based on their closeness to the ideal optimal point, which is an unfeasible point holding the best values of the optimization objectives. The percentage of variation (equivalent to an uncertainty index) of the selected Pareto optimal points to their ideal optimal point are listed in Table . The idea behind this selection is to provide an Euclidean-like distance from the ideal optimal, with a reasonable variation. The selected points aim to establish operational windows. However, these operating conditions require further dynamic evaluations and an operability assessment[61,62] to guarantee its final applicability.

Table 7

Percentage (%) of Variation from the Ideal Optimal Points for Selection

algorithm	NSGA-III	UNSGA-III	MOEAD
unconstrained results	10	5	0.5
constrained results	6	5	0.5

Unconstrained Results

The points observed in Figure represent the PF solutions obtained by the MOO algorithms. The blue pentagons, gold stars, and purple squares represent the PF set achieved by the NSGA-III, UNSGA-III, and MOEAD algorithms, respectively. The colored squared surface denotes the separation between positive and negative CF. In the case of the NSGA-III and UNSGA-III, a broader exploration in the solution domain is observed. These results have a higher diversity and draw a clear surface. On the other hand, the MOEAD results favor the CF. Indeed, the results show only positive and higher values if compared with the results achieved by the other MOO genetic algorithms. Moreover, 2-D projections are illustrated in Figure . These projections depict a nonlinear relationship between the CAPEX and the ef, which might be influenced by the heat recovery unit in the compression section. A strong concavity point is observed when the ef reaches approximately 6.4%. The other two projections CF vs ef, and CF vs CAPEX appear to be scattered without a solid trend.

Figure 3

PF points using NSGA-III (blue pentagon), UNSGA-III (gold star), and MOEA/D (purple square). The violet semiplane delimits a boundary between the negative and the positive CF values.

Figure 4

2D projections of the Pareto obtained by the evolutionary algorithms NSGA-III (blue pentagon), UNSGA-III (yellow star) and MOEAD (purple square). The red vertical lines (lefter and center panels) delimits both negative and the positive cash flow values. The gray colored points represent the ideal points for their respective method.

Even though MOEAD does not provide a wide range of CF, CAPEX and ef, it clearly improves the CF. Such results define highly profitable operating conditions. On the other hand, this algorithm fails to minimize the CAPEX and maximize the ef when compared with the other two algorithms. The NSGA-III and UNSGA-III generate a set of results with a wider distribution and diversity, which achieve similar optimal results for CAPEX and ef. Nevertheless, these algorithms fail to improve the CF when compared to the MOEAD.

Tables –Table list the selected PF for each MOO algorithm, including the decision variables. Table corresponds to the unconstrained results of NSGA-III. The minimum CAPEX is 43.08 MM $, followed by a maximum ef of 6.4% and CF of 15.08 MM $/Y. The UNSGA-III unconstrained results are in Table . The highest ef, 6.73%, is obtained. In contrast, a higher CAPEX is observed. If comparing it to NSGA-III, the results are similar for CF (15.08 MM $/Y). The top-5 results of MOEAD are in Table . An annualized CF of 17.66 MM $/year, an ef of 6.47%, and a minimum CAPEX of 43.06 MM $ are observed. Notice that in all cases the decision variables Lx tend to reach xUB, and F and ΔT are closer to xLB. Moreover, P and P show a comparable range. On the other hand, R1/2 vary depending on the algorithm.

Table 8

Selected Unconstrained Results from NSGA-III

objectives	case 1	case 2	units
CF	14.68	15.08	MM $/Y
ef	6.40	6.05	-
CAPEX	43.06	42.98	MM $

Table 10

Selected Unconstrained Results from MOEAD

objectives	case 1	case 2	case 3	case 4	case 5	units
CF	17.66	17.66	17.66	17.65	17.65	MM $/Y
ef	6.47	6.47	6.47	6.47	6.49	-
CAPEX	43.06	43.06	43.06	43.06	43.07	MM $

Table 9

Selected Unconstrained Results from UNSGA-III

objectives	case 1	case 2	units
CF	15.08	15.00	MM $/Y
ef	6.73	6.42	-
CAPEX	43.38	43.05	MM $

The QIs for the unconstrained evaluation are shown in Table . In terms of cardinality, MOEAD exhibits the highest error ratio as some of the obtained solutions are not part of the overall nondominated solutions. On the other hand, the uniformity of this solution set has the lowest spacing score, showing a high uniformity of results among the MOO algorithms. NSGA-III and UNSGA-III have high spread meaning that offers a broad range for the objectives. These assertions are observed in the 2D projections of the PFs when comparing each solution set.

Table 11

Quality Indicators from Unconstrained Solution Sets

quality indicator	NSGA-III	UNSGA-III	MOEAD
maximum spread	51.14	44.08	0.51
error ratio	0	0	0.33
spacing	2.12	3.57	0.02

Table shows the performance of the MOO genetic algorithms in terms of the C-indicator. As observed, the elements of the NSGA-III and UNSGA-III solution sets fail to dominate the results of MOEAD. In addition, 22% of NSGA-III results and 26% of UNSGA-III results are dominated by MOEAD. Regarding to the UNSGA-III C-indicator, it shows a dominance of 26% over the solution set from NSGA-III. On the other hand, 3% of NSGA-III set dominates UNSGA-III solution.

Table 12

C-Indicators from Solution Sets

algorithms	NSGA-III	MOEAD	UNSGA-III
NSGA-III	1	0	0.03
MOEAD	0.22	1	0.26
UNSGA-III	0.26	0	1

Constrained Results

The constrained PF solutions correspond to the points in Figure . As in the previous section, the algorithms NSGA-III, UNSGA-III, and MOEAD hold the same visual representation. The outcomes achieved in this setup favor the CF results because only profitable points (positive) are considered. NSGA-III and UNSGA-III demonstrate a broader exploration of the solution domain and improve all optimization objectives.

Figure 5

Pareto points obtained by the evolutionary algorithms NSGA-III (blue pentagon), UNSGA-III (yellow star), and MOEAD (purple square) with positive cash flow constrain.

Influenced by the constraint, an overall lower bound CF of 8 MM $/Y was obtained by NSGA-III. The overall maximum CF is 17.66 MM $/year, and it is obtained by MOEAD. However, this result is achieved at a CAPEX of 43.07 MM $, meaning that higher investments will be required. In general, MOEAD offers a narrow set of optimal solutions in the Pareto surface, and it fails to explore lower CAPEX and higher ef. On the other hand, both UNSGA-III and NSGA-III generate a set of more diverse points. Figure illustrates more detailed solutions with positive CF. In addition, a similar exploration for both algorithms is observed. The UNSGA-III optimal set achieves the highest ef (6.55%) of the heat recovery cycle, and the NSGA-III provides the lowest CAPEX (43.09$ MM).

Furthermore, the 2-D projections of the Pareto surface are shown in Figure . No correlation is observed between ef vs CF and CAPEX vs CF. However, the CAPEX vs ef subplot shows two inflection points at 5.6% and 6.4%. The conflict between these variables, and its nonlinear relationship, remains in the constrained evaluation. When comparing these results with the unconstrained, a gap between 5.8% and 6.1% is observed. A reason for this empty range might be that unfeasible operating conditions are neglected during the optimization routine. The possibility of these solutions not satisfying the CF restriction is also considered.

Figure 6

2D projections of the Pareto obtained by the evolutionary algorithms NSGA-III (blue pentagon), UNSGA-III (yellow star), and MOEAD (purple square) with positive cash flow constrain. The gray colored points represent the ideal points for their respective method.

Tables –Table list the constrained results of each MOO algorithm. The top 3 results of NSGA-III are in Table . A minimum CAPEX of 43.09 MM $, and a maximum efficiency of 6.53% and CF of 15.96 MM $/Y are observed. In contrast to the NSGA-III unconstrained results, higher values of CF are achieved along with a higher CAPEX. Table has the top 2 constrained results of UNSGA-III, denoting an efficiency of 6.55%, a CF of 15.78 MM $/Y, and a CAPEX of 43.14 MM $. Regardless the CF constrain, UNSGA-III has a similar trend for Lx and ΔTdew as the previous results, Lx approximates to xUB and ΔTdew approaches to xLB. When compared to the results of NSGA-III, UNSGA-III offers a narrower operational window for the pump P-101, the heat exchangers E-101/102, and the ratio of coolant flow R1. Table depicts the PF of MOEAD. It achieves a CF of 17.66 MM $/Y, an ef of 6.47%, and a minimum CAPEX of 43.06 MM $. As observed, MOEAD has the narrowest operational window for P-101 and the lowest R1 when compared to the previous algorithms. The decision variables and the objective functions are bounded by the same operational windows as the unconstrained results, meaning that the constraint does not have a strong effect on MOEAD.

Table 13

Selected Constrained Results from NSGA-III

objectives	case 1	case 2	case 3	units
CF	15.94	15.49	15.96	MM $/Y
ef	6.34	6.53	6.40	-
CAPEX	43.09	43.26	43.13	MM $

Table 15

Selected Constrained Results from MOEAD

objectives	case 1	case 2	case 3	case 4	case 5	units
CF	17.66	17.66	17.66	17.65	17.65	MM $/Y
ef	6.47	6.47	6.47	6.47	6.49	-
CAPEX	43.06	43.06	43.06	43.06	43.07	MM $

Table 14

Selected Constrained Results from UNSGA-III

objectives	case 1	case 2	units
CF	15.78	15.56	MM $/Y
ef	6.42	6.55	-
CAPEX	43.14	43.15	MM $

Even though the MOEAD results for CF show the highest profitable scenarios, the set of decision variables construct a narrow operating window which is not appropriate in practical terms. These results get trapped and do not offer a broader perspective of the problem nor show an extensive exploration of the solution domain. Such characteristic could be overcome by re adjusting the hyper-parameters of the optimization algorithm. However, default setups were considered in this work, and the mentioned analysis could be considered for future work.

A broader range of operating conditions, as those attained by the NSGA-III, will permit process adjustments while maintaining the ASU profitability above its threshold. This characteristic will offer the flexibility of changing the evaluated decision variables avoiding failure and other detrimental process conditions.

The QI for each MOO evolutionary algorithm are described in Table . In this section, negative values of the cash flow objective function are neglected, leading to a considerable reduction of the spacing and spread values. Nevertheless, since MOEAD does not provide negative CF values in the previous section, its QIs are not affected by this change. In addition, the error ratio of all the solution sets remain the same.

Table 16

Quality Indicators

quality indicator	NSGA-III	UNSGA-III	MOEAD
maximum spread	8.51	7.98	0.51
error ratio	0	0	0.33
spacing	0.37	0.35	0.02

The C-indicators denoted in Table illustrate the augmentation of MOEAD dominance over the other solution sets. In this context, 40% of the elements of NSGA-III and 32% of UNSGA-III solution sets are dominated by MOEAD. No dominance is presented in the other way. Regarding the performance of UNSGA-III and NSGA-III, 10% of NSGA-III solution set dominates UNSGA-III. However, 29% of UNSGA-III results dominate the solution set from NSGA-III.

Table 17

C-Indicators from Solution Sets

algorithms	NSGA-III	MOEAD	UNSGA-III
NSGA-III	1	0.0	0.10
MOEAD	0.40	1	0.32
UNSGA-III	0.29	0	1

Conclusions

In this work, a modular multifunctional framework was proposed for the MOO of an ASU. The integration between the simulation package PRO/II and a Python environment was performed using the Python COM Interface. The continuous monitoring of the simulation status and the replacement of the original simulation before each iteration evaluation assured reliability of the results. Three objective functions were implemented and optimized. The ef of the heat recovery unit and the annualized CF were maximized while the CAPEX minimized. The results displayed a well-defined 3D Pareto surface, in which the trade-off between the CAPEX and the ef was remarked.

Although the results carried out by the MOEAD algorithm achieved the best values in terms of the annualized CF, the NSGA-III and UNSGA-III offered more diversified solutions, confirmed by the spread and spacing QIs. Based on the obtained QIs, the UNSGA-III showed the overall most competitive performance (maximum spread, error ratio, spacing, and dominance).

The unconstrained results showed that some combinations of decision variables lead to unprofitable results (negative CF) at the PF. For this reason, the constrained formulation explored only solutions with a positive annualized CF. The nondominated solutions for the constrained formulation exhibited improvements in some objective functions at the expense of the spread and spacing QIs, which resulted in a poorer performance. However, the overall benefits of better allocating computational resources for MOO with only feasible scenarios. An advantage of the proposed framework is that other objective functions, such as products composition, greenhouse gas emissions, and water use, could be easily incorporated. Such formulation could bring broader insights regarding the compromise from operational, techno-economic, and environmental perspectives.

For future work, the study of fine-tuning the parameters of the evolutionary MOO algorithms could be considered. These parameters control several functional aspects of the algorithms, such as the probability of the mutation process, alternative crossover techniques, and the number of required parents to create offspring.