A Novel Optimization Method via Box-Behnken Design Integrated with Back Propagation Neural Network – Genetic Algorithm on Hydrogen Purification

1. Introduction

Hydrogen, as a highly concerned clean energy in the world, can be used for fuel cells. However, the current production of hydrogen is mainly obtained from fossil fuels currently, resulting in a mixture gas mainly contains impurities such as carbon monoxide, carbon dioxide and nitrogen [1]. There are several methods for hydrogen purification and separation, such as pressure swing adsorption (PSA) [2], temperature swing adsorption (TSA) [2], vacuum swing adsorption (VSA) [2], metal hydride [3,4], organic frameworks [5], and so on.

PSA is an innovative gas separation and adsorption technology that can achieve adsorption and desorption processes through variations in pressure. The conventional process and simple process flow of PSA make it suitable for gas separation and purification in industry and laboratories [6,7,8]. Freund et al. compared the separation technologies of PSA and palladium membrane for hydrogen separation in a polygeneration systems of fuel-rich operated HCCI engines. The results showed that PSA is a superior technology for hydrogen purification [9]. Zhang et al. proposed a two-stage hydrogen purification process based on PSA and CO selective methanation, achieving an effective process for hydrogen purification [10]. Lin et al. developed a modeling and simulation PSA process for CO₂ in natural gas. The results showed that the PSA process can be used for natural gas separation, with methane purity exceeding 95% [11].

It is complicated to achieve good purification performance of PSA process. Including selecting adsorbents, designing the structure of adsorption bed, and determining the number of cyclic operation steps. To assist in the design of PSA process, Isabella et al. simulated the column size and optimized the operating conditions of the process [12]. Yavary et al. compared the adsorption capacity of activated carbon and 5A zeolite granules for multi-component gas mixtures (CO₂, CO, CH₄, and H₂). The adsorption capacity of activated carbon was in the order of CO₂>CH₄>CO>H₂, while the adsorption capacity of zeolite 5A was in the order of CO₂>CO>CH₄>H₂ [13]. Lopes et al. provided the adsorption equilibrium and kinetics of the gases exiting a methane reformer, a layered PSA unit for hydrogen purification was designed and simulated [14]. In addition, they conducted multi-component breakthrough experiments and VPSA cycle experiments under different operating conditions [15]. Ribeiro et al. presented an eight steps PSA process using an activated carbon as an adsorbent, different feed flow rates, humid/dry feed and adiabatic/nonadiabatic operation were considered [16]. Mariem et al. evaluated the adsorption performance of activated carbon and zeolite, the results indicated that the main advantage of activated carbon is its complete reversibility after pressure regulation regeneration, while zeolite exhibits a decrease and stability in adsorption capacity during repeated processes [17]. Reference [18] discussed the purification performance of single- and double-bed systems, and within a certain range, the PSA performance of double-bed system is much better than that of single-bed system. Ribeiro et al.proposed a PSA unit for studying the behavior of both single column and four column systems with layered activated carbon/zeolite beds. The results showed that the introduction of the zeolite layer improved the purity and recovery of the process. In order to investigated the influence of operation conditions on PSA, the effects of feed flow rate, purge to feed ratio and lengths of both adsorbent layers on system performance were also evaluated [19]. Casas et al. found that separation performance depends on operating temperature, adsorption pressure, and desorption pressure [20]. Han et al. proposed a multi-objective optimization method of Bayesian to design steam methane reforming (SMR) reactors and maximize the efficiency of the hydrogen production process [21].

