Simultaneous Integration of the Methanol-to-Olefin Separation Process and Heat Exchanger Network Based on Bi-Level Optimization

Abstract

The separation section of the methanol-to-olefin (MTO) process is energy-intensive, and the optimization and heat integration can enhance energy efficiency and reduce costs. A bi-level optimization model framework is proposed to optimize the separation process and simultaneously integrate the heat exchanger network (HEN). The upper level employs a data-driven BP neural network proxy model instead of the mechanism model for the separation process, while the lower level adopts a stage-wise superstructure for the HEN without stream splits. The interaction between the two systems is realized effectively through information exchange. A bi-level particle swarm algorithm is employed to optimize complex problems and determine the optimal operational parameters for the distillation system and HEN. Compared with the typical sequential synthesis method, the optimization by the proposed approach reduces the total annual cost by $1.4293\times {10}^6$ USD/y, accounting for 4.76%.

1. Introduction

Ethylene and propylene, recognized as the “foundation of the chemical industry” [1], can be applied for synthesizing various crucial derivatives [2]. The demand for low-carbon olefins in the modern chemical industry is steadily increasing, while the supply in the global market is insufficient. Numerous promising pathways have been developed for olefin production. Economic evaluations have been conducted for 20 olefin production routes, encompassing feedstocks such as coal, petroleum, natural gas, and biomass (including steam cracking, propane dehydrogenation, methanol-to-olefins, etc.). The results indicate that the methanol-to-olefin (MTO) route utilizing fossil fuels as a source of methanol is cost-competitive [3]. The MTO reactor’s effluent is a gaseous mixture, and pure ethylene and propylene are separated with light (CO, H₂, CH₄) and heavy components (C₂H₆, C₃H₆, and C₄₊) removed. The recovery process for light components is complex and energy-intensive. Traditional methods, such as cryogenic separation, incur high investment and operational costs because of expensive refrigeration cycles [4,5]. Recognizing the lower light fraction content in MTO process products, a prefractionation and oil absorption separation process was proposed to avoid the application of cryogenic distillation. The lower the light fraction content in the MTO reactor’s effluent, the higher the benefits derived from absorption technology [6,7].

Light olefin recovery consumes a substantial amount of energy. A rational heat integration can partially recover heat and thereby reduce production costs. In the late 1970s, pinch technology was proposed for synthesizing heat exchanger networks (HENs), which developed rapidly [8]. Additionally, mathematical programming methods based on hierarchical superstructures were developed and employed to design HENs [9]. The corresponding optimization can be either a nonlinear programming (NLP) or mixed-integer linear/nonlinear programming (MILP/MINLP) problem, depending on the specific system and target. The superstructure model can cover all possible solutions and is widely applied in process systems engineering, and algorithms have been developed to solve models [10,11].

Sequential synthesis is generally applied in chemical process design, i.e., reaction and separation sections are optimized first, followed by HEN synthesis. This method reduces the difficulty in problem-solving by decomposing the problem into two sub-problems and results in suboptimal solutions because it neglects the correlation between the process and the HEN [12]. For a practical process, simultaneous optimization of chemical processes considering HEN integration is crucial for achieving significant economic and environmental benefits.

The core of optimizing chemical processes is integrating different types of units and sections, considering the energy requirements, and configuring the process to meet these requirements. Since Papoulias and Grossmann [13] proposed the simultaneous optimization strategy, multiple studies have been conducted on synthesizing chemical processes with heat integration. Zhang et al. [14] proposed a method that combines pinch analysis with mathematical programming for systematically integrating the reactor and threshold HEN, determining the optimal conversion rate of the reactor, reactor temperature, and minimum heat transfer temperature difference. Based on the trans-shipment-based HEN model, a MINLP model has been introduced for the simultaneous optimization of chemical processes and the HEN with unclassified cold/hot process streams. The stream’s inlet/outlet temperatures are divided into “dynamic” temperature intervals so that the heat load at each interval can be appropriately calculated [15,16]. Considering the computational difficulties in integrating chemical processes and the HEN, artificial neural networks are used as a computationally efficient alternative to training mechanism models with complex dynamics, aiming to improve computational efficiency [17].

For distillation systems, the operating parameters of each column (pressure and reflux ratio) affect not only its separation performance but also the temperatures and loads of the condenser and reboiler, as well as the outlet products’ temperatures. Integrating the separation system with the background process significantly affects energy consumption. Adjusting the pressure or reflux ratio can change the distillation column’s relative position to the background process’s ground composite curve, achieving better heat integration performance [18]. Zhang and Liu [19] analyzed the impact of column pressure on the composite curve and proposed a graphical method and rules for HEN integration considering changes in distillation column pressure. However, only the distillate and bottom products’ temperature variations are considered, while the changes in the condenser and reboiler’s load are overlooked. On this basis, Duan et al. [20] incorporated the loads of the condenser and reboiler into the composite curve and further analyzed the impact of column pressure variations on the utility consumption of the entire plant. The studies mentioned above took utility consumption as the evaluation index and did not consider the variations in heat exchange area and capital costs. A superstructure optimization method [21] was applied to optimize the olefin distillation separation system, including separation sequence synthesis, column pressure optimization, and heat integration, and the harmony algorithm was proposed to solve the model. A T-Q diagram was used to target the match between the condenser and reboiler, but it could not guarantee that the structure was the most cost-effective.

For bi-level optimization [22], separators’ parameters should be determined on the upper level considering reasonable process plans and production requirements, and they affect the temperatures and flow rates of streams to be integrated into the HEN. In lower-level optimization, the HEN structure is optimized to maximize heat recovery or minimize costs for each set of process parameters given at the upper level. The upper-level parameters play a decisive role in optimization, and the optimization of the lower-level HEN is also crucial. With minimizing the entire system’s total annual cost (TAC) taken as the objective function, the simultaneous optimization of two levels is essential in determining the system’s optimal performance and reducing total annual costs.

