## 3 Methodology
### Mathematical Model  
The PSA process is complex because it operates cyclically and reaches a cyclic steady state. The rigorous model is difficult to control. Our solution involves using the input and  
Figure 3: Implemented methodology.  
Figure 2: (**a**) Purities of the products obtained in the PSA process. (**b**) Magnification of the \(\mathrm{H_{2}}\) composition. (**c**) Magnification of the \(\mathrm{CH_{4}}\) composition.  
output data from the rigorous model to create a simplified model using system identification techniques [42]. The model used is a nonlinear Hammerstein-Wiener model, which consists of a linear component and two nonlinear components. The data used in this study were sourced from the work referenced in [23; 43]. Figure 4a shows the identified model and Figure 4b the inverse functions.  
The H-W model combines a linear dynamic model with two static nonlinear blocks at the input and output, designed to capture the rigorous PSA model's dynamics accurately. Unlike linear models, which fail to reproduce the PSA's dynamics correctly, the H-W model was chosen due to its capability to overcome this limitation. In the controller design, the linear dynamic model was applied in conjunction with the inverse functions of the Hammerstein and Wiener blocks, ensuring smooth integration with the detailed model (see Figure 4b). The parameters of the static nonlinear functions and their inverses are detailed in Tables 4 and 5, while the linear dynamic model is described in Equation (9). In Tables 4 and 5, \(X\) and \(Y\) represent the breakpoints of the static nonlinear functions in the Hammerstein-Wiener (H-W) model. These nonlinear functions are used to model the system's nonlinear behavior at the input and output.  
\[\frac{\bar{q}(k)}{\bar{q}(k)}=\frac{-0.3410z^{-1}+sz^{-2}-0.9803z^{-3}+0.3212z^ {-4}}{1-2.9300z^{-1}+2.  
The inverse functions for the H-W blocks were derived using the parameters provided in Tables 4 and 5, as shown in Equations (10) and (11). Moreover, Equation (9), which defines the transfer function, was reformulated into a discrete-time state-space representation, as depicted in Equation (12).  
\[\tilde{u}_{k}=\begin{cases}3.00u(k)+310.20,&55.70<u(k)\leq 56.60,\\ 2.30u(k)+309.10,&54.50<u(k)\leq 55.70,\\ 2.30u(k)+308.90,&53.20<u(k)\leq 54.50,\\ 2.30u(k)+307.70,&51.90<u(k)\leq 53.20,\\ 2.60u(k)+306.60,&50.80<u(k)\leq 51.90,\\ 3.50u(k)+305.20,&49.90<u(k)\leq 50.80.\end{cases} \tag{10}\]  
\[y(k)=\begin{cases}160.55g(k)-182.0,&0.97<g(k)\leq 0.99,\\ 0.019(k)-11.5,&0.96<g(k)\leq 0.97,\\ -539.80g(k)+512.0,&0.88<g(k)\leq 0.96,\\ -302.48g(k)+302.7,&0.76<g(k)\leq 0.88.\end{cases} \tag{11}\]  
The following state-space representation describes the discrete linear system:  
\[\begin{array}{c}x(k+1)=Gx(k)+H\tilde{u}(k)\\ g(k)=Cx(k)\end{array} \tag{12}\]  
\[G=\begin{bmatrix}2.930&-1.446&0.504&0.085&-0.334&0.337&-0.124&0.007&0.039&-0. 113\\ 2&0&0&0&0&0&0&0&0\\ 0&1&0&0&0&0&0&0&0\\ 0&0&1&0&0&0&0&0&0\\ 0&0&0&1&0&0&0&0&0\\ 0&0&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0.500&0&0\\ 0&0&0&0&0&0&0&0.250&0\end{bmatrix},\]  
\[H=\begin{bmatrix}1&0&0&0&0&0&0&0&0\\ \end{bmatrix}^{T},\]  
\[C=\begin{bmatrix}-0.34&0.50&-0.49&0.16&0&0&0&0&0\\ \end{bmatrix},\]  
Let \(\mathbf{x}\in\mathbb{R}^{n\times 1}\) represent the state vector, \(\tilde{\mathbf{u}}\in\mathbb{R}^{m\times 1}\) represent the input vector, and \(\tilde{\mathbf{y}}\in\mathbb{R}^{r\times 1}\) represent the output vector. Consequently, let \(\mathbf{G}\in\mathbb{R}^{n\times n}\) denote the state matrix, \(\mathbf{H}\in\mathbb{R}^{n\times m}\) denote the input matrix, and \(\mathbf{C}\in\mathbb{R}^{n\times r}\) denote the output matrix. In this representation, the sampling period \(k\) is assumed to be \(1\) s.