With the increasing maturity and accuracy of PSA technology, as well as the growing demand for high-performance PSA separation, it is crucial to develop optimization strategies for PSA systems. Due to the fact that the PSA process is a periodic dynamic system, optimization remains a challenging task. There are different mathematical models and optimization theories attempting to optimize the PSA process. Rebello et al. proposed a new method for multi-objective optimization of PSA units with optimal Pareto front, which faces a conflict between maximizing productivity and minimizing energy consumption [22]. Subraveti et al. developed a novel steam-purge PSA cycle for pre-combustion CO₂ capture, using genetic algorithm (GA) optimization to minimize energy and maximize productivity [23]. Sankararao et al. used a improved multi-objective simulated annealing (MOSA) and two jumping-gene (JG) approach to study the multi-objective optimization for air separation PSA unit operation. A set of non-dominated solutions was obtained for the PSA unit [24,25,26]. Uebbing et al. developed an optimization strategy that combines a simplified PSA model based on equilibrium theory and can be used to optimize various cyclic adsorption processes [27]. Wang et al. developed a graphical method for optimizing the purification recovery and feed, the results showed that the method can accurately and conveniently determine the optimal hydrogen recovery, optimal purification feed purity, optimal hydrogen utility savings and the pinch points [28].

Several studies have explored the optimization and prediction of hydrogen production based on response surface methodology (RSM). Han et al. proposed a novel hydrogen production process using coke oven gas and optimized it through RSM to achieve maximum annual profit. The cost of hydrogen production process is 1.56$/kg [29]. Tang et al. developed a regression model based on RSM and accurately predicted the optimization objectives [30]. Kombe et al. developed a three-phase numerical simulation model for a fixed-bed gasifier, and the optimal operating conditions can be obtained through RSM analysis [31]. In the preliminary research work, we attempted to optimize the PSA model using Box-Behnken design (BBD) method based on RSM, which can obtain quadratic regression equations for hydrogen purity, recovery, and productivity, thereby significantly saving computation time [32]. However, in extreme states, due to complex nonlinear situations, using only RSM quadratic polynomials in BBD may not achieve the expected prediction accuracy.

As machine learning is widely applied in engineering and science to solve complex problems, the use of artificial neural network (ANN) models can reduce the computational time required for process optimization, and can be used to perform more complex optimization tasks. To explore the performance of ANN models, Tahir et al. compared ANN algorithms, Bayesian regularization, and Scaled Conjugate gradients. The results indicated that the ANN model is most suitable for validation and comparison with process model, demonstrating its accuracy and potential for optimizing biomass gasification processes, with an R₂ value as high as 0.99 [33]. Yu et al. developed an ANN model to optimize the performance of PSA process in producing high-purity hydrogen from steam methane reforming (SMR) gas mixtures. The results showed that ANN can approximate the performance and dynamic behavior of PSA process with extremely high accuracy [34]. Rebello et al. implemented a deep neural network (DNN) model to predict key process performance indicators for CO₂ capture in PSA units. Then, the decision variables corresponding to the Pareto front were tested to verify the effectiveness of the optimization process, and the results showed that DNN achieved convergence more effectively [35]. Streb et al. used ANN as an alternative model for optimizing VPSA, achieving the integration of H₂ purification and CO₂ capture. The results showed good performance in the constraint optimization of H₂ separation [36].

Although ANN has the following advantages: it does not require prior specified appropriate fitting functions; it has universal approximation ability and can approximate almost all kinds of nonlinear functions, including quadratic functions. However, ANN still has some limitations as it cannot guarantee a global optimal solution. To address this limitation, researchers have attempted to combine ANN with other methods to achieve better performance. Kani et al. implemented a hybrid approach in his work [37], using Bayesian and genetic algorithms (GA) to optimize the architecture of ANN models to find the optimal operating point for heat pipes. The results showed that the proposed model can reduce operational and maintenance costs. Samad et al. first developed an ANN model to predict exergy efficiency, and then used the ANN model as an alternative to GA and particle swarm optimization environments to achieve higher exergy efficiency [38]. Numerous studies have shown that optimizing ANN models using heuristic algorithms (such as genetic algorithm) and statistics algorithm (such as Bayesian algorithm) can improve the accuracy and precision of predictions. However, the calculation process of ANN is a “black box” process, which makes it impossible to observe intermediate results. There are several models of ANN. Considering the limitations of RSM and ANN methods, a novel method is proposed in this work. Firstly, the commonly used ANN model of back propagation neural network (BPNN) was chose in this work. A genetic algorithm optimized BPNN structure (BPNN-GA) was designed for the PSA process. Then, a novel optimization model based on BBD RSM was integrated with BPNN-GA (BBD-BPNN-GA) to optimize the hydrogen purification performance of PSA process.