Bi-level optimization involves two nested problems with nonlinear, non-convex, and discontinuous characteristics. Heuristic algorithms commonly search a given complex space to target optimal or satisfactory solutions. Commonly used heuristic algorithms include genetic algorithms [21], particle swarm algorithms [23], simulated annealing algorithms [24], and differential evolution algorithms [25]. With improved computer performance, nested evolution algorithms are more prevalent in solving complex bi-level optimization problems. Two evolution algorithms were combined to optimize a non-isothermal reactor network, using a simulated annealing algorithm at the upper level and a particle swarm algorithm at the lower level [26]. The excellent performance of heuristic algorithm-based bi-level optimization in solving HEN synthesis problems has been reported [27,28,29], and the basic idea is to optimize binary integer variables at the upper level and continuous variables, such as the heat load and split ratio, at the lower level. It is worth mentioning that single-level optimization algorithms can be used for the HEN without stream splits [30]. The main difficulty in the simultaneous optimization of the process and HEN design lies in the changes in process parameters. Bi-level optimization can cope with this challenge well. For bi-level optimization problems, the lower-level model is optimized with the upper-level parameters fixed. The optimal solution, together with the upper-level variables, is a feasible solution for the bi-level optimization problem. Among different optimization techniques, bi-level optimization can consider the separation and HEN parameters, characteristics, and their interaction, allowing for better focus on the specific requirements of each problem.

In this study, the optimization of distillation columns and the HEN of the MTO process will be studied, aiming to reduce energy consumption and cost. A method based on the bi-level optimization model will be proposed to synthesize the separation process and the HEN simultaneously. The separation process will be optimized at the upper level and the HEN structure at the lower level. Both processes’ parameters and detailed HEN structures will be obtained simultaneously. This manuscript is organized as follows: the separation system of MTO and the problem to be solved will be introduced in the second section. In the third section, the BP neural network proxy model of the MTO separation system and the stage-wise superstructure model of the HEN without stream splits will be introduced, and the method will be presented for solving the bi-level optimization problem. The application of the proposed method and a comparison of the results with the traditional sequential synthesis approach will be presented in the fourth section, and the last section will contain the conclusions.

2. Methanol-to-Olefin Separation System and Problem Statement

2.1. Methanol-to-Olefin Separation System

The methanol-to-olefin (MTO) separation system separates the reactor effluent, yielding ethylene and propylene products with a certain purity. The reaction gas primarily contains H₂, C₁, C₂, C₃, C₄, C₅₊, and others. The pre-depropanization separation system is employed, which mainly consists of distillation columns, as shown in Figure 1. The reaction gas is partially condensed and separated into liquid and gas streams, which are then introduced into the high-pressure depropanizer (T101H) from different positions. To prevent the polymerization and fouling of diolefins components at the bottom of the high-pressure depropanizer, the bottom product, which is separated by the low-pressure depropanizer (T101L) for further removal of C₄₊ components, is allowed to contain a certain amount of C₃ components. The overhead product of the low-pressure depropanizer is sent to the high-pressure depropanizer as reflux. The top product of the high-pressure depropanizer is pressurized and sent to the demethanizer (T102), which adopts a prefractionation and oil absorption technology to remove methane. In terms of raising the cryogenic temperature of the demethanizer’s condenser, the overhead product is allowed to contain a certain proportion of ethylene, which is absorbed by propane in the absorption column (T103) and then sent back to the demethanizer. Methane and other light components are removed at the top of the absorption column. The bottom product of the demethanizer is sent to the de-ethanizer (T104), which separates C₂ from C₃. The ethylene and ethane gas from the top of the de-ethanizer is sent into the ethylene rectification column (T105). The bottom product of the de-ethanizer contains propylene and propane and is sent to the propylene rectification column (T106). The liquid ethylene product is withdrawn as a side product of the ethylene rectification column and transported in gaseous form after passing through the ethylene vaporizer. For the propylene rectification column, the liquid propylene is obtained at the top; part of its bottom product, propane, is sent to the absorption column as an absorbent, while another part is discharged from the system.

2.2. Problem Statement

In the distillation system shown in Figure 1, the column’s operating pressure and reflux ratio influence the temperatures and heat capacity flow rates of the top/bottom products, the reboiling/condensing streams, and the optimal HEN. The methane removal rate of T102 directly affects the operation of T105. If the methane content at the bottom of T102 is too high, methane will concentrate at the top of T105. Since methane is lighter than ethylene, increasing the reflux in T105 does not affect the methane content in the ethylene product. To ensure the ethylene purity and temperature at the top of T105, the operation that could be performed on T105 is to increase the flow rate of the gas discharged from the top, and this will cause a decrease in the recovery of ethylene. Increasing the distillate flow of T102 can reduce the methane content at the column’s bottom. However, this also leads to an increase in ethylene content at the top of T102. To absorb additional ethylene, more propane absorbent is demanded in absorption column T103, thus increasing the loads of T104 and T106. Additionally, the feed streams of the distillation column may need to pass through a cooler, pressure-reducing valve or pump to ensure the stable operation of the distillation column.

Optimizing a single distillation column to target the best operating conditions in a highly interconnected distillation system is challenging, as it is necessary to consider all distillation columns simultaneously. Selecting the operating pressure of each distillation column (p_i) and the distillate-to-feed ratio (α) of the demethanizer as decision variables is crucial for optimizing the entire distillation system. Under different operating pressures, the reflux ratio can be adjusted to meet separation requirements.

Heat integration is crucial for enhancing energy recovery. However, a tricky issue is that the parameters of the process streams (initial temperature, target temperature, and heat capacity flow rate) to be integrated are uncertain before the results come out. In some cases, the properties of some streams (source or sink) are unknown. For example, when the bottom product of the demethanizer needs to be introduced into the de-ethanizer, a pump, cooler, or throttle valve may be required. The properties of this stream are influenced by the bottom pressure of the demethanizer and the feed pressure of the de-ethanizer.

For a distillation sequence, it is necessary to identify the optimal operational parameters for each column (with fixed tray numbers and feed tray positions) and the most efficient HEN, including detailed specifications for structure, heat exchange areas, and utility consumption. The aim is to minimize the annualized total cost of the HEN. In the heat integration, the condensing and reboiling streams of the distillation columns are taken as hot and cold streams, respectively.

3. Bi-Level Optimization Model

