Reduced-Order Modeling and Control of Heat-Integrated Air Separation Column Based on Nonlinear Wave Theory

Abstract

The process of low-temperature air separation consumes a significant amount of energy. Internal heat-integrated distillation technology has considerable energy-saving potential. Therefore, the combination of low-temperature air separation and heat-integrated distillation technology has led to the development of a heat-integrated air separation column (HIASC). Due to the heat integration and the inherent complexity of air separation, the modeling and control of this process poses significant challenges. This paper first introduces the nonlinear wave theory into the HIASC, derives the expression for the velocity of the concentration distribution curve movement and the curve describing function, and then establishes a nonlinear wave model. Compared to the traditional mechanical models, this approach greatly reduces the number of differential equations and variables while ensuring an accurate description of the system characteristics. Subsequently, based on the wave model, a model predictive control scheme is designed for the HIASC. This scheme is compared with two conventional control schemes: PID and a general model control. The simulation results demonstrate that MPC outperforms the other control schemes from the response curves and performance metrics.

1. Introduction

Air separation is the process of separating the main components of air based on the differences in their physical properties. The oxygen, nitrogen, chlorine, and other gases produced by air separation equipment have wide applications in the national economy [1,2,3,4]. The cryogenic air separation is a primary approach for obtaining high-purity gases like nitrogen, oxygen, and argon. The raw material for cryogenic air separation units is atmospheric air, and their primary energy consumption comes from electricity, accounting for over 70% of the production cost [4,5,6,7]. Particularly, as environmental concerns intensify, the demand for energy efficiency and emission reduction is increasing [8,9]. Therefore, reducing energy consumption in air separation units has become especially crucial.

The primary working principle of cryogenic air separation technology involves compressing and liquefying air, removing impurities, and cooling it. The desired products are obtained through a specialized distillation column, known as the air separation column, by utilizing the different boiling points of various components. This air separation process is essentially a specific form of cryogenic distillation. In the cryogenic air separation process, the distillation column serves as a crucial mass transfer unit operation and represents the primary energy-consuming unit. Most of the energy in the air separation process is utilized for the distillation and separation processes, with only a small portion used to produce specific gas products. Therefore, researching energy-saving solutions for this distillation process holds significant economic and societal importance [10].

Among numerous distillation energy-saving technologies, heat-integrated distillation technology has garnered extensive attention due to its significant potential for energy savings [11,12,13,14]. In heat-integrated distillation technology, the rectifying section and the stripping section are two separate columns, referred to as the rectifying column and the stripping tower, respectively [15,16,17,18,19]. These columns are connected by a compressor and a throttle valve to regulate the pressure between them. Heat exchangers are placed between the rectifying and stripping columns to achieve thermal coupling. This arrangement significantly enhances energy utilization, leading to energy savings of over 30% compared to the conventional distillation processes. The energy-saving effect of internal heat integration in the distillation process is pronounced, and since air separation is a specific form of cryogenic distillation, introducing heat integration technology into the air separation process holds vital significance. Additionally, traditional cryogenic distillation in air separation involves a high-pressure tower and a low-pressure tower, with existing single-plate coupling between the bottoms of the high-pressure column and the tops of the low-pressure column. Therefore, the practical applicability of heat integration technology in air separation is further evident [20,21,22]. Currently, some scholars have made progress in the design and optimization of HIASCs [23,24,25,26]. However, due to the unique thermal coupling structure, challenges arise in the dynamic modeling, control scheme design, and other aspects [27,28].

The dynamic modeling and corresponding control scheme design are crucial and challenging aspects of the distillation and air separation processes. The complexity increases when dealing with systems that incorporate the unique structure of internal heat coupling and the demand for high-purity and ultra-high-purity products. In recent years, wave theory has played a significant role in nonlinear dynamic modeling and control scheme design for separation processes, capturing the attention of scholars [29,30]. In the concept of wave theory for separation processes, certain physical quantities, such as the concentration and temperature distributions of the components, exhibit a constant waveform that propagates towards new steady-state positions during the dynamic separation processes. This phenomenon serves as a tool to explain and describe the nonlinear and dynamic characteristics of complex separation processes simply and accurately. For distillation and air separation processes, researchers like Luyben [31,32,33] and Hwang [34,35,36,37,38] have conducted extensive studies that demonstrated how different steady states of a distillation column can be described using different constant waveform positions. When the setpoint of the distillation column changes or external disturbances disrupt the system’s steady state, the constant waveform will propagate and move within the column until the distillation column reaches a new steady state. The waveform propagation within the column offers a fresh perspective on the nonlinear dynamic process of distillation columns, opening up new possibilities for establishing nonlinear dynamic models for distillation and air separation processes. Moreover, several control scheme performances have been enhanced to some extent based on the wave model. Some researchers have explored the integration of control schemes and wave theory, achieving certain levels of success in this endeavor [39,40,41].