In this work, a dynamic adsorption model for a five-component mixture (H₂/CH₄/CO/N₂/CO₂=56.4/26.6/8.4/5.5/3.1 mol%) flowing through a layered bed packed packed with activated carbon and zeolite 5A was developed in Aspen Adsorption software. Breakthrough curves were validated by comparing with the experimental data [46]. A model for a six-step PSA cycle was implemented, then the optimal parameters of purification performance including the adsorption time, the pressure equalization time and the feed flow rate of PSA process were obtained. Using BBD method, we establish a direct mathematical relation between hydrogen purification performances (hydrogen purity and productivity) and operating parameters, the quadratic regression equations predicted by BBD method are obtained. However, in the limit state, the desired prediction accuracy often cannot be satisfied by only using the BBD quadratic polynomial because of the complex nonlinear situation. In this work, the first objective is to obtain the best neural network architecture based on BPNN-GA model. The second goal is to compare the better approach for optimization condition. The experimental runs generated from BBD are integrated with the BPNN-GA model (BBD-BPNN-GA) for the prediction and optimization of PSA process. Thereafter, optimization of purification conditions is conducted through BPNN-GA and BBD-BPNN-GA models in this work.

2. Materials and Methods

2.1. Mathematical Model for PSA Process

To describe the adsorption process, mass conservation equations, energy conservation equations, momentum equation and adsorption isotherm equation should be written. The mathematical models refer to reference [39].

The mass conservation for the bulk phase in the adsorption column contains tow parts: for the component and overall mass conservation.

$$-{\mathrm{D}}_{\mathrm{L}}{\displaystyle \frac{{\partial }^2{\mathrm{y}}_{\mathrm{i}}}{\partial {\mathrm{z}}^2}}+{\displaystyle \frac{\partial {\mathrm{y}}_{\mathrm{i}}}{\partial \mathrm{t}}}+{\mathrm{u}}_{\mathrm{z}}{\displaystyle \frac{\partial {\mathrm{y}}_{\mathrm{i}}}{\partial \mathrm{z}}}+{\displaystyle \frac{\mathrm{R}\mathrm{T}}{\mathrm{p}}}{\displaystyle \frac{1-{\mathsf{\epsilon }}_{\mathrm{b}}}{{\mathsf{\epsilon }}_{\mathrm{b}}}}{\mathsf{\rho }}_{\mathrm{p}}\left({\displaystyle \frac{\partial {\mathrm{q}}_{\mathrm{i}}}{\partial \mathrm{t}}}-{\mathrm{y}}_{\mathrm{i}}{\sum }_{\mathrm{j}=1}^{\mathrm{N}}{\displaystyle \frac{\partial {\mathrm{q}}_{\mathrm{j}}}{\partial \mathrm{t}}}\right)=0,\mathrm{i}=1,···,\mathrm{N},$$

(1)

$$-{\mathrm{D}}_{\mathrm{L}}{\displaystyle \frac{{\partial }^2\mathrm{p}}{\partial {\mathrm{z}}^2}}+{\displaystyle \frac{\partial \mathrm{p}}{\partial \mathrm{t}}}+\mathrm{p}{\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{z}}}+{\mathrm{u}}_{\mathrm{z}}{\displaystyle \frac{\partial \mathrm{p}}{\partial \mathrm{z}}}-{\displaystyle \frac{\mathrm{p}}{\mathrm{T}}}\left(-{\mathrm{D}}_{\mathrm{L}}{\displaystyle \frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}}+{\displaystyle \frac{\partial \mathrm{T}}{\partial \mathrm{t}}}+{\mathrm{u}}_{\mathrm{z}}{\displaystyle \frac{\partial \mathrm{T}}{\partial \mathrm{z}}}\right)+{\displaystyle \frac{1-{\mathsf{\epsilon }}_{\mathrm{b}}}{{\mathsf{\epsilon }}_{\mathrm{b}}}}{\mathsf{\rho }}_{\mathrm{p}}\mathrm{R}\mathrm{T}{\sum }_{\mathrm{j}=1}^{\mathrm{N}}{\displaystyle \frac{\partial {\mathrm{q}}_{\mathrm{j}}}{\partial \mathrm{t}}}=0,$$