Bi-level optimization is a hierarchical approach involving two interrelated levels. The upper-level problem is responsible for the global optimization of the entire system, while the lower-level problem optimizes local details given the parameters determined in the upper level. As shown in Figure 2, the upper level aims to minimize costs while ensuring product quality and output by adjusting operating conditions. The objective of the lower level is to reduce the HEN’s cost, which involves the costs of utilities and heat exchangers. The objective functions of the two levels are interconnected; parameter changes in the upper level will directly affect the scale and cost of the HEN to be optimized at the lower level. Lower-level optimization can reduce energy costs by improving the HEN’s efficiency and providing feedback to the upper level. This allows the upper level to adopt more cost-effective operational decisions, thereby reducing overall costs. Such a bi-level optimization model facilitates the coordination between two levels, with the upper level’s decision providing the lower level’s constraints and the lower level’s performance providing the upper level with a way to minimize the total cost.

3.1. Upper-Level Optimization Model

3.1.1. BP Neural Network Proxy Model

The process shown in Figure 1 is simulated using Aspen Plus software (version: V12.1), which is embedded with a wealth of physical property methods and unit modules and can carry out efficient and accurate process simulation. The Aspen Plus’s ActiveX interface is used to establish a connection with MATLAB programming software (version: R2021a) [31]. However, solving the simulation model is time-consuming as repeated simulation is demanded by the optimization design and sensitivity analysis.

Directly analyzing a complex system can be highly time-consuming or challenging in numerous engineering and scientific applications. In these instances, proxy models, designed to alleviate problem complexity, play a pivotal role and are widely utilized in the analysis and optimization of complex systems, as their simplified versions can effectively reduce computational costs and complexity while retaining key characteristics and behaviors [32]. Constructing a proxy model involves gathering experimental or simulated data about the original system and using these data to train an approximate model. These models can be physically simplified or purely mathematical–statistical models, such as regression models, neural networks, or machine learning algorithms. The BP neural network is a nonlinear model that excels at capturing the nonlinear relationships within a system, demonstrating strong flexibility when approximating the behavior of the original system [33]. Therefore, the data-driven BP neural network proxy model is used to replace the strict mechanism model to reduce the time of invoking Aspen Plus software and solving the optimization problem.

The model of the process shown in Figure 1 is established in Aspen Plus and used to simulate the process under different operating conditions. Two thousand sets of simulation data, including each stream’s data and each device’s energy consumption, are obtained and used to train the neural network. In the neural network model, the input variables include the pressure of each column and the distillate-to-feed ratio of T102, and the reflux ratio is built into each column’s module as a design specification to meet the product requirements under different operating conditions. In addition to the concentration of key components, the output variables also include the inlet and outlet temperature of each heater and cooler and the energy consumption of the equipment; the total number of output variables is 52. In order to improve the prediction accuracy, the output variables are classified and embedded into four neural network models, and the operation parameters of the associated columns are taken as inputs. The first model predicts outputs for T101 and T101L, the second for T102, T103, and T104, the third for T105 and T106, and the fourth for other heat exchangers.

In the neural network training process, two thousand sets of data are divided into training sets, verification sets, and test sets according to the ratio of 0.8, 0.1, and 0.1, and the mean square error (MSE) is used as the loss function to adjust the weight and threshold. Figure 3 shows the performance curves of the training set, verification set, and test set in the training process of the three neural networks. The consistency trend in the three curves indicates the stability of the training process. It suggests that the model can successfully learn the data’s characteristics in the training process and generalize the unknown samples well instead of overfitting the training data.

In order to analyze the prediction performance of the neural network clearly, the trained neural network is used to predict 500 sets of data and compare them with the actual value. Figure 4 shows that the R² of each output data set is greater than 0.998, indicating a high performance of the fitted model. The trained neural network models can be used in bi-level optimization to realize the nonlinear mapping of input to output.

3.1.2. Objective Function

The process consumes electricity and heating and cooling utilities. The upper-level model takes the total production cost as the objective function, including the power cost identified in the upper level and the HEN cost determined in the lower level, as shown in Equation (1).

$$\min \hspace{0.17em}\hspace{0.17em}{F}_u(x_u,x_l)=Cos{t}^{ELE}+F_l(x_u,x_l)=AOT\cdot P^{ELE}\cdot C^{ELE}+F_l(x_u,x_l)$$

(1)

where $x_u$ is the upper-level decision variable. $x_l$ is the lower-level decision variable. $P^{ELE}$ and $C^{ELE}$ are the power and price of electricity, respectively. $F_u$ and $F_l$ are the upper and lower objective functions, respectively. AOT is the annual operating time.

3.1.3. Constraint Conditions

Propylene is used as the refrigerant throughout the prefractionation and oil absorption separation process, and its lowest temperature is −40 °C. To guarantee the heat transfer temperature difference, the minimum temperature of all streams should not be lower than −37 °C. The purity requirements of the target products (ethylene and propylene) are 99.95% and 99.6%, respectively, as shown in Equation (2).

