Simultaneous Optimization of Hydrogen Network with Pressure Swing Absorption Based on Evolutionary Response Surface Method

Abstract

The simultaneous optimization of complex process units and hydrogen networks is a significant challenge in refinery hydrogen network integration. To address this, an evolutionary response surface-based collaborative optimization method is proposed, enabling the concurrent optimization of pressure swing adsorption (PSA) and the hydrogen network. This method develops a mechanistic model for PSA and alternates between random sampling and evolutionary response surface-based hydrogen network optimization to obtain diverse sampling points and potential optimal solutions. The PSA mechanistic model is then used to compute the accurate output parameters for the sampled points, and these parameters are incorporated into the hydrogen network optimization to obtain precise objective function values. An efficient optimization framework is presented to streamline the process. The proposed method is applied to a refinery hydrogen network integration case study, comprehensively considering both PSA costs and hydrogen utility costs. The results demonstrate that the method is computationally efficient and effectively reduces the refinery’s total annual costs. The accuracy of the optimization results is significantly improved compared to traditional methods, providing an effective solution for the collaborative optimization of the refinery hydrogen network and PSA.

1. Introduction

According to the International Energy Agency, global hydrogen demand reached 97 Mt in 2023, marking a 2.5% increase compared to 2022. Demand is expected to further rise to 100 Mt in 2024 [1]. Reducing hydrogen consumption in refineries is critical for both cost savings and environmental sustainability. Hydrogenation reactions in refineries often require additional hydrogen to achieve high conversion rates. Excess hydrogen in the product stream can be recovered through gas–liquid separation processes. However, hydrogen recovered from low-pressure separators is typically of low concentration and contains impurities, limiting its usability and leading to hydrogen losses. To address these challenges, researchers have explored various strategies for designing and optimizing hydrogen networks that incorporate PSA technology.

PSA is a dynamic separation process that is highly sensitive to both capital and operating parameters. Current research predominantly adopts simplified models or surrogate models for process units. Knaebel [2] introduced a simplified PSA model based on linear adsorption isotherms, which expresses hydrogen recovery as a function of adsorbent selectivity, operating pressure, and feed concentration. This approach simplifies the process while retaining key performance characteristics. Researchers have also applied simplified PSA models to inter-plant hydrogen network optimization [3] and hydrogen network retrofit [4]. However, due to the simplified PSA model’s neglect of correlations between product and design parameters, its accuracy is limited.

To balance computational cost and model accuracy, Kriging models can be employed as surrogate models for PSA, incorporating additional design parameters such as adsorption time, pressure, and feed conditions. Li et al. [5] developed a hydrogen network optimization method integrating a Kriging-based PSA model, deriving relationships among hydrogen-saving capacity, hydrogen storage material, adsorption time, and purification feed concentration. This approach enables the identification of the optimal purification feed, maximum hydrogen savings, required hydrogen storage material, and adsorption time. However, surrogate models inevitably introduce errors. When modeling key equipment, even minor errors can significantly impact the final optimization results.

In addition to PSA-based optimization, other researchers have explored innovative approaches to hydrogen network design and integration. Yang et al. [6] proposed a two-stage strategy for hydrogen network synthesis, incorporating a group cascade layout of compressors to address pressure uncertainties and improve energy efficiency. Similarly, Zhou et al. [7] developed a mixed-integer nonlinear programming (MINLP) model for refinery hydrogen network synthesis that integrates detailed performance models for centrifugal and reciprocating compressors to minimize annual costs. These approaches highlight the importance of compressor configurations in hydrogen networks but often assume fixed operating conditions, limiting their flexibility in dynamic scenarios.

Hydrogen network optimization has also been integrated with light hydrocarbon recovery (LHR). Yang et al. [8] employed a hybrid method that combines pinch analysis with rigorous process simulation to optimize hydrogen networks integrated with LHR, significantly reducing fresh hydrogen consumption and improving economic returns. However, this study did not optimize the configuration of LHR processes or fully integrate them with hydrogen network design. Deng et al. [9] proposed a systematic retrofit method that incorporates LHR into hydrogen network optimization using pinch analysis and simulation, achieving significant economic benefits and hydrogen savings. While effective, their method assumed fixed configurations for LHR processes, limiting its adaptability to multi-period and dynamic scenarios.

Efforts to enhance hydrogen network performance through thermodynamic principles have also been explored. Zhang et al. [10] proposed a thermodynamic principle-based optimization method for hydrogen networks, considering sulfur content variations in hydrorefining feed oil to minimize total exergy consumption. This approach demonstrated the potential to reduce hydrogen consumption and compression work but focused only on sulfur contaminants, excluding other reaction categories. Additionally, Chang et al. [11] introduced a globally optimal design framework for refinery hydrogen networks using pressure discretization to linearize compression power constraints. However, this method assumed steady-state operating conditions and ignored dynamic fluctuations in hydrogen demand.