(2)

where ${\mathrm{D}}_{\mathrm{L}}$ is the axial dispersion coefficient; ${\mathrm{u}}_{\mathrm{z}}$ is the axial physical velocity; ${\mathrm{y}}_{\mathrm{i}}$ and ${\mathrm{q}}_{\mathrm{i}}$ represent the molar fraction and the adsorbed phase concentration of species $\mathrm{i}$ respectively; ${\mathsf{\epsilon }}_{\mathrm{b}}$ is the interparticle void fraction; ${\mathsf{\rho }}_{\mathrm{p}}$ is the pellet density of the adsorbent; $\mathrm{R}$ is the universal gas constant, $\mathrm{T}$ and $\mathrm{P}$ correspond to the temperature and pressure in the adsorption bed; $\mathrm{t}$ and $\mathrm{z}$ respectively represent the time since the beginning of the sorption process and the axial position in the bed.

The energy conservation of the PSA system also includes two parts: the energy conservation between the gas and solid phase in the adsorption, the energy conservation at the wall of the adsorption bed.

$$-{\mathrm{K}}_{\mathrm{L}}{\displaystyle \frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}}+\left({\mathsf{\epsilon }}_{\mathrm{t}}{\mathsf{\rho }}_{\mathrm{g}}{\mathrm{C}}_{\mathrm{p}\mathrm{g}}+{\mathsf{\rho }}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{p}\mathrm{s}}\right){\displaystyle \frac{\partial \mathrm{T}}{\partial \mathrm{t}}}+{\mathsf{\rho }}_{\mathrm{g}}{\mathrm{C}}_{\mathrm{p}\mathrm{g}}{\mathsf{\epsilon }}_{\mathrm{b}}{\mathrm{u}}_{\mathrm{z}}{\displaystyle \frac{\partial \mathrm{T}}{\partial \mathrm{z}}}-{\mathsf{\rho }}_{\mathrm{b}}{\sum }_{\mathrm{i}}^{\mathrm{N}}{\mathrm{Q}}_{\mathrm{i}}{\displaystyle \frac{\partial {\mathrm{q}}_{\mathrm{i}}}{\partial \mathrm{t}}}+{\displaystyle \frac{2{\mathrm{h}}_{\mathrm{i}}}{{\mathrm{R}}_{\mathrm{b}\mathrm{i}}}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{w}}\right)=$$

(3)

$${\mathsf{\rho }}_{\mathrm{w}}{\mathrm{C}}_{\mathrm{p}\mathrm{w}}{\mathrm{A}}_{\mathrm{w}}{\displaystyle \frac{\partial {\mathrm{T}}_{\mathrm{w}}}{\partial \mathrm{t}}}=2\mathsf{\pi }{\mathrm{R}}_{\mathrm{b}\mathrm{i}}{\mathrm{h}}_{\mathrm{i}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{w}}\right)-2\mathsf{\pi }{\mathrm{R}}_{\mathrm{b}\mathrm{o}}{\mathrm{h}}_{\mathrm{o}}\left({\mathrm{T}}_{\mathrm{w}}-{\mathrm{T}}_{\mathrm{a}\mathrm{t}\mathrm{m}}\right)$$

(4)

$$A_w=\pi \left(R_{bo}^2-R_{bi}^2\right)$$

(5)

where ${\mathrm{K}}_{\mathrm{L}}$ is the thermal axial dispersion coefficient; ${\mathrm{C}}_{\mathrm{p}\mathrm{g}}$ is the heat capacity of the gas phase; ${\mathrm{C}}_{\mathrm{p}\mathrm{s}}$ is the specific heat capacity of the adsorbent; ${\mathsf{\epsilon }}_{\mathrm{t}}$ is the total void fraction of the bed; ${\mathsf{\epsilon }}_{\mathrm{b}}$ is the interparticle void fraction; ${\mathrm{Q}}_{\mathrm{i}}$ is the heat of adsorption of species $\mathrm{i}$; $h_{\mathrm{i}}$ and $h_{\mathrm{o}}$ are respectively the heat transfer coefficient with inner and outer wall of column; ${\mathrm{T}}_{\mathrm{w}}$ is the wall temperature and ${\mathrm{T}}_{\mathrm{a}\mathrm{t}\mathrm{m}}$ is the ambient temperature; ${\mathrm{R}}_{\mathrm{b}\mathrm{i}}$ is the inside radius of adsorption bed; ${\mathrm{R}}_{\mathrm{b}\mathrm{o}}$ is the outside radius of adsorption bed.