$$\left\{\begin{array}{l}{t}_{T_i}^{top}\ge -37\\ x_{C_2{H}_4}\ge 0.9995\\ x_{C_3{H}_6}\ge 0.996\end{array}\right.$$

(2)

3.2. Lower-Level Optimization Model

3.2.1. Stage-Wise Superstructure of the HEN without Stream Splits

A stage-wise superstructure (SWS) is used to design the HEN. A HEN with NH hot streams and NC cold streams is divided into NK stages and NK = max(NH, NC), as illustrated in Figure 5. The heat exchange matching within each stage is performed, and utility heaters and coolers are placed at the low- and high-temperature ends of the hot and cold streams, respectively. The appropriate utility can be selected automatically based on investment and operating costs. The optimization variables in this model are the heat loads of all possible heat exchangers ($q_{ijk})$ with a total of $NH\times NC\times NK$ units. The assumptions of the model are as follows:

(1)

The inlet and outlet temperatures of the heating/cooling utility are fixed and can only exchange heat with one cold/hot stream;

(2)

The heat capacity flow rates of steams are constants, and the phase change is not taken into account;

(3)

All heat exchangers, including the heater and cooler, are countercurrent, and the overall heat transfer coefficient is constant.

3.2.2. Objective Function

The performance of the HEN directly impacts the total cost of the entire production process. Minimizing the HEN’s cost is taken as the optimization objective, as shown by Equation (3).

$$\min \hspace{1em}{F}_l(x_u,x_l)={\displaystyle \sum _{i=1}^{NH}{\displaystyle \sum _{j=1}^{NC}{\displaystyle \sum _{k=1}^{NK}{C}_{ijk}^{ex}\cdot z_{ijk}+}}}{\displaystyle \sum _{i=1}^{NH}{C}_{cu}^{ex}\cdot z_i^{cu}}+{\displaystyle \sum _{j=1}^{NC}{C}_{hu}^{ex}\cdot z_j^{hu}+C^{op}}$$

(3)

where $z_{ijk}$, $z_i^{cu}$ and $z_j^{hu}$ represent 0 or 1 based on whether heat exchangers, heaters, or coolers exist, respectively. $C_{ijk}^{ex}$, $C_{cu}^{ex}$, and $C_{hu}^{ex}$ are their investment costs. $C^{op}$ is the operating cost.

(1)

Investment cost

The investment cost of the heat exchanger is primarily estimated based on the heat exchange area, which is determined by the heat load, overall heat transfer coefficient, and temperature difference for heat transfer, as shown by Equation (4).

$$C^{ex}=7269\cdot {\left[Q/(U\cdot LMTD)\right]}^{0.65}/n$$

(4)

where $C^{ex}$ denotes the investment cost of the heat exchange device. $Q$ is the heat load. U is the overall heat transfer coefficient. LMTD is the average heat transfer temperature difference for hot and cold streams. n is the equipment depreciation period.

For a countercurrent heat exchanger, when the temperature difference between the hot and cold stream at the two ends of the heat exchanger is equal, the logarithmic average temperature difference cannot be applied. In this work, LMTD is calculated by the Chen approximation method [34], as shown in Equation (5).

$$\begin{array}{c}LMTD={\left[\left(d{t}_1\cdot d{t}_2\right)\cdot \left(d{t}_1+d{t}_2\right)/2\right]}^{1/3}\\ d{t}_1=t_{out}^h-t_{in}^c\\ d{t}_2=t_{in}^h-t_{out}^c\end{array}$$

(5)

where $d{t}_1$ and $d{t}_2$ represent the temperature differences at two sides of the heat exchanger. $t_{in}^h$ and $t_{out}^h$ are the inlet and outlet temperatures of the hot stream. $t_{in}^c$ and $t_{out}^c$ are the inlet and outlet temperatures of the cold stream.

(2)

Operating cost

The annual operating cost depends on the energy consumption and the price of utilities and can be calculated using Equation (6).

$$C^{op}=AOT\cdot \left[{\displaystyle \sum _{i=1}^{NH}{C}_i^{cu}{q}_i^{cu}}+{\displaystyle \sum _{j=1}^{NC}{C}_j^{hu}{q}_j^{hu}}\right]$$

(6)

where $q_i^{cu}$ and $q_j^{hu}$ are the heat loads of cooler i and heater j, respectively. $C_i^{cu}$ and $C_j^{hu}$ denote the unit price of the corresponding cooling and heating utilities.

3.2.3. Constraints

For a HEN, the constraints primarily consist of the heat balance equations for cold and hot streams, heat exchangers, and those within each stage, the non-negativity constraint of heat load, the maximum heat load limits, the heat transfer temperature difference constraint, and the structural constraints of the HEN, indicated as follows by Equations (7)–(14):

(1)

Heat balance for cold and hot streams.

$$\begin{array}{l}C{P}_i^h\left(T_{i,in}^h-T_{i,out}^h\right)={\displaystyle \sum _{k=1}^{NK}{\displaystyle \sum _{i=1}^{NH}{q}_{ijk}}}+q_i^{cu}\\ C{P}_j^c\left(T_{j,out}^c-T_{j,in}^c\right)={\displaystyle \sum _{k=1}^{NK}{\displaystyle \sum _{j=1}^{NC}{q}_{ijk}}}+q_j^{hu}\end{array}$$

(7)

(2)

Heat balance within each stage.

$$\begin{array}{l}C{P}_i^h\left(t_{i,k}^h-t_{i,k+1}^h\right)={\displaystyle \sum _{j=1}^{NH}{q}_{ijk}}\\ C{P}_j^c\left(t_{j,k}^{\mathrm{c}}-t_{j,k+1}^{\mathrm{c}}\right)={\displaystyle \sum _{i=1}^{NC}{q}_{ijk}}\end{array}$$

(8)

(3)

Heat balance for heat exchanger.

$$\begin{array}{l}C{P}_i^h\left(t_{ijk,in}^h-t_{ijk,out}^h\right)=q_{ijk}\\ C{P}_j^c\left(t_{ijk,out}^c-t_{ijk,in}^c\right)=q_{ijk}\end{array}$$

(9)

(4)

Heat transfer temperature difference constraints.

$$\begin{array}{l}\mathrm{d}{t}_{ijk,1}\ge \mathrm{\Delta }{T}_{\min }\\ \mathrm{d}{t}_{ijk,2}\ge \mathrm{\Delta }{T}_{\min }\end{array}$$

(10)

(5)

The non-negativity of heat loads (lower limit of variable).

$$\begin{array}{l}{q}_{ijk}\ge 0\\ q_i^{\mathrm{cu}}\ge 0\\ q_j^{hu}\ge 0\end{array}$$

(11)

(6)

The maximum heat load (upper limit of variable).

Not every pair of cold and hot streams can be matched. According to the heat transfer temperature difference constraint, when the initial temperatures of a pair of cold and hot streams satisfy Equation (12), these two streams cannot exchange heat within any stage of the HEN. In other words, the maximum heat loads for these heat exchangers can only be 0.

$$T_{i,in}^h\le T_{j,in}^c+\mathrm{\Delta }{T}_{\min }$$

(12)

Heat exchange potential exists for a pair of cold and hot streams that violate Equation (12). The heat load of the match does not exceed the two streams’ minimum heat load, as shown in Equation (13).

$$\begin{array}{l}{q}_{ijk}\le \min \left[C{P}_i^h\left(T_{i,in}^h-T_{i,out}^h\right),C{P}_j^c\left(T_{j,out}^c-T_{j,in}^c\right)\right]\\ q_i^{cu}\le C{P}_i^h\left(T_{i,in}^h-T_{i,out}^h\right)\\ q_j^{hu}\le C{P}_j^c\left(T_{j,out}^c-T_{j,in}^c\right)\end{array}$$

(13)

(7)

Structural constraints.

In each stage, a maximum of one match can exist for each stream, as represented by Equation (14).

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{NC}{z}_{ijk}\le 1},\hspace{1em}i=1,\hspace{0.17em}2,\dots ,\hspace{0.17em}NH,\hspace{1em}k=1,\hspace{0.17em}2,\dots ,\hspace{0.17em}NK\\ {\displaystyle \sum _{i=1}^{NH}{z}_{ijk}\le 1},\hspace{1em}j=1,\hspace{0.17em}2,\dots ,\hspace{0.17em}NC,\hspace{1em}k=1,\hspace{0.17em}2,\dots ,\hspace{0.17em}NK\end{array}$$