Incorporating surplus hydrogen into broader energy systems has also been investigated. Rezaie et al. [12] developed a MINLP-based optimization model for retrofitting refinery hydrogen networks, utilizing surplus hydrogen for electricity generation via fuel cells and turbines, as well as ammonia production, to enhance economic performance. Although innovative, this study relied on simplified surplus hydrogen scenarios and did not consider dynamic fluctuations in hydrogen availability or energy integration strategies. Furthermore, Zhang et al. [13] introduced a simulation and modeling-based optimization methodology for refinery hydrogen network integration, incorporating process risk analysis to minimize total exergy consumption. However, this study focused primarily on flash pressure as a risk factor and assumed fixed operating conditions, limiting its applicability to more complex industrial scenarios.

The problem investigated in this study involves resource network optimization with process unit models, which is a common challenge in chemical engineering. Some studies have addressed the simultaneous optimization of complex processes and hydrogen networks using mechanistic models. For instance, Yang et al. [14] proposed a genetic algorithm-based simultaneous optimization method for light hydrocarbon recovery units (including cryogenic separation, membrane separation, oil-based adsorption, and PSA) and hydrogen networks. This study employed simulation software to establish rigorous light hydrocarbon recovery models and optimized the entire process using genetic algorithms, reducing refinery hydrogen consumption and total annual costs (TAC). However, this work only analyzed the effects of feed stream parameter variations at the equipment inlets without fully considering the relationships between operational and dimensional parameters of reaction and separation units and the hydrogen network. Huang et al. [15] utilized the Bayesian method to optimize a hydrogen network integrated with a PSA mechanistic model, further enhancing refinery-wide benefits. However, their method incorporated the hydrogen network optimization process into the surrogate model, resulting in low iterative search efficiency. Wu et al. [16] employed polynomial response surface techniques to establish a low-complexity surrogate model for hydrogen sulfide removal and simultaneously optimized it with a hydrogen network. They did not account for how the surrogate model’s accuracy is affected by the number of decision variables and sample size.

This study addresses the challenge of optimizing hydrogen networks integrated with PSA, which is critical for improving refinery efficiency and reducing operational costs. The primary motivation is to overcome the limitations of existing optimization approaches that either rely on simplified surrogate models, which may compromise accuracy, or use rigorous mechanistic models, which often lead to high computational costs.

The integration of PSA into hydrogen networks plays a vital role in enhancing hydrogen recovery and reducing losses, especially in refineries where hydrogen is a key resource for multiple processes. Optimizing such networks not only minimizes hydrogen consumption and operational expenses but also improves sustainability by reducing emissions and waste. The evolutionary response surface method (RSM) is particularly suitable for addressing these challenges, as it effectively balances computational efficiency with modeling accuracy, enabling scalable solutions for complex industrial systems.

To tackle this issue, we propose an evolutionary response surface-based method for the simultaneous optimization of PSA and hydrogen networks. The method begins by developing a PSA mechanistic model to generate an initial sample set, which is subsequently used to construct an RSM for PSA output parameters. The RSM is then employed as a surrogate model to facilitate hydrogen network optimization and serves as the acquisition function for iterative point selection. To ensure accuracy, the PSA mechanistic model is used to compute precise output parameter values, refining the RSM-based optimization results and improving iterative efficiency. The proposed approach is validated through a practical case study, demonstrating its effectiveness in reducing computational costs while maintaining high accuracy, thereby addressing the trade-off between computational efficiency and model precision.

2. Background

This study aims to optimize the refinery hydrogen network considering PSA. An evolutionary RSM-based optimization algorithm is proposed. The approach comprises mechanistic modeling, hydrogen network optimization, and RSM development.

2.1. PSA Mechanistic Model

The PSA operates cyclically, alternating between high-pressure adsorption and low-pressure desorption of target gases. This study adopts an operational cycle derived from Li et al. [17], encompassing six phases: adsorption (AD), depressurizing pressure equalization (DPE), blowdown (BD), purge (PG), pressurizing pressure equalization (PPE), and pressurization (PR) (Figure 1). During AD, the feed gas is pressurized into bed I at p₁. Part of the hydrogen-rich gas exits as a product, while another portion purges bed II. The operation duration of this stage is set according to Li et al. [17], which prevents impurity concentration breakthrough. In DPE, pressure in bed I drops to p₂, transferring gas to bed II and raising its pressure to p₃. BD reduces bed I’s pressure concurrently to p₄ (near atmospheric), expelling desorbed waste gas. In PG, residual impurities are purged from bed I with hydrogen from bed II at p₄. Concurrently, PPE pressurizes bed I to p₃ using gas from bed II. Finally, PR raises bed I’s pressure to its adsorption level using feed gas, completing the cycle.