Based on the work of Sereno et al. [40], the model neglects the kinetic energy change in the mechanical energy balance. Momentum in the adsorption bed is determined through Ergun’s equation:

$$-{\displaystyle \frac{\mathrm{d}\mathrm{p}}{\mathrm{d}\mathrm{z}}}=\mathrm{a}\mathsf{\mu }{\mathsf{\upsilon }}_{\mathrm{z}}+\mathrm{b}\mathsf{\rho }{\mathsf{\upsilon }}_{\mathrm{z}}\left|\overrightarrow{\mathsf{\upsilon }}\right|$$

(6)

where the coefficients $\mathrm{a}$ and $\mathrm{b}$ are determined by the following equations:

$$\mathrm{a}={\displaystyle \frac{150}{4{\mathrm{R}}_{\mathrm{p}}^2}}{\displaystyle \frac{{\left(1-{\mathsf{\epsilon }}_{\mathrm{b}}\right)}^2}{{\mathsf{\epsilon }}_{\mathrm{b}}^3}}$$

(7)

$$\mathrm{b}=1.75{\displaystyle \frac{\left(1-{\mathsf{\epsilon }}_{\mathrm{b}}\right)}{2{\mathrm{R}}_{\mathrm{p}}{\mathsf{\epsilon }}_{\mathrm{b}}^3}}$$

(8)

In these equations, $\mathsf{\mu }$ is the dynamic viscosity; ${\mathsf{\upsilon }}_{\mathrm{z}}$ is Darcy’s velocity and ${\mathrm{R}}_{\mathrm{p}}$ corresponds to the particle radius.

The extended Langmuir-Freundlich model is used in this study to express the adsorption isotherms of multi-component gas:

$${\mathrm{q}}_{\mathrm{i}}^{\mathrm{*}}={\displaystyle \frac{{\mathrm{q}}_{\mathrm{m}\mathrm{i}}{\mathrm{B}}_{\mathrm{i}}{\mathrm{p}}_{\mathrm{i}}^{{\mathrm{n}}_{\mathrm{i}}}}{1+{\sum }_{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{B}}_{\mathrm{j}}{\mathrm{p}}_{\mathrm{j}}^{{\mathrm{n}}_{\mathrm{j}}}}},\mathrm{i}=1,\cdots ,\mathrm{N}$$

(9)

In this equation, $q_i^{*}$ is the equilibrium absorbed phase concentration; $q_{mi} \mathrm{a}\mathrm{n}\mathrm{d} B_i$ are the extended Langmuir-Freundlich isotherm parameters and $p_i$ represent the partial pressure of species $i$. The isotherm parameters are functions of temperature:

$${\mathrm{q}}_{\mathrm{m}}={\mathrm{k}}_1+{\mathrm{k}}_2\mathrm{T},\mathrm{B}={\mathrm{k}}_3{\mathrm{e}}^{{\displaystyle \frac{{\mathrm{k}}_4}{\mathrm{T}}}}$$

(10)

$$\mathrm{n}={\mathrm{k}}_5+{\displaystyle \frac{{\mathrm{k}}_6}{\mathrm{T}}}$$

(11)

2.2. BBD Method for PSA Process

RSM is a statistical method used to study the relationship between multi factors and response variables, such as Box-Behnken Design (BBD), Central Composite Design (CCD). Within the range of experimental design, the optimal combination of the factors and responses can be found. Under the same factors, the number of BBD experimental design groups is less than that of CCD, making it more economical. BBD is commonly used for predicting nonlinear models [41,42].