(14)

where $z_{ijk}$ is an integer variable, and its values 0 and 1 denote the presence and absence of the heat exchanger, respectively.

3.2.4. Initialization of the Head Loads and Restructuring of the Infeasible Solution

Initialization of the Heat Loads

In the superstructure of the HEN, there is a strong correlation between heat exchangers. The loads of the upstream heat exchangers directly influence the inlet temperatures of the downstream ones and the maximum number of heat exchangers, while these correlations are ignored in Section 3.2.3. The key to optimizing the HEN is determining the matches between process streams and the heat exchange load. According to the superstructure diagram shown in Figure 5, the heat load of each heat exchanger is initialized sequentially in the optimization from the high- to the low-temperature end as follows:

(1) Input the initial data of cold and hot streams, including initial and target temperatures, heat capacity flow rates, and the minimum heat transfer temperature difference.

(2) Initialize the heat loads. From the high- to the low-temperature end, the heat load of each heat exchanger, $q_{ijk}$, depends on the remaining maximum heat exchange load, $q_{ijk,\max }$, and the value of $z_{ijk}$. $q_{ijk,\max }$ and $z_{ijk}$ are determined as follows.

An NC × NH dimensional matrix of the binary parameter, $z_{ij}$, is generated. Each row corresponds to a cold stream, and each column is a hot stream. $z_{ij}$ is set to 0 when the hot stream i and cold stream j satisfy Equation (12). This implies that streams i and j cannot match in principle. The remaining elements are 0 or 1 randomly, but the sum of each column and each row should not exceed 1 to satisfy structural constraints and ensure that there is only one heat exchanger for each stream in each stage.

For a heat exchanger, when the left-side heat transfer temperature difference is less than the allowed minimum value, the heat exchanger does not exist and $z_{ij}=0$. Otherwise, the heat exchanger probably exists, and its maximum load is related to the two streams’ heat capacity flow rates. The inlet temperature of the downstream heat exchanger will be affected by the heat load of the upstream heat exchanger, which will then affect the maximum heat load of the downstream heat exchanger. Therefore, the heat load of the heat exchanger should be initialized one by one from the high-temperature to the low-temperature end. The corresponding calculation formulas are shown in Equations (15) and (16).

$$\begin{array}{l}i{q}_{ijk,\max }=0\\ \mathrm{w}.\mathrm{r}.\mathrm{t}.\hspace{1em}{t}_{ijk,in}^h-t_{ijk,out}^c<\mathrm{\Delta }{T}_{\min }\\ q_{ijk,\max }=\min \left[C{P}_i^h\left(t_{ijk,in}^h-T_{i,out}^h\right),\hspace{0.33em}C{P}_j^c\left(t_{ijk,out}^c-T_{j,in}^c\right)\right]\\ \mathrm{w}.\mathrm{r}.\mathrm{t}.\hspace{1em}{t}_{ijk,in}^h-t_{ijk,out}^c\ge \mathrm{\Delta }{T}_{\min }\hspace{1em}\&\hspace{1em}C{P}_i^h\ge C{P}_j^c\\ q_{ijk,\max }=\min \left[C{P}_i^h\left(t_{ijk,in}^h-T_{i,out}^h\right),\hspace{0.33em}C{P}_j^c\left(t_{ijk,out}^c-T_{j,in}^c\right),\hspace{0.17em}\frac{t_{ijk,in}^h-t_{ijk,out}^c-\mathrm{\Delta }{T}_{\min }}{1/C{P}_i^h-1/C{P}_j^c}\right]\\ \mathrm{w}.\mathrm{r}.\mathrm{t}.\hspace{1em}{t}_{ijk,in}^h-t_{ijk,out}^c\ge \mathrm{\Delta }{T}_{\min }\hspace{1em}\&\hspace{1em}C{P}_i^h<C{P}_j^c\end{array}$$

(15)

$$\begin{array}{l}{q}_{ijk}=rand\cdot q_{ijk,\max }\cdot z_{ijk}\\ t_{ijk,out}^h=t_{ijk,in}^h-q_{ijk}/C{P}_i^h\\ t_{ijk,out}^c=t_{ijk,in}^c-q_{ijk}/C{P}_j^c\end{array}$$

(16)

During the initialization process outlined above, it is crucial to identify the inlet temperature of the hot stream and the outlet temperature of the cold stream in the first stage. In the first stage, the inlet temperature of the hot stream corresponds to the initial temperature. However, there is a heater on the cold stream with an unknown heat load, and it is imperative to assign a random heat load to this heater. As a result, the inlet temperature of the hot stream and the outlet temperature of the cold stream in the first stage can be shown by Equation (17).

$$\begin{array}{l}{q}_j^{hu}=rand\cdot C{P}_j^c\left(T_{j,out}^c-T_{j,in}^c\right)\\ t_{i,1}^h=T_{i,in}^h\\ t_{i,1}^c=T_{j,out}^c-q_j^{hu}/C{P}_j^c\end{array}$$

(17)

Some heat exchangers with small loads will inevitably be generated in this process. These heat exchangers can be redefined with a load of 0 to reduce the complexity of the network.

(3) Calculate the heating and cooling utility consumptions. Suppose a random load is initially allocated to the heater at the end of the cold stream. It is challenging to ensure that the inlet temperature of the cold stream in the last stage equals the initial temperature. A straightforward adjustment is necessary to eliminate this equation constraint. For any given hot stream, if the total heat exchange load determined in Step 2 is less than the demanded load, a cooler can be set to achieve the final target temperature, and its duty, $q_i^{cu}$, can be calculated based on the energy balance. In order to calculate $q_j^{hu}$ without treating it as a variable, the calculation order is adjusted. $q_j^{hu}$ is first initialized to zero, and then $q_{ijk}$ is initialized. Subsequently, $q_{ijk}$ is kept unchanged, and the stream’s temperatures after traversing each heat exchanger are solved incrementally. For any intermediate heat exchanger, the calculation of the cold and hot streams’ outlet temperature is based on the initial temperature; that is, the hot streams’ outlet temperature is calculated from the high-temperature to the low-temperature end, and the cold stream, vice versa. The temperatures of the hot and cold streams outlet from the heat exchanger with a heat load of $q_{ijk}$ are calculated by Equation (18) and Equation (19), respectively.