The adsorption tower model is described by a series of governing equations [17]. These include mass balance equations for both individual components (Equation (1)) and the overall system (Equation (2)), a linear driving force model (Equation (3)) for adsorption kinetics, and energy balance equations (Equations (4) and (5)) for both gas and solid phases. The Ergun equation (Equation (6)) is used to model momentum balance, while an extended multi-site Langmuir model (Equation (7)) captures the adsorbent’s uptake capacity. Key parameters include molar concentrations, heat transfer coefficients, adsorption heat, and specific heat, providing a comprehensive representation of adsorption and transport phenomena.

$${\displaystyle \frac{\partial \left(v_{\mathrm{gas}}m{c}_k\right)}{\partial z}}+{\epsilon }_{\mathrm{B}}{\displaystyle \frac{\partial m{c}_k}{\partial t}}+{\rho }_{\mathrm{ad}}{\displaystyle \frac{\partial w_k}{\partial t}}=0$$

(1)

where ν_gas represents the superficial velocity; mc_k represents molar concentration of k; ε_B represents the voidage of the adsorption bed; ρ_ad represents the density of adsorbent particles; w_k represents the loading of k.

$${\displaystyle \frac{\partial \left(v_{\mathrm{gas}}m{c}_{\mathrm{gas}}\right)}{\partial z}}+{\epsilon }_{\mathrm{B}}{\displaystyle \frac{\partial m{c}_{\mathrm{gas}}}{\partial t}}+{\rho }_{\mathrm{ad}}{\displaystyle \sum _k\frac{\partial w_k}{\partial t}}=0$$

(2)

where mc_gas represents the molar concentration of gas phase.

$${\displaystyle \frac{\partial w_k}{\partial t}}=MT{C}_{s,k}\left(w_k^{\ast }-w_k\right)$$

(3)

where $w_k^{\ast }$ represents the equilibrium adsorbed concentration of k; MTC_s,k represents the linear driving force coefficient.

$$C_{v,\mathrm{gas}}{v}_{\mathrm{gas}}m{c}_{\mathrm{gas}}{\displaystyle \frac{\partial T_{\mathrm{gas}}}{\partial z}}+{\epsilon }_{\mathrm{B}}{C}_{v,\mathrm{gas}}m{c}_{\mathrm{gas}}{\displaystyle \frac{\partial T_{\mathrm{gas}}}{\partial t}}+p{\displaystyle \frac{\partial v_{\mathrm{gas}}}{\partial z}}+HTC{a}_P\left(T_{\mathrm{gas}}-T_{\mathrm{S}}\right)+{\displaystyle \frac{4H_{\mathrm{w}}}{D_{\mathrm{B}}}}\left(T_{\mathrm{gas}}-T_{\mathrm{amb}}\right)=0$$

(4)

where C_ν,g represents gas-specific heat; T represents temperature; p represents pressure; HTC represents heat transfer coefficient between gas and solid; a_P represents specific surface per unit volume bed; H_w represents heat transfer coefficient between gas and wall of adsorption bed; T_amb represents ambient temperature.

$${\rho }_{\mathrm{ad}}{C}_{p,\mathrm{ad}}{\displaystyle \frac{\partial T_{\mathrm{S}}}{\partial t}}+{\rho }_{\mathrm{S}}{\displaystyle \sum _{k=1}^n\left(C_{p,a,k}{w}_k\right)}{\displaystyle \frac{\partial T_{\mathrm{S}}}{\partial t}}+{\rho }_{\mathrm{S}}{\displaystyle \sum _{k=1}^n\left(\Delta H_k\frac{\partial w_k}{\partial t}\right)}-HTC{a}_{\mathrm{P}}\left(T_{\mathrm{g}}-T_{\mathrm{S}}\right)=0$$

(5)

where C_p,ad represents solid specific heat; C_p,a,k specific heat of k; ∆H_k represents the adsorption heat of k.

$${\displaystyle \frac{\partial p}{\partial z}}=-\left({\displaystyle \frac{1.5\times {10}^{-3}{\left(1-{\epsilon }_1\right)}^2}{{\left(2r_P\psi \right)}^2{\epsilon }_1^3}}\mu v_{\mathrm{gas}}+1.75\times {10}^{-5}Mm{c}_{\mathrm{g}}{\displaystyle \frac{\left(1-{\epsilon }_1\right)}{2r_P\psi {\epsilon }_1^3}}{v}_{\mathrm{gas}}^2\right)$$

(6)

where r_P represents particle radius; ψ represents the particle shape factor; μ represents dynamic viscosity; M represents molecular weight.