The one factor design was used firstly to determine the range of independent factors for BBD method, each independent optimization parameter is placed at one of three equally spaced values, usually coded as “-1”, “0” and “+1” (the factors). In the BBD method, desirability is an objective function that ranges from zero (outside of the limits of the performance objectives) to one (where the performance objectives are met). A value of one represents the ideal case, and a zero indicates that one or more responses fall outside desirable limits. The overall desirability (D) and the desirability for each response (di) are defined as follows:

$$\mathrm{D}={\left({\mathrm{d}}_1\times {\mathrm{d}}_2\times \cdots \times {\mathrm{d}}_{\mathrm{n}}\right)}^{{\displaystyle \frac{1}{\mathrm{n}}}}={\left(\prod _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{d}}_{\mathrm{i}}\right)}^{{\displaystyle \frac{1}{\mathrm{n}}}},$$

(12)

$${\mathrm{d}}_{\mathrm{i}}={\left({\displaystyle \frac{{\mathrm{y}}_{\mathrm{i}}-\mathrm{L}}{\mathrm{T}-\mathrm{L}}}\right)}^{\mathrm{w}},$$

(13)

where n is the number of responses in the measure, T and L represent the maximum and minimum possible values for the responses, respectively, w is the weight of each response, and yi is the optimum value of each response as determined by the BBD method.

2.3. BPNN-GA for Optimization

ANN is a machine learning method that simulates the structure of the human brain’s nervous system. It is usually used to estimate or approximate functions depend on large number of inputs and generally the functions are unknown. In this work, back propagation neural network (BPNN) is used for prediction and optimization of PSA process, the structure of BPNN is showed in Figure 1. A tangent sigmoid transfer function (tansig) is used for hidden layer, with a log-sigmoid transfer function (logsig) is used at output layer. For training the designed networks, the Levenberg-Marquardt backpropagation (trainlm) is selected. The calculation formula for forward propagation of BPNN can be expressed as follows:

$${\mathrm{X}}_{\mathrm{i}\mathrm{j}}={\sum }_{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{i}-1}}{\mathrm{Y}}_{\left(\mathrm{i}-1\right)\mathrm{k}}\times {\mathrm{W}}_{\left(\mathrm{i}-1\right)\mathrm{k}\mathrm{j}},$$

(14)

$${\mathrm{Y}}_{\mathrm{i}\mathrm{j}}={\mathrm{f}}_{\mathrm{s}}\left({\mathrm{n}\mathrm{e}\mathrm{t}}_{\mathrm{i}\mathrm{j}}\right)={\displaystyle \frac{1}{1+\mathrm{e}\mathrm{x}\mathrm{p}\left[-\left({\mathrm{X}}_{\mathrm{i}\mathrm{j}}-{\mathsf{\theta }}_{\mathrm{i}\mathrm{j}}\right)\right]}},$$

(15)

with ${\mathrm{X}}_{\mathrm{i}\mathrm{j}}$ is the total input of the j-th neuron in the i-th layer; ${\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}$ is the connection weight from the j-th neuron in the i-th layer to the i+1-th layer to the k-th neuron; ${\mathrm{Y}}_{\mathrm{i}\mathrm{j}}$ is the output of the j-th neuron in the i-th layer; ${\mathsf{\theta }}_{\mathrm{i}\mathrm{j}}$ is the threshold of the j-th neuron in the i-th layer; ${\mathrm{N}}_{\mathrm{i}}$ is the number of neurons in the i-th layer.

Although BPNN has strong nonlinear mapping ability, it does not guarantee the global optimal solution. Consequently, the genetic algorithm (GA) is introduced to solve the optimization problem in this work, and it is also cost-effective and less time-consuming technique. The procedure of the BPNN optimized by GA is showed in Figure 2.

3. Results

3.1. Subsection

The PSA process model was build in the previous work [32], the material properties and PSA cycle obtained from the works [43,44,45,46]. A simulation PSA process is build in Aspen Adsorption software as Figure 3 shows.