However, when it comes to the HIASC, the complexity of the thermal coupling structure combined with the non-ideal nature of multicomponent air mixtures results in even more unique dynamic characteristics. Strong nonlinearities and other complex features make capturing the dynamic behavior even more challenging. Moreover, applying wave theory to the context of HIASCs requires a re-derivation of the relevant wave equations. There are relatively few reports regarding MPCs for a HIASC based on the wave model. We analyze the main reason for this is that, since many parameters need real-time estimation through live data, the wave model, although more accurate, exhibits a relatively unstable error direction. This can lead to control instability in error-compensating control schemes like MPCs. Furthermore, compared to other control schemes, MPCs exhibit less reliance on the model structure and is more suitable for solving constrained control problems.

A re-modeling study of the HIASC is conducted based on the wave theory in this work. This involves discarding the assumptions of constant molar flow and waveform velocity typically used in traditional distillation process models. This approach aims to align the model more accurately with the characteristics of the heat-coupled air separation process. Simultaneously, by leveraging the constructed wave model, a nonlinear model predictive control scheme is proposed. This scheme is compared with traditional control approaches like PID and general model control. The simulation results indicate that the proposed control scheme outperforms the other two approaches in terms of the control performance.

2. Wave Modeling of the HIASC

2.1. The Mechanism Model of the HIASC

The HIASC structure used in this study is illustrated in Figure 1. The feed compressed air (AIR) and the gas products generated from the air separation undergo heat exchange in the main heat exchanger (MH), allowing the feed stream to be pre-cooled. The compressed and cooled stream F is introduced as the feed at the bottom of the high-pressure column. Oxygen-rich liquid air (LAIR) flows through the heat exchanger (H1) and enters the top of the low-pressure column through a throttle valve (V1), providing a downward liquid phase flow. Oxygen-rich vapor air (GAIR) at the top of the low-pressure column follows the opposite path, flowing into the high-pressure column and providing an upward gas phase flow. A low-purity liquid nitrogen product (WLN) side stream is withdrawn near the argon-rich stage. A compressor (C2) and throttle valve (V1) are installed between the high-pressure and low-pressure columns, allowing the rectifying and stripping sections to operate at different pressures. The rectifying section operates at a higher pressure and temperature, allowing the LAIR stream from the high-pressure column to flow into the low-pressure column and the GAIR stream from the low-pressure column to flow into the high-pressure column. Internal heat coupling is achieved through the full coupling between the high-pressure and low-pressure columns. A certain amount of heat is transferred from each tray of the high-pressure column to the low-pressure tower, providing a downward liquid phase flow for the high-pressure column and an upward gas phase flow for the low-pressure column. This fully heat-coupled structure closely approximates the ideal heat-integrated distillation column configuration with reduced loads on the condenser and reboiler. As a consequence, the fully heat-coupled design offers superior energy efficiency. Specific operating conditions are provided in

Table 1. The operating conditions of the HIASC.

Operating Condition	Value	Operating Condition	Value
Flow rate F (kmol/s)	10	Feed composition (O2)	0.2095
Feed temperature Tf (K)	102.2	Feed composition (N2)	0.7812
Feed stage	20	Pressure of the rectifying column (MPa)	0.69
Stage number	40	Pressure of the stripping column (MPa)	0.13
Feed thermal condition	0.696	Nitrogen product purity	0.9817
Side stream stage	26	Oxygen product purity	0.9929
Side stream flow rate (kmol/s)	0.08

. The flow rate represents the flow rate of the feed. The feed stage is the location at which the feed is introduced. The stage number denotes the total number of trays in the column. The side stream stage refers to the location at which the side stream is extracted. Nitrogen and oxygen product purity represent the nitrogen and oxygen concentrations at the column’s two ends under initial steady-state conditions.

According to the principles of material conservation and energy conservation:

$$\frac{d{M}_j{x}_{i,j}}{dt}=L_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j z_{i,j}-(L_j+S{L}_j)x_{i,j}-(V_j+S{V}_j)y_{i,j};$$

(1)

$$\frac{d{M}_j{H}_j{}^L}{dt}=L_{j-1}{H}_{j-1}{}^L+V_{j+1}{H}_{j+1}{}^V+F_j{H}_j{}^F-(L_j+S{L}_j)H_j{}^L-(V_j+S{V}_j)H_j{}^V+Q_j$$