$$w_k^{*}={\displaystyle \frac{w_{s,k}{b}_k{p}_k}{1+{\displaystyle \sum _k\left(b_k{p}_k\right)}}}={\displaystyle \frac{I{P}_{1k}{p}_k}{1+{\displaystyle \sum _k\left(I{P}_{2k}{p}_k\right)}}}$$

(7)

where w_s,k represents specific saturation adsorption capacity of k; b_k represents the affinity parameter of k; IP represents the isotherm parameter.

Following the establishment of the PSA mechanistic model, it is important to assess the associated costs to evaluate its economic viability. The annual cost of the PSA is related to the capital and operating parameters. The costing method and parameters are adopted from Turton et al. [18] and Susarla et al. [19].

$$A{C}_{\mathrm{PSA}}=ACC+AOC$$

(8)

$$ACC=0.0815\times TPC$$

(9)

$$TPC=1.716\times \left(BM{C}_{\mathrm{B}}+BM{C}_{\mathrm{cool}}+BM{C}_{\mathrm{com}}+BM{C}_{\mathrm{motor}}\right)$$

(10)

$$AOC=CA+613.2H{P}_{\mathrm{com}}+0.055\times TPC$$

(11)

$$CA={\displaystyle \frac{2}{3}}\times UC\times \rho \times \left(1-{\epsilon }_{\mathrm{B}}\right)\times N\times V$$

(12)

where AC_PSA represents the annual cost of the PSA; ACC represents the annualized capital cost; AOC represents annualized operating cost; TPC represents total plant cost; BMC represents bare module cost; B represents heat exchanger area index; com represents compressor; CA represents cost of adsorbent; UC represents unit cost of adsorbent; ρ represents density; ε represents voidage; N represents number of adsorption bed; V represents volume.

The bare module costs for various major equipment in the capture plant are shown in the following equations:

$$BM{C}_{\mathrm{B}}=4.65\times {\displaystyle \frac{586}{397}}\times M\times N\times P{C}_{\mathrm{B}}$$

(13)

$$\log \left(P{C}_{\mathrm{B}}\right)=3.4974+0.4485\log \left(V\right)+0.1074{\left[\log \left(V\right)\right]}^2$$

(14)

$$BM{C}_{\mathrm{cool}}=3.3\times {\displaystyle \frac{586}{397}}\times P{C}_{\mathrm{cool}}$$

(15)

$$\log \left(P{C}_{cool}\right)=3.2138+0.2688\log \left(Area\right)+0.07961{\left[\log \left(Area\right)\right]}^2$$

(16)

$$BM{C}_{\mathrm{com}}=2.5\times {\displaystyle \frac{586}{397}}\times P{C}_{\mathrm{com}}$$

(17)

$$\log \left(P{C}_{\mathrm{com}}\right)=2.2897+1.3604\log \left(E{C}_{\mathrm{com}}\right)-0.1027{\left[\log \left(E{C}_{\mathrm{com}}\right)\right]}^2$$

(18)

$$BM{C}_{\mathrm{motor}}=2.5\times {\displaystyle \frac{586}{397}}\times P{C}_{\mathrm{motor}}$$

(19)

$$\log \left(P{C}_{\mathrm{motor}}\right)=2.1206+0.9545\log \left(E{C}_{\mathrm{com}}\right)-0.0661{\left[\log \left(E{C}_{\mathrm{com}}\right)\right]}^2$$

(20)

where PC represents purchase cost; EC represents energy consumption.

2.2. Hydrogen Network Superstructure Optimization

In refineries, multiple hydrogen-consuming units operate concurrently. Hydrogen streams exiting these units serve as hydrogen sources, while the incoming streams represent hydrogen sinks. These hydrogen sources not only include hydrogen recovered from process units but also hydrogen produced by dedicated hydrogen generation facilities and externally purchased hydrogen. All available hydrogen sources can be routed to meet the demand for hydrogen sinks, forming an integrated hydrogen network, as illustrated in Figure 2. Unused hydrogen sources are routed to the fuel gas system.

To optimize the hydrogen network, a superstructure-based modeling approach is adopted. This superstructure accounts for all possible configurations of hydrogen connections within the refinery, including both direct reuse of hydrogen streams and purification options, such as PSA. In a basic hydrogen network superstructure optimization, variations in equipment parameters are not considered. The primary objective of this optimization is to minimize hydrogen consumption while maintaining sufficient hydrogen quality for the consuming units. The model focuses on key factors such as flow rates and hydrogen concentrations. The resulting optimization framework facilitates the identification of optimal matches between hydrogen sources and hydrogen sinks, with economic objectives taken into consideration.

Without considering equipment parameter changes, the hydrogen network optimization process can be formulated as a nonlinear programming (NLP) problem:

$$\begin{array}{l}\min f\left(F_{i,j}\right)\\ s.t.\\ g_a\left(F_{i,j}\right)=0,a=1,2,\dots ,N_e\\ h_a\left(F_{i,j}\right)\le 0,a=1,2,\dots ,N_i\end{array}$$

(21)

where F_i,j represents the flow rate of hydrogen from source i (i ∈ I) to sink j (j ∈ J); g represents equality constraints, and h represents inequality constraints; N_e and N_i correspond to the number of equality and inequality constraints, respectively. The equality constraints ensure mass conservation between hydrogen sources and sinks, while the inequality constraints ensure that the flow rate from a single hydrogen source to sinks does not exceed the source’s maximum capacity. Given known hydrogen flow rates and concentrations, solving this problem yields the optimal hydrogen network.

2.3. Development of the Response Surface Model

The RSM is a statistical and mathematical tool used to determine the relationship between input and output parameters. In this study, RSM is employed to approximate the complex relationships between PSA parameters and performance outcomes, enabling more efficient optimization. The polynomial response surface expression is presented in Equation (22), while Figure 3 illustrates the construction process of the RSM. The mechanistic model of PSA was developed to evaluate sample points, which were gathered through random sampling. Based on these results, a polynomial response surface was fitted to capture the underlying relationships between the input and output parameters of PSA.

$$y=f\left(x\right)=f_0+{\displaystyle \sum _{k=1}^N{f}_k\left(x_k\right)}+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{k^{\prime }=k+1}^K{f}_{k{k}^{\prime }}\left(x_k,x_{k^{\prime }}\right)}}+\cdots +f_{12\cdots K}\left(x_1,x_2,\cdots ,x_K\right)$$

(22)

where x_k represents the input parameters, y represents the output parameter, f_kk’(x_k, x_k’) represents the second-order terms.

As shown in Figure 4, the input parameters of PSA are length (L) and diameter (D) of the adsorption bed, adsorption time (t_ad), adsorption pressure (p), feed flowrate at the adsorption stage (F_ad), H₂ feed concentration ($c_{{\mathrm{H}}_2}^{\mathrm{in}}$). The output parameters are average feed flowrate (F_PSAin), H₂ product flowrate (F_PSApro), and H₂ product concentration (c_PSApro).

3. Problem Formulation

This study addresses an NLP problem involving a black-box function, specifically focused on hydrogen network planning within the context of a PSA, as outlined in Equation (23). The formulation considers the influence of input parameters on the operation of the PSA. By integrating the PSA model, the approach aims to analyze and optimize the hydrogen network to achieve efficient resource allocation and minimize operational costs.

$$\begin{array}{l}\min f\left(x,F_{i,j}\right)\\ s.t.\\ g_a\left(x,F_{i,j}\right)=0,a=1,2,\dots ,N_e\\ h_a\left(x,F_{i,j}\right)\le 0,a=1,2,\dots ,N_i\end{array}$$

(23)

where x represents the input parameters of the PSA. Compared to Equation (21), Equation (23) introduces additional PSA input parameters as decision variables, and the objective function accounts for PSA cost variations. In Equation (21), the flow rates and concentrations of hydrogen sources and sinks are treated as constants. However, in Equation (23), these parameters for PSA-related hydrogen sources and sinks are functions of the input parameters x. This dependency cannot be described using simple constraint equations, significantly increasing the complexity of solving Equation (23).

To address this challenge, this study proposes an evolutionary RSM. This approach leverages multiple RSMs to approximate the PSA mechanistic model. Equation (23) is formulated as an NLP problem that adds constraints on the input parameters x to Equation (21). After solving Equation (23), the PSA mechanistic model computes precise outputs corresponding to the optimization results, and the RSMs are iteratively updated with newly generated sample points. In this process, the hydrogen network optimization based on the RSM serves as a sampling framework, targeting the identification of potential optimal solutions. The true results are derived from the PSA mechanistic model and are reintegrated into Equation (21) for further refinement. The following section provides a comprehensive explanation of solving Equation (23) using the evolutionary RSM.

4. Methodology

This study aims to solve the hydrogen network optimization problem involving PSA. To improve the efficiency and accuracy of solving this problem, two sample point collection methods (SCM) are proposed in this paper.

As illustrated in Figure 5, the first method, referred to as SCM1, utilizes random sampling to diversify PSA sample points, thereby enhancing the accuracy of the RSMs. The second method, SCM2, optimizes the hydrogen network based on the current RSMs, where the optimal solution obtained serves as the sampling point for further refinement. In SCM2, PSA input–output parameters are approximated using polynomial RSMs, which are integrated as constraints into the hydrogen network optimization problem, forming Equation (23). The decision variables in this model include PSA input parameters and the flow allocation between hydrogen sources and sinks. The primary objective is to identify potential optimal PSA input parameters based on the RSMs. To maintain computational efficiency, the model is solved using the local solver Conopt4 within the GAMS (33, GAMS Development Corporation, Washington, DC, USA) framework.