To validate the simulation PSA model, the adsorption isotherms and breakthrough curves are simulated and compared with the experimental in Figure 4 and Figure 5. Figure 4 shows that the model parameters properly represent the experimental data (parameters and data from [46]) which is an indication that the Langmuir-Freundlich adsorption isotherm model can be effectively used in this work to predict the breakthrough curves. The simulated breakthrough curves were compared with experimental data from Ahn et al [46]. Figure 5 shows the experimental and simulated breakthrough curves of the layered bed with a 0.65 carbon-to-zeolite ratio, a feed flow rate of 8.6 L/min and a pressure of 10 atm during the adsorption step (step I of Table 1). Figure 5 shows overall good agreement between the simulations and the experimental data, validating the model and the use of the extended Langmuir-Freundlich isotherm.

After validating the models, a six-step PSA cycle for a layered bed packed with activated carbon and zeolite 5A was implemented in Aspen Adsorption. The cyclic sequences of the layered bed PSA process are given in Table 1. The simulated PSA process for H₂ purification involved the following cyclic sequence: Ⅰ. High-pressure adsorption (AD); Ⅱ. depressurizing pressure equalization (DPE); Ⅲ. concurrent depressurization (DP); Ⅳ. purge with H2 product (PG); Ⅴ. pressurizing pressure equalization (PPE); and Ⅵ. pressurization with the feed gas (FP) [46]. Figure 6 shows the evolution of the H₂ molar fraction during a PSA cycle. The H₂ product is obtained mainly from the adsorption step.

3.2. Prediction by BBD Method

In this work, the BBD method was implemented in Design Expert™ to correlate the performance of a PSA system to its operating parameters in order to find the optimal operating conditions for hydrogen purification.The adsorption time, the pressure equalization time and the feed flow rate were chosen as independent optimization parameters, hydrogen purity and hydrogen productivity were chosen as the responses of the system. The one factor design was used firstly to determine the range of independent factors for BBD method, each independent optimization parameter is placed at one of three equally spaced values, usually coded as “-1”, “0” and “+1” (the factors), as Table 2 shows. In addition, “-1” as the low value of the independent optimization parameter, “+1” as the high value of the independent optimization parameter and “0” as the center point. With respect to the experimental data as well as the statistical parameters, the BBD method is used to establish quadratic regression equations.

Quadratic regression equations of two responses (hydrogen purity and productivity) were obtained from the analysis:

$$\begin{gathered} \mathrm{p}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{y}=85.26+0.15\mathrm{A}-0.09\mathrm{B}+0.70\mathrm{C}-7.11\times {10}^{-5}\mathrm{A}\mathrm{B}-3.47\times {10}^{-3}\mathrm{A}\mathrm{C}+0.01\mathrm{B}\mathrm{C} \\ -3.59\times {10}^{-4}{\mathrm{A}}^2+7.97\times {10}^{-4}{\mathrm{B}}^2-0.03{\mathrm{C}}^2 \end{gathered}$$

$$\begin{gathered} \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}=9.86-0.21\mathrm{A}+0.23\mathrm{B}+2.43\mathrm{C}+3.28\times {10}^{-4}\mathrm{A}\mathrm{B}+6.98\times {10}^{-4}\mathrm{A}\mathrm{C} \\ -0.03\mathrm{B}\mathrm{C}+5.73\times {10}^{-4}{\mathrm{A}}^2-2.53\times {10}^{-3}{\mathrm{B}}^2-0.05{\mathrm{C}}^2, \end{gathered}$$

where A represents adsorption time, B represents the pressure equalization time and C represents the feed flow rate.

In Figure 7, 3D response surface plots were illustrated by showing the influences of adsorption time, pressure equalization time, and feed flow rate on response variables of purity (a) and productivity (b).

3.3. Prediction by BPNN-GA Model

The BPNN-GA topology utilized for predicting of PSA process. The adsorption time, pressure equalization time, and feed flow rate as the three inputs, with the purity and productivity of hydrogen as the two outputs of BPNN-GA model. The flow chart in Figure 2 illustrates the procedure of BPNN-GA, with the black frame on the right representing the flow chart of BPNN. According to the flow chart, the procedures of the BPNN-GA can be categorized into three parts. Firstly, the structure of BPNN is determined by the inputs and outputs, the weights and bias values of neural network are acquired. Subsequently, a GA is introduced to optimize the values of weights and bias for BPNN, and the optimal initial weights and bias values are obtained. Finally, training and prediction using existing database are conducted for BPNN-GA. The operation process of BPNN-GA primarily involves population initialization, fitness function calculation, and genetic operator execution. The objective of using BPNN-GA is to obtain improve the values of initial weights and bias.