(2)

where the subscript “i” represents the component index, with “i = 1” denoting nitrogen, “i = 2” denoting argon, and “i = 3” marking oxygen. The subscript “j” indicates the stage number, where “j = 1” corresponds to the column top and “j = n” corresponds to the column bottom. The subscript “V” represents the gas phase, “L” denotes the liquid phase, and “F” represents the feed. The variables “M” represent the holdup of liquid on a tray, “L” and “V” represent the molar flow rates of gas and liquid phases, respectively, “SL” and “SV” represent the molar flow rates of liquid and gas phases for the side stream, respectively. The variable “F” stands for the feed flow rate, “H” represents the enthalpy value, “x” is the liquid phase component concentration, “y” is the gas phase component concentration, “z” is the molar fraction of the feed component, and “Q” represents the coupling heat between the rectification and stripping sections. The material conservation and energy conservation equations, combined with the gas–liquid equilibrium equation, concentration summation equation, heat coupling equation, and gas–liquid flow rate equation, constitute the mechanism model of HIASC. The mechanism model employs modeling methods from the literature [23] as a foundation for subsequent derivation of the wave model.

2.2. Wave Modeling of the HIASC

In the wave theory, we treat the concentration distribution curve as a whole to study its behavior. When operational conditions change, the concentration distribution on the trays will also change, resulting in the propagation of the concentration waves through the distillation column. As shown in Figure 2, the original Profile 1 shifts in response to the changing operational conditions, moving within the column. This movement signifies the propagation of the concentration wave within the distillation column, ultimately coming to rest at the position of Profile 2. The location of the waveform in Profile 2 represents the new steady state. The process of Profile 1 moving towards Profile 2, which is the propagation of the concentration wave, represents the dynamic response of the entire system to changes in the operational conditions.

For simplicity, we often choose a specific point on the concentration distribution curve to represent the entire curve, describing its position and movement. This point is referred to as the inflection point and is typically located where the curve has the steepest slope. Given the complexity of the components in the air separation process and the involvement of internal heat coupling, the wave modeling of the HIASC becomes more challenging. This requires us to derive the wave model specifically for the HIASC.

First, according to the material conservation equation Equation (1), we have:

$${\displaystyle \underset{j_1}{\overset{j_2}{\int }}\frac{d{M}_j{x}_{i,j}}{dt}dj=}(L_{j_1-1}{x}_{i,}{}_{j_1-1}+V_{j_2+1}{y}_{i,j_2+1}+{\displaystyle \sum _{k=j_1}^{j_2}{F}_k{z}_k})-(L_{j_2}{x}_{i,}{}_{j_2}+V_{j_1}{y}_{i,j_1}+{\displaystyle \sum _{k=j_1}^{j_2}(S{L}_{j_1}{X}_{j_1}+S{V}_{j_1}{y}_{j_1}}))$$

(3)

where $j_1$ and $j_2$ are the two ends of the integration, representing the tray positions at the two ends of the integration. Figure 3 provides a simplified schematic diagram of waveform propagation in the stripping section. Therefore, the integral on the left side of Equation (3) can be approximated as follows:

$${\displaystyle \underset{j_1}{\overset{j_2}{\int }}\frac{d{M}_j{x}_{i,j}}{dt}dj=}\frac{(M_{j_1}{x}_{i,j_1}-M_{j_2}{x}_{i,j_2})dj}{dt}$$

(4)

Combining Equations (3) and (4), the wave propagation velocity expression is obtained as follows:

$${\left.\frac{dS}{dt}\right|}_{\Delta x_{i,j}}=\frac{(L_{j_1-1}{x}_{i,}{}_{j_1-1}+V_{j_2+1}{y}_{i,j_2+1}+{\displaystyle \sum _{k=j_1}^{j_2}{F}_k{z}_k})-(L_{j_2}{x}_{i,}{}_{j_2}+V_{j_1}{y}_{i,j_1}+{\displaystyle \sum _{k=j_1}^{j_2}(S{L}_k{x}_k+S{V}_k{y}_k}))}{M_{j_1}{x}_{i,j_1}-M_{j_2}{x}_{i,j_2}}$$

(5)

where “S” represents the inflection point position of the integrated segment waveform, indicating the location of the integrated segment waveform.

Similarly, according to Equation (2), we can derive the following expression for the wave velocity:

$${\left.\frac{dS}{dt}\right|}_{\Delta H_j{}^L}=\frac{(L_{j_1-1}{H}_{j_1-1}{}^L+V_{j_2+1}{H}_{j_2+1}{}^V+{\displaystyle \sum _{k=j_1}^{j_2}{F}_k{H}_k{}^F})-(L_{j_2}{H}_{j_2}{}^L+V_{j_1}{H}_{j_1}{}^V+{\displaystyle \sum _{k=j_1}^{j_2}(S{L}_k{H}_k{}^L+S{V}_k{H}_k{}^V}))}{M_{j_1}{H}_{j_1}{}^L-M_{j_2}{H}_{j_2}{}^L}$$

(6)

By combining Equations (1) and (3), extending the integration from a local position to the entire waveform, and performing the appropriate algebraic manipulations, the formula for the wave velocity can be derived as follows:

$${\left.\frac{dS}{dt}\right|}_{x_{i,j}}=\frac{(L_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_j)-(L_j{x}_{i,j}+V_j{y}_{i,j}+S{L}_j{X}_j+S{V}_j{y}_j))}{M_{j-1}{x}_{i,j-1}-M_j{x}_{i,j}}$$

(7)

After deriving the wave velocity, we can comprehend the situation of component concentration waveforms. To eventually establish the nonlinear wave model, obtaining the actual values of the component concentrations from their waveforms is necessary. A new function is introduced to link the waveform positions with the concentration values. By observing the concentration distribution curve, the overall shape of the curve exhibits an “S” shape. In other words, there is an inflection point in the middle of the curve, and the concentrations at both ends of the curve asymptotically approach a certain straight line. Hence, empirical functions can be used for fitting. Empirical functions with substantial verification in the separation process are introduced for the waveforms of the high-pressure and low-pressure columns [42]:

$${\widehat{x}}_{i,j}=X_{i,H\_\min }+\frac{X_{i,H\_\max }-X_{i,H\_\min }}{1+e^{-{\gamma }_{i,H}(j-S_{i,H})}},j=1,2,\dots ,n/2$$

(8)

$${\widehat{x}}_{i,j}=X_{i,L\_\min }+\frac{X_{i,L\_\max }-X_{i,L\_\min }}{1+e^{-{\gamma }_{i,L}(j-S_{i,L})}},j=n/2+1,\dots ,n$$

(9)

where ${\widehat{x}}_{i,j}$ represents the estimated liquid concentration obtained through the waveform positions and empirical relationships, $X_{i,H\_\min }$ and $X_{i,H\_\max }$ denote the asymptotic limits at the extended ends of the high-pressure tower waveform, $X_{i,L\_\min }$ and $X_{i,L\_\max }$ represent the asymptotic limits at the extended ends of the low-pressure tower waveform, ${\gamma }_{i,H}$ and ${\gamma }_{i,L}$ are the slopes at the inflection points of the waveforms for the high-pressure column and low-pressure column, respectively, $S_{i,H}$ and $S_{i,L}$ are the inflection point positions of the waveforms for the two columns.

The expression for the gas–liquid equilibrium equation is as follows:

$$y_{i,j}=k_{i,j}{x}_{i,j}$$

(10)

where k is the gas–liquid equilibrium coefficient.

The heat coupling equation describes the heat transfer relationship between the rectifying section and the stripping section, and the expression is as follows:

$$Q_j=U_{ov}A\Delta T_j$$

(11)

where $U_{ov}$ is the heat transfer coefficient, $A$ is the heat transfer area, and $\Delta T_j$ is the temperature difference between a pair of heat-coupled trays.

The calculation equation for the gas–liquid equilibrium coefficient is an implicit expression of the tray temperature, pressure, and equilibrium molar fraction, with the specific relationship obtainable through the Peng-Robinson state equation [43]. The expression for the gas–liquid equilibrium coefficient is as follows:

$$k_{i,j}=k_{i,j}(x_{i,j},y_{i,j},P_j,T_j)$$

(12)

Due to the involvement of internal heat coupling, the molar flow rates of gas and liquid are no longer assumed to be constant but are considered to vary with tray position. The detailed flow rate relationships are as follows:

$$L_j=\frac{Q_j}{\lambda }+L_{j-1}+F_j{q}_j-U_j$$

(13)

$$V_j=\frac{Q_j}{\lambda }+V_{j+1}+F_j(1-q_j)-G_j$$

(14)