$$t_{ijk,out}^h=\left\{\begin{array}{l}{T}_{i,in}^h-{\displaystyle \sum _{k=1}^k{\displaystyle \sum _{j=1}^j{q}_{ijk}/C{P}_i^h,\hspace{1em}k=1}}\\ T_{i,in}^h-\left({\displaystyle \sum _{k=1}^{k-1}{\displaystyle \sum _{j=1}^{NC}{q}_{ijk}+}}{\displaystyle \sum _{k=k}^k{\displaystyle \sum _{j=1}^j{q}_{ijk}}}\right)/C{P}_i^h,\hspace{1em}1<k\le NK\end{array}\right.$$

(18)

$$t_{ijk,out}^c=\left\{\begin{array}{l}{T}_{j,in}^c+{\displaystyle \sum _{k=NK}^{NK}{\displaystyle \sum _{i=i}^{NH}{q}_{ijk}}}/C{P}_j^c,\hspace{1em}k=NK\\ T_{j,in}^c+\left({\displaystyle \sum _{k=k}^k{\displaystyle \sum _{i=i}^{NH}{q}_{ijk}}}+{\displaystyle \sum _{k+1}^{NK}{\displaystyle \sum _{i=1}^{NH}{q}_{ijk}}}\right)/C{P}_j^c,\hspace{1em}1\le k<NK\end{array}\right.$$

(19)

If the heat transfer temperature differences at both sides of a heat exchanger are less than $\mathsf{\Delta }T$, the heat exchanger’s load is set to zero. The calculation is reiterated until all heat exchangers comply with the minimum heat transfer temperature difference constraint. Consequently, the outlet temperature $t_{i,NK+1}^h$ for the hot stream in the last stage and the outlet temperature $t_{j,1}^c$ for the cold stream in the first stage can be determined. The heat loads for the cooler and heater are subsequently obtained by Equation (20).

$$\begin{array}{l}{q}_i^{\mathrm{cu}}=C{P}_i^h(t_{i,NK+1}^h-T_{i,out}^h)\\ q_j^{hu}=C{P}_j^c(T_{j,out}^c-t_{j,1}^c)\end{array}$$

(20)

Initializing $q_j^{hu}$ to zero and adjusting the calculation process to remove equality constraints can improve the quality of the solution and benefit heat recovery. Meanwhile, the adjustment also provides convenience for restructuring the infeasible solution of the following text.

Restructuring of the Infeasible Solution

In HEN optimization, q_ijk undergoes random changes based on the update rules, and the newly generated one is prone to violating constraints and, hence, is infeasible. These infeasible solutions are often dealt with by the penalty approach, which diminishes the impact of infeasible solutions by adding a significant value to the calculated objective value. The penalty approach is efficient for optimization problems with relatively weak constraints and can be applied with the penalty factors set in the optimization. For a complex HEN, numerous constraints are involved, and many infeasible solutions may arise. Applying the penalty approach could diminish the quality of the obtained optimal solution. This issue becomes more pronounced with an increased number of streams, variables, and constraints.

The constraints most prone to be violated are the non-negativity of the utility consumptions and the structural constraints of the HEN. The reason is that a process stream’s cumulative heat exchange load exceeds its practical value. Based on this, this paper proposes a method to restructure the infeasible solutions to enhance the quality of the optimal solution.

(1)

Restructure ${\mathit{q}}_{\mathit{i}}^{\mathit{c}\mathit{u}}$ and ${\mathit{q}}_{\mathit{j}}^{\mathit{h}\mathit{u}}$: For the cold and hot streams violating the constraint, the heat loads of heat exchangers are eliminated or reduced in ascending order. Firstly, the utility consumption, $\mathrm{\Delta }{q}_u$, is calculated. If $\mathrm{\Delta }{q}_u<0$ and the minimum heat load ($q_{ijk,\min }$) among all heat exchangers is less than $\mathrm{\Delta }{q}_u$, the corresponding heat exchanger can be eliminated, and its heat load is reset to 0. The new utility consumption, $\mathrm{\Delta }{q}_u^{*}$, equals $\mathrm{\Delta }{q}_u-q_{ijk,\min }$ and is still negative. If the new $q_{ijk,\min }^{*}$ is greater than $\mathrm{\Delta }{q}_u^{*}$, reducing the heat load of the corresponding heat exchanger by $\mathrm{\Delta }{q}_u^{*}$ is adequate to ensure that the utility consumption is non-negative. Based on these steps, the heat exchange loads are restructured.

(2)

Restructure the structural constraints: For a stream matching with two or more streams at each stage, the heat exchanger in the ascending order of heat load is deleted to ensure that each stream exchanges heat only once and satisfies the constraint on the number of heat exchange occurrences. With the heat loads initialized and the infeasible solutions restructured based on the above steps, the obtained solution satisfies the constraints.

When the HEN is integrated simultaneously with the separation process, the streams’ temperatures are unknown. The heating utility is employed if a cold stream has a high temperature and cannot be matched with the hot stream. Similarly, a cooling utility is implemented if a hot stream has a relatively low temperature and cannot be matched with cold streams. This approach can reduce the number of cold and hot streams participating in the HEN synthesis, thereby improving calculation efficiency. The particle swarm algorithm can be used to optimize the HEN, and the optimization procedure is shown in Figure 6.

3.3. Solution Method

The classical and heuristic algorithms represent two distinct optimization approaches [35]. The classical method, which has its foundations in the principles of convex optimization theory, leverages the principles of Lagrange duality and Karush–Kuhn–Tucker conditions to efficiently solve problems with convex objectives and constraints [36,37]. However, its applicability becomes limited when dealing with non-convexity, non-differentiability, or complex constraint conditions, where it may fail to converge to global optima or even provide feasible solutions. In contrast, heuristic algorithms, such as particle swarm optimization (PSO), provide a general alternative for optimization problems that break the assumptions of convexity and differentiability. Unlike KKT methods, heuristic algorithms are not constrained by strict mathematical requirements and can be adapted to a wide range of problem structures, including those with non-convex and non-smooth characteristics [38]. This adaptability makes them particularly suitable for scenarios characterized by complex, irregular solution spaces that are difficult to analyze using classical methods. The nested method is a common method for solving bi-level optimization problems based on heuristic methods. In this method, the decision variables of the upper problems are solved first, and for each of them, the lower problems are solved. Therefore, the overall quality of a bi-level solution is evaluated by considering upper- and lower-level decision variables. It overcomes the limitations associated with KKT methods and enhances the robustness of problem-solving. Considering its global search characteristics and versatility in handling complex problems, its application in chemical processes and HEN synthesis is reasonable and effective.

Particle swarm optimization (PSO) is an evolutionary algorithm inspired by the behavior of flocks of birds [39]. The collective intelligence of the swarm is powerful. It locates the optimal or approximately optimal solutions within a search space based on simulating the foraging behavior of biological populations.

In a swarm containing N particles, particle i (where $i=1,2,\dots ,N$) represents a potential solution to the problem. Each particle possesses the following two main attributes: position in the search space represented by vector $X_i=\left\{x_1,x_2,\dots ,x_d\right\}$ and velocity represented by vector $V_i=\left\{v_1,v_2,\dots ,v_d\right\}$. d denotes the problem’s dimensionality. These particles roam within a given search space to seek the optimal solution. Each particle tracks and records its best-known position (personal best solution, pbest_i) and the global best-known position in the entire swarm (global best solution, gbest). Initially, the positions and velocities of all particles are randomly initialized. In each iteration, each particle updates its velocity and position based on pbest_i and gbest. The updated equations are shown by Equations (21) and (22).

$$V_i^{t+1}=\omega \cdot V_i^t+c_1\cdot ran{d}_1(pbes{t}_i -X_i^t)+c_2\cdot ran{d}_2(gbest-X_i^t)$$

(21)

$$X_i^{t+1}=X_i^t+V_i^{t+1}$$

(22)

The update process is influenced by three factors as follows: the particle’s velocity (inertia term), the particle’s historical experience (individual cognition), and the entire group’s collective experience (social cognition). The particles iterate continuously, updating their velocity and position to approach better solutions gradually. The search process terminates after a certain number of iterations or reaching a specific condition.

4. Optimization of the Olefin Separation System

For the olefin separation process of a practical methanol-to-olefin plant, the data of the mixed olefins sent to the separation system are listed in

Table 1. The data for the mixed olefins.

Components	Mole Fraction (%)
	Vapor	Liquid
H2	4.51	0.1
N2	0.38	0.01
CO	0.34	0.02
CH4	5.4	0.67
C2H4	57.15	25.85
C2H6	1.16	0.75
C3H6	26.41	48.51
C3H8	2.02	4.18
C4H8	2.45	16.04
C5H10	0.18	3.87
Flowrate (kmol·h−1)	2414	1128
Temperature (°C)	9.0
Pressure (kPa)	1720

. The target products are ethylene and propylene; their purities are 99.95% and 99.6%, respectively, and their flow rates are 1665 kmol·h⁻¹ and 1187 kmol·h⁻¹. The objective is to optimize the pressures of each column and implement heat integration to minimize the total annual cost (TAC). The TAC comprises heating and cooling utility costs, heat exchanger investments, and electricity costs. Only the electricity consumption of compressors pressurizing the vapor product of the high-pressure depropanizer is considered, as its energy consumption is significantly higher than that of pumps. The compressor’s efficiency is taken as 0.8. The depreciation life of the heat exchangers is set at ten years, and the annual operating time is 8000 h.

The streams that need to be cooled and heated in the MTO process are listed in

Table 2. Hot and cold streams of the HEN.

Stream Properties	Name	ID
Hot stream	T101H’s top vapor to be condensed	H1
	T101L’s top vapor to be condensed	H2
	T102’s top vapor to be condensed	H3
	T103’s top vapor to be condensed	H4
	T104’s top vapor to be condensed	H5
	T105’s top vapor to be condensed	H6
	T106’s top vapor to be condensed	H7
	T101L’s feed	H8
	T102’s feed	H9
	T103’s propane absorbent	H10
	T103’s mid-pump around	H11
	T104’s feed (uncertain)	H12
	T106’s feed (uncertain)	H13
Cold stream	T101H’s bottom liquid to be reboiled	C1
	T101L’s bottom liquid to be reboiled	C2
	T102’s bottom liquid to be reboiled	C3
	T104’s bottom liquid to be reboiled	C4
	T105’s bottom liquid to be reboiled	C5
	T106’s bottom liquid to be reboiled	C6
	T105’s liquid ethylene product	C7

. Because of the unknown pressure relationships among T102, T104, and T106, some properties of streams H9, H12, and H13 are uncertain. When the pressure at the bottom product (H12) of T102 is lower than the feed pressure of T104, only a pump is required to boost it. Otherwise, a cooler and throttle valve are necessary to keep the stream liquid, preventing impacts on pipes and T104. In this case, H9 is a hot stream, and the same logic applies to the bottom product of T104 (H13). Since low temperatures can enhance the performance of the absorption column, the target temperatures of the gaseous feed (H3) and absorbents (H4, H10) of T103 are all set at −37 °C. The data on the utilities to be used are listed in

Table 3. The utility prices.

Utility	Inlet Temperature (°C)	Outlet Temperature (°C)	Unit Price (USD∙GJ−1)
Propylene refrigerant	−40	−40	12.17
Propylene refrigerant	−20	−20	7.89
Chilled water	5	15	2.43
Cooling water	32	42	0.354
LP steam	120	120	7.40
Electricity	-	-	16.8

.

4.1. Simultaneous Synthesis

When the distillation sequence and HEN are optimized simultaneously, the optimal TAC is $2.8625\times {10}^7$ USD/y. The upper-level variables and partial parameters of the distillation columns are presented in

Table 4. Column parameters obtained based on the simultaneous optimization method.

	T101H	T101L	T102	T104	T105	T106
Pressure (kPa)	1200	384	3330	1789	1628	2173
Condenser temperature (°C)	−3.6	−12.9	−12.2	−32.6	−36.8	51.9
Reboiler temperature (°C)	54.5	41.2	23.9	46.5	−12.9	64.6
Reflux ratio	0.172	1.500	0.644	0.880	1165.0	19.1
Condenser load (kW)	2430	2498	1256	4000	13,419	75,834
Reboiler load (kW)	5917	2690	7489	10,410	8762	75,797

, and the optimal HEN ($\mathrm{\Delta }{T}_{\min }=5$ °C) is illustrated in Figure 7. Note that the identified $\mathrm{\Delta }{T}_{\min }=5$ °C is less than less than the general one, which is about 10 °C. The reason is that decreasing $\mathrm{\Delta }{T}_{\min }$ can significantly reduce the refrigerant cost, and the reduced cost can compensate for the increase in capital cost.

4.2. Sequationuential Synthesis

The results obtained by the sequential optimization method are compared with those obtained by the proposed simultaneous integration method to validate its effectiveness. When the distillation sequence is optimized without considering heat integration and the utility is utilized to fulfill the energy demands of each stream, the TAC is $3.8298\times {10}^7$ USD/y, the optimal decision variables and partial parameters of the distillation columns are listed in

Table 5. Column parameters obtained based on the sequential optimization method.

	T101H	T101L	T102	T104	T105	T106
Pressure (kPa)	1650	456	2178	1836	1733	979
Condenser temperature (°C)	6.8	−7.7	−19.7	−31.7	−36.8	18.0
Reboiler temperature (°C)	70.0	47.1	3.4	48	−10.7	31.8
Reflux ratio	0.223	1.500	0.178	0.970	1110.4	10.3
Condenser load (kW)	2916	2453	464	4370	13,494	55,394
Reboiler load (kW)	6858	2648	4416	10,850	8931	55,342

. With the minimum temperature difference taken as 5 °C and the HEN integrated using the model introduced in Section 3.2, the TAC of the system is reduced to $3.0054\times {10}^7$ USD/y, and the optimal matching structure is illustrated in Figure 8.

4.3. Discussion of the Results

The utility consumptions, total heat recovery, and corresponding costs of the optimal systems identified by two optimization approaches are compared in

Table 6. Comparison of utility consumption and cost identified by the simultaneous and sequential optimization methods.

	Sequential Optimization	Simultaneous Optimization
Cooling water consumption (kW)	0	54,007.7
Chilled water consumption (kW)	50,781.1	800.2
−20 °C propylene refrigeration (kW)	2711.8	4410.8
−40 °C propylene refrigeration (kW)	24,259.1	22,232.9
LP steam consumption (kW)	75,133.3	77,725.6
Total heat recovery (kW)	19,241.1	38,761.4
Number of heat exchangers	28	27
HEN cost (USD/y)	2.9609 × 107	2.7010 × 107
Electricity cost (USD/y)	0.4456 × 106	1.6156 × 106
TAC (with heat integration) (USD/y)	3.0054 × 107	2.8625 × 107
TAC (without heat integration) (USD/y)	3.8298 × 107	4.2288 × 107

. The results indicate that the HEN structure generated using the proposed method has a greater heat recovery; the optimal HEN’s cost is reduced by $2.5992\times {10}^6$ USD/y, accounting for 8.78%.

Figure 9 shows the heat load accumulation diagram of each stream obtained by the two methods. It visually shows the total heat load demand, heat exchanger amount, and the utility consumed by each stream. T106’s pressure obtained by the simultaneous optimization method is 2.173 MPa, and H7 and C6 are the condensing and reboiling stream of the propylene rectification column, respectively. Although the high pressure leads to greater energy consumption, the outlet steam at the top has a significantly elevated temperature and can provide energy to the reboiler of the de-ethanizer. The propylene product (vapor stream) can be switched to the more cost-effective cooling water.

This reduction in TAC also benefits from the increased pressure in the depropanizer, which causes an increment in the compressor’s outlet temperature and further benefits the heat recovery and heat transfer temperature difference. Hence, the energy of stream H9 is fully utilized. The consumption of −40 °C propylene refrigerant decreases significantly, indicating efficient utilization of the process streams’ energy. However, the increased pressure difference across the compressor generates additional electricity consumption. The TAC eventually decreases by $1.4293\times {10}^6$ USD/y.

For the separation system identified using the simultaneous optimization approach, if all heating and cooling demand is satisfied by the utility, the TAC is $4.2288\times {10}^7$ USD/y, which is even greater than that obtained using the sequential synthesis method ($3.8298\times {10}^7$ USD/y). Such results indicate that although the upper-level solution may not be optimal, it leads to improved performance after optimizing the overall system. Therefore, when complex interactions exist among the various components of the system, it is crucial to consider the optimization of each level comprehensively.

The analysis shows that an integrated hierarchical optimization can be obtained with the heat integration model embedded into the separation process optimization. The upper -level process’s superior parameters enable the synthesis of a high-quality HEN. Conversely, the energy consumption of the lower-level HEN is transferred to the upper level to guide the optimization direction of process parameters. This cyclic information transmission between two levels benefits the optimization.

5. Conclusions

In this work, a bi-level optimization model framework based on the particle swarm algorithm is proposed for simultaneously optimizing the separation section and the HEN of the MTO process. The established data-driven BP neural network proxy model can enhance the efficiency of the distillation system optimization. For the stage-wise superstructure model of the HEN without stream splits, initializing heat loads and renovating infeasible solutions can ensure that each solution satisfies the constraints and improves the quality of the optimal solution. Embedding the heat integration model into the process parameter optimization process can realize the information transformation between the separation section and the HEN and effectively promote optimization. Compared with traditional sequential synthesis methods, the optimal operating parameters of the distillation columns and the HEN identified using the proposed model decreased the TAC by $1.4293 \times {10}^6$ USD/y, accounting for 4.76%.

In the chemical process, the suitable number of trays and feed positions may vary with column pressures; their inappropriate selection affects the separation efficiency and cost of the distillation column. Neglecting these two parameters might lead to the design deviating from the system’s optimum performance. In further research, incorporating the number of trays and feed positions into the optimization will be crucial to obtain a more comprehensive and accurate evaluation of the distillation column’s performance, providing more reliable guidance for future process design and operations.