130 sets of data from PSA model are used for the training of the neural network. The datasets is randomly divided into three parts of training set, validation set, and test set with the proportion of 80%: 10%: 10%. The training set is mainly used to enable the neural network to learn and fit the model. The validation set is used to prematurely terminate training when the performance of the network fails to improve or remains stagnant for several epochs. Meanwhile, the test set is employed to evaluate the generalization performance of the network without any effect on training process. Then the test error serves as an evaluation metric for assessing the fit of BPNN.

Figure 8 shows the variation of mean square error (MSE) of training set, validation set, and test set relative to epochs. MSE decreases rapidly at the beginning, and as the number of epochs iterations increase, the rate of decrease slows down. At epoch 5, the minimum validation of MSE is 0.0035, at which point the training stops. The comparison between the actual Aspen values and the predicted BPNN-GA values of purity and productivity is shown in Figure 9. It can be seen that there is a certain deviation between the actual value and predictive value, especially in terms of purity. The R-value is used to measure the correlation between two variables, also known as the correlation coefficient. The closer the correlation coefficient is to 1, the higher the accuracy of BPNN-GA prediction. Figure 10 shows the R-values between the predicted outputs of the BPNN-GA and actual outputs of the Aspen model on the training set, validation set, test set, and overall set. The results indicate that BPNN-GA can be used for the prediction and optimization of PSA process, although there are still some errors.

3.4. Rediction by BBD-BPNN-GA Model

To optimize the predictive performance of the BPNN-GA model, the BBD method of the RSM is introduced in this work. 100 data sets of PSA process are used to predict BBD-BPNN-GA model.

Figure 11 shows the MSE performance of BBD-BPNN-GA mdoel. The best validation performance of MSE is 0.0005 at epoch 12, which is less than the MSE of BPNN-GA model of 0.0035, indicating that the BBD-BPNN-GA model has better prediction performance for the PSA model. Then, the comparison between the actual values of the Aspen model and the prediction values of BBD-BPNN-GA model is shown in Figure 12. It can be seen that compared with the prediction results of the BPNN-GA model, the BBD-BPNN-GA model has better prediction performance and the prediction results are closer to the actual values, which verifies that BBD-BPNN-GA can accurately predict PSA performance. The R-value of the BBD-BPNN-GA model is shown in Figure 13. As we can see, the R-values of the training set, validation set, test set, and overall set all exceeded 0.99. The results also indicate that the predictive performance of the BBD-BPNN-GA model is superior to that of the BPNN-GA model.

4. Conclusions

In this work, breakthrough curves of a five-component mixture (H₂/CH₄/CO/N₂/CO₂ = 56.4/26.6/8.4/5.5/3.1 mol%) flowing through a layered bed packed with activated carbon and zeolite 5A were designed and implemented. The model predictions were found to be in good agreement with experimental data. A 6-step PSA cycle based on the layered bed was implemented in Aspen Adsorption to estimate its performance for hydrogen purification using purity, recovery, and productivity as response functions. To simplify the calculation process, the quadratic regression equations between hydrogen purification performances (hydrogen purity and productivity) and operating parameters (adsorption time, pressure equalization time, and feed flow rate) were established using BBD method. Considering the limitation of BBD method, this study proposed BBD-APNN-GA combined method to optimize the purification performance. Firstly, the best neural network architecture based on BPNN-GA model was obtained. Then, the experimental runs generated from BBD are integrated with the BPNN-GA model (BBD-BPNN-GA) for the prediction and optimization of PSA process. Comparing the optimization results of the two approaches, the results indicated that the BBD-BPNN-GA model have a better performance with the MSE of 0.0005 and R-values are much closer to 1, which is better than that With the MSE of BPNN-GA model is 0.0035. It is also demonstrated that the BBD-BPNN-GA model can be effectively used in the prediction of PSA purification performance.