The SCM2 focuses on exploring the optimal solution, but it may lead to the clustering of PSA sample points, which can affect the effective fitting of the RSM. To balance these two methods, this paper alternates between SCM1 and SCM2 for 20 and 30 iterations, respectively, and incorporates duplicate point checks. If the input parameters obtained by SCM2 are duplicates of existing sample points, the process shifts to SCM1. The total iteration limit is set to 500 iterations.

After obtaining the sample points, the output parameters are computed using the PSA mechanistic model. The resulting input–output parameters are then integrated into the hydrogen network for optimization. The decision variables in this optimization model include only the flow allocation between hydrogen sources and sinks. A global solver, Baron, is used to solve the Equation (21) to obtain accurate optimization results. At the end of each iteration, the optimization results are added to a dynamic result set, from which the optimal solution is extracted after the iterations are completed.

4.1. Random Sampling-Based Hydrogen Network Optimization

As shown in Figure 5, random values for the PSA input parameters are selected for m ≤ 20. Each selected point is checked against the existing sample set to ensure no duplication; if a duplicate is found, a new point is randomly chosen. The selected input parameters are then fed into the PSA mechanistic model to compute the corresponding accurate output parameters. Both the input and output parameters are stored in the sample set. Subsequently, the input parameters obtained from this random sampling and the corresponding output parameters calculated by the PSA mechanistic model are integrated into the hydrogen network. This specific set of parameters is used to optimize the network by solving Equation (21). The optimization results are stored in the dynamic result set.

4.2. Hydrogen Network Optimization with RSMs

For m > 20, the current PSA input–output sample set is used to fit RSMs that approximate the relationship between input and output parameters. In this study, a polynomial RSM, represented by Equation (22), is employed. The polynomial RSMs are then incorporated into the hydrogen network optimization Equation (21), transforming it into Equation (23). Once the optimal input parameters are obtained, a check is performed to ensure that they are not duplicates of points already in the sample set. If duplicates are found, the process reverts to the random sampling procedure described in Section 4.1. If no duplicates are found, the optimal PSA input parameters from this iteration are used in the PSA mechanistic model to compute accurate output parameters. These newly computed input and output parameters are then integrated into Equation (21) to determine the accurate optimal objective function.

Throughout the iterative process, the polynomial RSM adapts and improves as additional PSA input–output samples are incorporated into the sample set, progressively enhancing its accuracy. This ensures that the optimization becomes increasingly precise with each iteration.

4.3. The Constraints of Hydrogen Network

4.3.1. Hydrogen Source and Hydrogen Sink

To fulfill the stipulated criteria for the flow rate (F_j) and hydrogen concentration (c_j) of SK_j, it is imperative that the balance of flow rate and hydrogen concentration, as depicted in Equations (24) and (25), is strictly adhered to. Furthermore, it is essential for each source to conform to the prescribed flow rate constraint detailed in Equation (26).

$${\displaystyle \sum _i{F}_{i,j}}=F_j$$

(24)

$${\displaystyle \sum _i{F}_{i,j}\cdot c_i}\ge F_j\cdot c_j$$

(25)

$$F{\min }_i\le F_i\le F{\max }_i$$

(26)

where Fmin_i and Fmax_i represents the upper and lower limit of sources, respectively.

4.3.2. Fuel Gas System

The hydrogen concentration within the PSA exhaust (F_ex) typically remains at a low level, necessitating its direct transmission to the fuel system. The incoming streams of the fuel system encompass the excess hydrogen from the source, as delineated in Equation (27). Those sources that cannot be aligned with sinks or are subject to purification are directed toward the fuel gas system (FGS), with the corresponding mass balance elucidated by Equation (28).

$${\displaystyle \sum _i{F}_{i,\mathrm{FGS}}}+F_{\mathrm{ex}}=F_{\mathrm{FGS}}$$

(27)

$$F_{i,\mathrm{FGS}}=F_i-{\displaystyle \sum _j{F}_{i,j}}$$

(28)

where F_i,FGS represents the flow rate of source i to the FGS.

4.3.3. PSA

For Equation (23), constraints based on the PSA RSMs are added in this study. The PSA output parameters include the input flow rate of the PSA (F_PSAin), the product flow rate (F_PSApro), and the product concentration (c_PSApro). After fitting RSMs, these parameters can be calculated using Equations (29)–(31), respectively. The PSA inlet stream corresponds to hydrogen sink set (PSAin ∈ J), and the outlet stream corresponds to hydrogen source (PSApro ∈ I). Therefore, Equations (29)–(31) influence the constraints defined in Equations (24)–(26).