where $\lambda$ is the latent heat of vaporization and $q$ represents the feed thermal condition. The liquid flow of each tray comprises several components: liquid generated by condensation resulting from heat release, liquid flowing down from the tray above, the liquid portion from the feed, and the liquid portion removed in the side draw. The gas flow exhibits a similar composition. Due to heat exchange, either vaporization or condensation, the gas and liquid phase flows on each tray cannot be considered constant.

The waveform description function characterizes the system’s state at a specific moment, while the wave velocity function describes the system’s dynamic process. Combined with the gas–liquid equilibrium equation, heat coupling equation, and gas–liquid molar flow rate equation, these components can describe both the steady-state and dynamic processes of the heat-integrated air separation process. The above relations establish the nonlinear wave model for the HIASC.

2.3. Model Test

A characteristic nonlinear dynamic behavior in distillation processes is asymmetry. It means that the system transitions rapidly from one steady state to another after a disturbance occurs. However, when reverting the disturbance, the system’s return to the original equilibrium state is much slower than the speed of reaching the new equilibrium state. We can illustrate and explain this nonlinear feature using the wave model. Figure 4 displays the dynamic asymmetry of the HIASC through changes in wave velocity. At the initial moment, with a 10% variation in feed thermal condition, the wave speed of the rectifying and stripping sections reaches zero after 3000 s, and the system reaches a new steady state. At around 8000 s, when the feed thermal condition is restored to its initial value, the wave velocity decays to zero again, and the system returns to its initial steady state. Throughout this process, the dynamic behavior exhibits clear asymmetry. Taking the stripping section as an example, the wave velocity of liquid nitrogen and liquid oxygen in the stripping section experiences a significant positive step disturbance in feed thermal condition, with a high average wave speed and rapid decay speed. The stripping section reaches a new steady state within 2000 s. However, after restoring the initial feed thermal condition, both waveforms in the stripping section have smaller wave velocities and slow decay rates. The stripping section returns to the initial steady state after 5000 s, a transition process much longer than the reverse process.

The asymmetry mentioned above is a qualitative analysis of the nonlinear dynamic behavior of the HIASC using the wave model. Other dynamic characteristics of the HIASC can also be matched using the wave model, such as reverse responses, sensitivity to disturbances, and more, which will not be delved into detailed descriptions here due to the space limit. Next, let us analyze the quantitative accuracy of the wave model.

Figure 5 and Figure 6, respectively, depict the comparison between the mechanism wave models in predicting the product concentrations under a step disturbance of +5% from the feed mole fraction. The wave model can accurately track the variations in the concentrations at the top and bottom of the column. The wave model closely matches the mechanism model, whether during the initial transient period following the disturbance or in the steady-state phase after reaching a new equilibrium. During the steady-state period, a small steady-state error exists, primarily due to the model reduction. Striking a balance between accuracy and simplicity is challenging when reducing the model’s complexity. Therefore, this error should be taken into account for online estimation. However, for online control design, the precision of the wave model is sufficient. On the one hand, the reduced-order model structure facilitates model-based online control design. On the other hand, the robustness of the control algorithms can easily mitigate such minor model errors.

In summary, the wave model significantly simplifies the model structure while maintaining a high model accuracy. In scenarios where real-time performance is a critical requirement, the application of the wave model can significantly enhance computational efficiency.

3. Control Scheme Design for the HIASC

3.1. Model Predictive Control Design Based on the Wave Model

The traditional control schemes based on data-driven models or other conventional models often face issues with model mismatch, resulting in a suboptimal control performance. Nonlinear model-based predictive control (MPC) strategies that utilize dynamic optimization methods to address the rolling optimization problems can mitigate the model mismatch by extensively utilizing mechanism models. However, the mechanism models for distillation processes often involve a large number of differential equations and variables, which poses significant challenges for online optimization. The wave model proposed in this study significantly reduces the number of differential equations in the distillation process, thereby facilitating the design of an MPC. Based on the wave model, an MPC scheme for HIASC is designed, with the rectifying section pressure (Pr) and feed thermal condition (q) as control variables. The objective function is as follows:

$$\varphi (k)={\displaystyle \sum _{i=1}^N{[Z(k+i\left|k\right.)-Z_s(k+i)]}^T{\theta }_p}[Z(k+i\left|k\right.)-Z_s(k+i)]+{\displaystyle \sum _{i=1}^M\Delta u^T(k+i\left|k\right.){\theta }_c\Delta u(k+i\left|k\right.)}$$

(15)