$$F_{\mathrm{PSAin}}=f^1\left(x\right)=f_0^1+{\displaystyle \sum _{k=1}^N{f}_k^1\left(x_k\right)}+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{k^{\prime }=k+1}^K{f}_{k{k}^{\prime }}^1\left(x_k,x_{k^{\prime }}\right)}}+\cdots +f_{12\cdots K}^1\left(x_1,x_2,\cdots ,x_K\right)$$

(29)

$$F_{\mathrm{PSApro}}=f^2\left(x\right)=f_0^2+{\displaystyle \sum _{k=1}^N{f}_k^2\left(x_k\right)}+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{k^{\prime }=k+1}^K{f}_{k{k}^{\prime }}^2\left(x_k,x_{k^{\prime }}\right)}}+\cdots +f_{12\cdots K}^2\left(x_1,x_2,\cdots ,x_K\right)$$

(30)

$$c_{\mathrm{PSApro}}=f^3\left(x\right)=f_0^3+{\displaystyle \sum _{k=1}^N{f}_k^3\left(x_k\right)}+{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{k^{\prime }=k+1}^K{f}_{k{k}^{\prime }}^3\left(x_k,x_{k^{\prime }}\right)}}+\cdots +f_{12\cdots K}^3\left(x_1,x_2,\cdots ,x_K\right)$$

(31)

4.3.4. Objective Function

The objective is to minimize the TAC (Equation (32)), which includes the hydrogen utility cost (AC_HP) and the annual cost of the PSA. The unit cost of hydrogen utility is 2.37 $/kmol, and the annual production time is assumed to be 8760 h. Therefore, the total hydrogen cost can be calculated using Equation (33). The annual cost of the PSA depends on its capital and operating parameters. The costing methodology and parameters are adopted from Li et al. [17] and are detailed in Equations (8)–(20).

$$TAC=A{C}_{\mathrm{HP}}+A{C}_{\mathrm{PSA}}$$

(32)

$$A{C}_{\mathrm{HP}}=8760\times 3600\times 2.37\times {10}^{-3}\times F_{\mathrm{HP}}$$

(33)

5. Case Study

This case is derived from the work of Elkamel et al. [20]. The hydrogen distribution system in the refinery is supplied by two sources: a catalytic reformer (CRU) and a hydrogen plant (HP). The hydrogen-consuming units within the refinery include a hydrocracker (HCU), a gas oil hydrotreater (GOHT), a residue hydrotreater (RHT), a straight-run naphtha hydrotreater (NHT), and a diesel hydrotreater (DHT). The source and sink stream data for this case are summarized in

Table 1. Source and sink streams of the hydrogen network [20].

Hydrogen Source	Hydrogen Composition	Flow Rate/mol·s−1	Hydrogen Sink	Hydrogen Composition	Flow Rate/mol·s−1
HCU	0.8000	1156.49	HCU	0.8670	753.46
GOHT	0.7500	1024.49	GOHT	0.8360	685.19
RHT	0.7500	485.43	RHT	0.8260	320.65
DHT	0.7000	154.59	DHT	0.7490	99.21
NHT	0.6500	70.97	NHT	0.7270	47.60
HP	0.9500	1106.89
CRU	0.8000	214.46

. In the initial study, hydrogen purification using PSA was not considered. The TAC was reported to be 82.729 M$/year, with a hydrogen consumption rate of 1106.89 mol·s⁻¹ from the HP. When PSA was integrated into the system, optimization using the method proposed by Li et al. [5] resulted in a TAC reduction to 15.490 M$/year. However, PSA simulations conducted in Aspen Dynamic revealed that the actual TAC, accounting for real concentration and recovery ratios, was 16.700 M$/year. Huang et al. [15] further refined the hydrogen network using a Bayesian method, with the constraints for the PSA of

Table 2. Upper and lower bounds of PSA parameters.

	L/m	D/m	tad/s	p/bar	Fad/mol·s−1·bar−1	cin
Lower bound	3.00	2.00	20.0	10.00	0.100	0.6500
Upper bound	7.20	4.00	80.0	20.00	5.500	0.7000

. This method achieved a TAC reduction to 12.146 M$/year, demonstrating significant improvements in both cost efficiency and optimization accuracy. However, the computation time is relatively long.

5.1. Evaluation of RSMs

The output parameters of the PSA include feed flow rate, product flow rate, and product concentration. Consequently, this study requires the construction of three distinct response surfaces to model these outputs effectively. With seven input parameters for the PSA, second-order, third-order, and fourth-order polynomial response surfaces are employed for optimization. A seven-variable second-order polynomial consists of 36 terms, a third-order polynomial includes 120 terms, and a fourth-order polynomial expands to 210 terms, reflecting increasing complexity with polynomial order.

During the optimization process, non-redundant points—those where polynomial response surfaces are uniquely fitted—total 69 for the second-order polynomial, 56 for the third-order polynomial, and 168 for the fourth-order polynomial. The relative errors between the RSM and the mechanistic model at these points are calculated, as shown in Figure 6.

Figure 6 shows that the relative errors for all three polynomial models (second-order, third-order, and fourth-order) decrease and approach zero as the response surface evolves in the later stages. This indicates that the evolved polynomial response surfaces accurately fit the PSA model, achieving high precision in their predictions.

Figure 7a highlights that the mean relative error of response surface fitting increases with polynomial order. This is mainly due to the increased number of terms in higher-order polynomials, which require a greater number of sample points to accurately estimate their coefficients. Consequently, third- and fourth-order polynomials exhibit higher mean relative errors, particularly in the earlier stages of iteration.

Figure 7b reveals that during the final 10 iterations, the mean relative errors of response surfaces decrease sequentially for second-, third-, and fourth-order polynomials. This observation suggests that while higher-order polynomials demand more samples for accurate fitting, they ultimately achieve greater accuracy in their final response surfaces.

Figure 7 underscores the trade-off associated with higher-order polynomials: although they exhibit larger initial errors and require more sample points to ensure accurate fitting, their errors diminish significantly in later iterations, resulting in highly precise response surfaces. This demonstrates the advantage of higher-order polynomials in capturing complex relationships when sufficient samples are available to support their evolution.

5.2. Optimization Results

Figure 8a presents the results of all iteration points, while Figure 8b shows the results of iteration points collected based on SCM2. During the iterative process, iteration points collected by SCM2 often exhibit a high number of duplicates, leading to an increased number of points randomly collected by SCM1. Consequently, in Figure 8a, the objective function values of the iteration points do not show a clear pattern. However, as seen in Figure 8b, for response surfaces of different polynomial orders, the total number of points collected by SCM2 varies. Nonetheless, in the later stages of iteration, the objective function values tend to stabilize.

Table 3. Optimal PSA parameters obtained by the developed method.

Method	Polynomial Order	L/m	D/m	tad/s	p/bar	Fad/mol·s−1·bar−1	cin	TAC/M$·Year−1	Time/s
RSMs	2	6.65	3.88	78.7	20.00	1.865	0.7294	12.099	22,134
	3	6.65	3.88	78.7	20.00	1.862	0.7294	12.085	24,854
	4	6.65	3.88	78.7	20.00	1.865	0.7294	12.099	23,885
Bayesian		6.54	3.88	78.7	20.00	1.865	0.7294	12.146	95,502

presents the hydrogen network optimization results based on second-, third-, and fourth-order evolutionary RSMs. The results obtained with second- and fourth-order polynomial response surfaces are identical, while the third-order polynomial response surface yields better results. Compared to the Bayesian algorithm, the optimization results of all three polynomial response surfaces outperform the Bayesian method, with significantly reduced computation times. Figure 9 illustrates the optimal hydrogen network matches for the three polynomial response surfaces.

6. Conclusions

This study proposes an evolutionary RSM for the simultaneous optimization of PSA and hydrogen networks, addressing the challenges associated with NLP optimization problems in complex systems. The proposed evolutionary RSM effectively approximates the PSA. In the later stages of iteration, the third- and fourth-order polynomial response surfaces demonstrate high accuracy, maintaining relative errors for feed flow rate, product flow rate, and product concentration within 1%. This ensures that the PSA process is accurately modeled, enabling precise optimization results. The method successfully tackles NLP optimization problems by guiding the selection of iteration points based on hydrogen network optimization results derived from the response surface. It efficiently identifies optimal solutions, outperforming traditional approaches. Compared to the Kriging method, the proposed method achieves a 27.6% reduction in TAC. Furthermore, it demonstrates superior computational efficiency over the Bayesian algorithm, reducing TAC by 0.5% while requiring only one-quarter of the computation time. The study also highlights that higher polynomial orders do not necessarily lead to better results. Although higher-order polynomials require more coefficients to be fitted and thus demand a larger number of sample points, their accuracy does not always justify the added complexity. In particular, the third-order polynomial response surface achieves better optimal solutions compared to both second- and fourth-order surfaces, balancing accuracy and computational efficiency.

In summary, this work validates the proposed method through practical applications, demonstrating its capability to improve computational efficiency and reduce costs while maintaining high accuracy. The findings provide a promising approach for optimizing hydrogen networks integrated with PSA, offering a practical and scalable solution for addressing complex engineering problems.