where N is the prediction horizon. M is the control horizon. Z represents the predicted values of the controlled vector (top liquid nitrogen concentration and bottom liquid oxygen concentration). Zs is the reference trajectory for the controlled vector. Δu represents the optimized values of the control vector, composed of the feed thermal condition and the rectifying section pressure. ${\theta }_p$ and ${\theta }_c$ are diagonal weighting matrices. Z(k + i|k) represents the predicted value of the controlled variable at time k in the i-th step ahead. Δu(k + i|k) represents the optimized control vector values from time k to the i-th sampling instant in the future. The reference trajectory for the controlled variable at time k, Zs(k), is defined by an iterative approximation curve as follows:

$$Z_s(k+i\left|k\right.)=\alpha Z(k+i-1\left|k\right.)+(1-\alpha )R(k)$$

(16)

where R(k) is the final setpoint for the controlled variable. It is often assumed to be a constant to simplify the optimization process. In the HIASC, the setpoint is typically the concentration of nitrogen or oxygen, and it is generally a constant value. Therefore, it is reasonable to assume that the setpoint remains fixed. $\alpha$ is a weighting coefficient between 0 and 1. A larger value of $\alpha$ results in a smoother reference trajectory curve but may lead to reduced robustness. To compensate for small errors arising from model mismatch and other issues due to the lack of integral action in NMPC, a modified reference trajectory is used based on Equation (16):

$$d(k)=Z(k)-Z(k\left|k\right.)$$

(17)

$${\widehat{Z}}_s(k+i\left|k\right.)=Z_s(k+i\left|k\right.)-d(k)$$

(18)

where Z(k) is the measured value of the controlled variable at time k. Z(k|k) is the predicted value of the controlled variable at time k. ${\widehat{Z}}_s(k+i\left|k\right.)$ is the modified reference trajectory. In model-based control, errors in the model are inevitably present. It is commonly assumed that these errors change very slowly or remain constant over a short period. Therefore, it is essential to leverage this characteristic. At time k, we obtain a concentration measurement representing the actual value. Predictions for the concentration at time k have been made in previous time steps. The difference between these two values characterizes the error of the model. Compensating for this error in the predictions for future time steps can reduce prediction errors and, consequently, enhance control accuracy.

The boundary conditions for the control variables are set as follows:

$$0.552\le P_r\le 0.828$$

(19)

$$0.6\le q\le 0.8$$

(20)

The online rolling optimization is carried out using the control variable parameterization method. In this study, the prediction horizon is 300 s, the control horizon is 300 s, and the sample time is 15 s. The prediction horizon is typically equal to or greater than the control horizon. The reason for not using a prediction horizon longer than the control horizon here is determined by the characteristics of the wave model. Some parameters in the model need to be updated based on data collected at each sampling instant. Therefore, if the prediction horizon is too long, it can lead to an increased model inaccuracy, affecting the reliability of the prediction results. Hence, we have chosen a relatively shorter prediction horizon in this work.

3.2. Comparative Control Scheme Design

To better validate the effectiveness of the MPC control scheme, two commonly used control schemes are designed as comparative schemes for the HIASC in this study.

The control objectives for the HIASC are the nitrogen concentration at the top of the rectifying column and the oxygen concentration at the bottom of the stripping column. In the control design process, the pressure of the rectifying column (Pr) and the feed heat condition (q) are chosen as control variables. First, a relatively straightforward approach is the PID control scheme, where a dual-loop PI control strategy is implemented using these two control variables to regulate the two control objectives.

The generalized model control (GMC) scheme typically requires a linearized model of the HIASC. Assuming the steady-state polynomial model between the controlled variables and the manipulated variables is as follows:

$$X_{1,1}{}^s=a_{11}{P}_r{}^2+a_{12}{P}_r+a_{13}$$

(21)

$$X_{3,n}{}^s=a_{21}{q}^2+a_{22}q+a_{23}$$

(22)

A polynomial approximation model can be obtained by employing the least squares method to identify the coefficients in the above equations.

The control law for GMC is as follows:

$$\frac{{{dX}_{1,1}}^{*}}{dt}=k_{11}(X_{1,1}{}^{*}-X_{1,1})+k_{12}{\displaystyle \underset{0}{\overset{t}{\int }}(X_{1,1}{}^{*}-X_{1,1}})dt$$

(23)

$$\frac{{{dX}_{3,\mathrm{n}}}^{*}}{dt}=k_{21}(X_{3,n}{}^{*}-X_{3,n})+k_{22}{\displaystyle \underset{0}{\overset{t}{\int }}(X_{3,n}{}^{*}-X_{3,n}})dt$$

(24)

After taking the derivatives of Equations (21) and (22), and substituting them into Equations (23) and (24), the analytical solution for the control law is obtained as follows:

$$P_r(k)=P_r(k-1)+[k_{11}(X_{1,1}{}^{*}-X_{1,1}(k-1))+k_{12}{\displaystyle \sum _{i=1}^k(X_{1,1}{}^{*}-X_{1,1}(i))}T]/(b_{11}{P}_r(k-1)+b_{12})$$

(25)

$$q(k)=q(k-1)+[k_{11}(X_{3,n}{}^{*}-X_{3,n}(k-1))+k_{12}{\displaystyle \sum _{i=1}^k(X_{3,n}{}^{*}-X_{3,n}(i))}T]/(b_{21}q(k-1)+b_{22})$$

(26)

By using Equations (25) and (26), the values of the control variables at each time step can be iteratively solved. The values of $b_{11}$, $b_{12}$, $b_{21}$, $b_{22}$ are related to the coefficients identified in the polynomial model and can be determined through derivation. The polynomial approximation model identified in this study is as follows:

$$X_{1,1}{}^s=-0.19014P_r{}^2+0.46165P_r+0.7538$$

(27)

$$X_{3,n}{}^s=1.703q-0.19512$$

(28)

4. Comparison of Control Scheme Effects

In the forthcoming simulation tests, we compare the control effectiveness of the three control schemes (PID, GMC, MPC) from both the servo control and conventional control viewpoints to observe the variations in control performance among these schemes.

4.1. Servo Control

Figure 7 illustrates the control results of the three control strategies when the setpoint for the top nitrogen concentration is changed from 0.9817 to 0.9819. It can be observed that all three control strategies successfully achieve the control objective. However, since changing setpoints is not commonly encountered in practical processes, this simulation study only serves comparative purposes. From Figure 7, it is evident that MPC exhibits a faster control speed and higher control accuracy. In the rectifying section, MPC enables the nitrogen concentration to reach the new setpoint more rapidly. In the stripping section, despite a slightly larger overshoot, the overall convergence speed is still very rapid in the case of MPC. It might be attributed to the fact that the rectifying section approaches the new setpoint quickly, leading to more intense adjustments of the control variables and consequently resulting in a higher overshoot. Furthermore, as the system approaches near-steady-state conditions, MPC demonstrates less oscillation than PID and GMC, thus achieving a smoother response. In the rectifying section, both GMC and PID show similar control performance. However, in the stripping section, GMC outperforms PID with a smaller overshoot and faster convergence. Figure 8 shows the responses of the control variables for serve control. Apparently, for MPC, the control variables exert a significantly greater influence during the initial stages, leading to a rapid stabilization of the rectifying section but with a relatively larger overshoot in the stripping section. It essentially means that a certain degree of stability is sacrificed for the sake of speed. This trade-off is justified because when the setpoint changes, the primary concern is to bring the product concentration as close to the setpoint as quickly as possible to minimize the production of non-conforming products.

4.2. Regulatory Control

Figure 9 depicts a comparison of the control performance of the three control strategies under a 10% increase in the feed flow rate. When facing the disturbance caused by this change in feed flow rate, the advantages of the MPC approach over the other two control strategies are quite evident. These advantages include a smaller overshoot, faster response time, and better stability. For the GMC and PID control strategies, it is evident that both the rectifying and stripping sections exhibit prolonged oscillations. It indicates that the mutual coupling between these sections poses challenges for the control strategies. In the initial phase of the stripping section, the oscillations in the PID response curve are irregular, implying that the system nonlinearity has a significant impact on the PID control. On the other hand, due to its basis in nonlinear modeling, MPC can adequately capture both the nonlinear and dynamic characteristics of the system. This advantage enables MPC to handle the coupling and nonlinear issues that the other two controllers lack the capability to address. Figure 10 shows the responses of the control variables in this case. From the changes in the control variables, we can also observe that the MPC scheme exhibits a relatively rapid change in control variables at the initial time. In contrast, the other two schemes show a more delayed response, resulting in significant differences in the control performance.

Figure 11 compares the control performance of the three control strategies under a 10% increase in the side stream flow rate. GMC and PID continue to exhibit prolonged oscillations in the rectifying section, while the stability of MPC remains significantly better than the other two methods. In the stripping section, GMC and MPC perform better than PID. It further reveals the relative weakness of PID when dealing with systems characterized by strong nonlinearity. Since the amount of side stream withdrawal is notably smaller than the feed flow rate, the overall fluctuations are considerably smaller than the feed disturbance. Figure 12 shows the responses of the control variables when a 10% increase in the side stream flow rate occurs. Similar to the previous situations, MPC exhibits a more intense response of control variables at the early stage, thus resulting in a faster overall stabilization time.

Comparing Figure 9 and Figure 11, we can observe that the data in Figure 9 is an order of magnitude smaller. It suggests that the system deviation from the original equilibrium state is not as severe, and the established steady-state polynomial model might be a good fit. This improved fit potentially contributes to the better performance of GMC compared to PID. Overall, MPC consistently demonstrates excellent performance in handling both types of disturbances, indicating the effectiveness of the established wave model.

4.3. Error Index Analysis

To enable a more rigorous comparative assessment of the control efficacy exhibited by the controller configurations, we engage in a quantitative appraisal of their performance, using two key performance metrics: the absolute integral error (IAE) and the square integral error (ISE). The mathematical formulations for these quantitative indices are delineated as follows:

$$IAE=\frac{1}{n_t}{\displaystyle \sum _{k=1}^{n_t}|x(k)-x^{*}|}$$

(29)

$$ISE=\frac{1}{n_t}{\displaystyle \sum _{k=1}^{n_t}{(x(k)-x^{*})}^2}$$

(30)

where n_t is the total sampling times, x(k) is the controlled variable, and x* corresponds to the set value of the controlled variable.

Table 2. Error index of PID.

	Servo Control	F + 10%	Sf + 10%
IAE_Y1	1.98 × 10−5	3.06 × 10−5	2.93 × 10−6
IAE_Xn	3.78 × 10−6	1.51 × 10−5	9.47 × 10−6
ISE_Y1	1.76 × 10−9	3.39 × 10−9	1.89 × 10−11
ISE_Xn	3.63 × 10−11	7.68 × 10−10	3.64 × 10−10

,

Table 3. Error index of GMC.

	Servo Control	F + 10%	Sf + 10%
IAE_Y1	1.85 × 10−5	1.88 × 10−5	2.69 × 10−6
IAE_Xn	8.93 × 10−7	2.43 × 10−6	1.84 × 10−6
ISE_Y1	1.54 × 10−9	1.34 × 10−9	1.56 × 10−11
ISE_Xn	1.77 × 10−12	1.89 × 10−11	1.18 × 10−11

and

Table 4. Error index of MPC.

	Servo Control	F + 10%	Sf + 10%
IAE_Y1	1.18 × 10−5	3.00 × 10−6	6.63 × 10−7
IAE_Xn	3.16 × 10−6	2.87 × 10−6	3.26 × 10−6
ISE_Y1	9.06 × 10−10	6.18 × 10−11	1.49 × 10−12
ISE_Xn	3.84 × 10−11	4.67× 10−11	4.30 × 10−11

provide a comprehensive assessment of the performance metrics for the three control schemes, PID, GMC, and MPC, under three distinct scenarios. It is apparent that GMC and MPC consistently exhibit significantly lower error metrics when contrasted with PID. In the comparison between MPC and GMC, MPC demonstrates a more balanced performance in the two-end control.

5. Conclusions

In this study, a reduced-order HIASC model based on the nonlinear wave theory is proposed to address the challenges faced by the traditional mechanism models in certain online scenarios. The traditional mechanism models involve material balance equations for each tray, resulting in a significant number of differential equations and variables. This complexity makes online solving difficult, especially for model-based control schemes where direct application of the mechanism models can be challenging.

The nonlinear wave model simplifies the model description by treating the concentration profile as a whole and studying its movement and distribution. In this wave model, only two differential equations for the rectifying and stripping sections are needed, eliminating the need for calculations on each tray. Furthermore, due to its mechanism foundation, this model offers higher accuracy than other simplified models, which is beneficial for online model-based control.

Subsequently, we proposed a model predictive control (MPC) strategy based on this model to address the online control problem of HIASC. For comparison purposes, we also employed the PID and GMC control strategies. PID utilized the common dual-loop PI control structure, while GMC was based on a polynomial model identified through data. Through simulation experiments, MPC demonstrated superior control performance compared to the other control strategies in both servo and regulatory control, which further validates the effectiveness of the proposed wave model. Using the wave model in the MPC strategy, the control variables can exhibit a faster response in the initial stages, resulting in an overall better system response.

Note that certain assumptions were made to simplify the derivation, which could introduce model mismatch errors. However, the rolling optimization inherent to MPC helps alleviate these errors. If the model is used for other online estimation tasks, careful consideration of the impact of these errors on the results would be necessary. In the future, further work should focus on establishing more accurate models to address scenarios requiring higher precision results.