© The Author(s), published by EDP Sciences, 2025

Abbreviations

BPP: Bipolar Plates

CL: Catalytic Layer

I_st: Input Current of stack

P_sm: Supply Manifold Pressure

P_ca: Cathode Pressure

P_O2: Oxygen partial pressure

P_N2: Nitrogen partial pressure

P_H2: Hydrogen partial pressure

GDL: Gas Diffusion Layer

FC: Fuel Cell

T_st: Stack Temperature

P_net: Net Power

V_st: Stack Voltage

P_rm: Return Manifold Pressure

MEA: Membrane Electrodes and Assemblies

MPC: Model Predict Control

SNFTSMC: Simplified Non-singular Fast Terminal Sliding Mode Control

NFTSMC: Non-singular Fast Terminal Sliding Mode Control

OER: Oxygen Excess Ratio

PID: Proportional Integral Derivative

PEM: Proton Exchange Membrane

PEMFC: Proton Exchange Membrane Fuel Cell

SMC: Sliding Mode Control

SMO: Sliding Mode Observer

TSMC: Terminal Siding Mode Control

1 Introduction

The IEA’s 2023 World Energy Outlook urges energy transitions by 2030, focusing on clean energy growth, fossil fuel reduction, and efficiency improvements. Hydrogen, abundant, efficient, and clean, is crucial for mitigating the greenhouse effect and driving China’s energy transformation [1–3]. Compared with traditional engine, fuel cells have the advantages of low emissions, silent operation, no pollution and so on [4, 5]. Compared to other fuel cells, PEMFCs are most prevalent due to their mature technology. They generate electricity via an electrochemical reaction between hydrogen and atmospheric oxygen, achieving high energy conversion efficiency [6, 7]. However, the cost of PEMFC is still high compared with other fossil energy sources, which leads to some obstacles in its application [8, 9]. The content of oxygen supplied to the system has a great impact on PEMFC stack. Achieving optimal oxygen content in real-time through precise control of the inlet gas is crucial for maintaining optimal output. Furthermore, coordinating air inlet and exhaust to maintain balanced cathode and anode pressures is crucial for preventing nitrogen penetration and enhancing PEMFC stability.

In our previous paper [10, 11], we studied air distribution in fuel cell stacks with varying channel structures using a 3D electrochemical-multiphysical coupling model. Additionally, we integrated this model with a thermomechanical model to investigate the impact of component structures on stack mechanical stress distribution [12, 13]. In Pukrushpan’s thesis, he systematically modeled the PEMFC stack and its subsystems using Simulink [14]. Damour et al. proposed a self-tuning PID controller based on a neural model for PEMFC oxygen supplement. Compared traditional PID, It demonstrates superior tracking capabilities, enhanced disturbance rejection in experiment [15]. Meanwhile, Baroud et al. adopted a hybrid fuzzy PID approach to further optimize fuzzy PID controller, based on a simplified mathematic model [16]. Afsharinejad et al. and Y. Zhu et al. designed linear and nonlinear optimal controller to regulate oxygen in PEMFC respectively [17, 18], and their process of model linearization is exceptionally meticulous. The results clearly demonstrate that the optimal controller has achieved significant advancements in comparison to traditional controllers. Ziogou’s et al. experiment successfully verified the operation of the MPC controller in PEMFC supply model [19]. Y. Wang et al. combined the MPC with a PID controller further, and the simulation results indicate that the performance is improved compared with MPC [20]. Napole et al. and Abbaker et al. adopted optimized terminal sliding mode control based on reaching law and delay estimation algorithm respectively for PEMFC air-feeding system, improving performance compared to terminal sliding mode [21, 22]. Apart from classical and modern control, intelligent control is also a current research hotspot, encompassing areas such as reinforcement learning [23, 24], neural networks [25], advanced fuzzy algorithm [26] and so on.

Multi-objective control displays vital importance in energy system [27–30]. Currently, there is a paucity of research on multi-objective coordinated control for PEMFC, with the majority of studies primarily focusing on oxygen excess ratio (OER) control. Coordinating the air inlet and exhaust to ensure consistent cathode and anode pressures is essential for balancing the internal stack pressure and preventing nitrogen penetration. Tang et al. employed PID and sliding mode controllers for the coordinated control of the anode supply [31]. X. Zhang et al. designed an adaptive robust predictive controller aimed at PEMFC coordinated management [32]. J. Li et al. and Song et al. adopted reinforcement learning approaches to construct a multi-objective controller for PEMFC provisioning [33, 34]. However, as presented in Table 1, the performance attained by the controllers proposed in these articles is relatively limited. Both PID and predictive control necessitate the linearization of the PEMFC system, which significantly elevates the complexity of controller design. Concurrently, intelligent controllers exhibit limited performance enhancements in comparison to modern control achieved and possess stringent performance requirements for actuators.

Based on these, dual Non-singular Fast Terminal Sliding Mode Controllers (NFTSMC) are employed to optimize oxygen content and pressure control under dynamic loading conditions. The NFTSM structure ensures finite-time convergence while mitigating potential singularity issues associated with Terminal Sliding Mode (TSM) [35].

The structure of this study is organized as follows: Initially, a mathematical model of the cathode supply system for the PEMFC is established. Subsequently, the design and detailed analysis of a dual no NFTSM controller are presented. Following this, an evaluation of the stability of the control system designed for cathode air management in PEMFCs is conducted. Furthermore, simulation experiments employing the proposed controller are designed under two distinct disturbances scenarios. Lastly, the results are rigorously validated, analyzed, and compared.

2 PEMFC model

The PEMFC stack primarily consists of several subsystems, including the supply system, reactor system, thermal management system and cooling system. As shown in Figure 1, these subsystems collaborate to accomplish energy conversion and control tasks. In this paper, it is ascertained that sufficient compressed hydrogen, optimal humidification, and appropriate temperature.

Figure 1 Sketch diagram of typical PEMFC stack.

2.1 Typical PEMFC output electromotive force

The output voltage of a single PEMFC can be mathematically denoted as,$$V_{\mathrm{Cell}}=E_{\mathrm{Nern}}-V_{\mathrm{Act}}-V_{\mathrm{Ohm}}-V_{\mathrm{Con}},$$(1)where E_Nern is theoretical maximum output voltage of a single fuel cell; V_Act represents the activation overpotential; V_Ohm denotes Ohmic polarization loss voltage, V_Con is the concentration polarization loss voltage. These parameters can be calculated as [36],$$\{\begin{array}{c}{E}_{\mathrm{n}}=1.229-a\left(T_{\mathrm{st}}-298.15\right)+b{T}_{\mathrm{st}}\left[\ln \left(P_{{\mathrm{H}}_2}\right)+\frac{1}{2}\ln \left(P_{{\mathrm{O}}_2}\right)\right]\\ V_{\mathrm{act}}=-\left[{\epsilon }_1+{\epsilon }_2{T}_{\mathrm{st}}+{\epsilon }_3{T}_{\mathrm{st}}\ln \left(C_{O_2}\right)+{\epsilon }_4{T}_{\mathrm{st}}\ln I_{\mathrm{st}}\right]\\ \begin{array}{c}{V}_{\mathrm{ohm}}=\frac{I_{\mathrm{st}}}{1000}A\left(R_C+R_M\right)\\ V_{\mathrm{con}}=c\exp \left(\mathrm{d}{I}_{\mathrm{st}}\right)\end{array}\end{array},$$(2)where a and b are constant and other parameters are summarized in Table 2. c is a temperature-dependent model parameter, and $C_{O_2}$ represents the dissolved concentration of oxygen at the gas–liquid interface, which can be determined using Henry’s law. They can be calculated by$$\{\begin{array}{c}c=\{\begin{array}{c}1.1\times 10^{-4}-1.2\times 10^{-6}\left(T_{\mathrm{st}}-273.15\right),\hspace{1em}\mathrm{if} T_{\mathrm{st}}>312.15\hspace{0.5em}\mathrm{K}\\ 3.3\times 10^{-3}-8.2\times 10^{-5}\left(T_{\mathrm{st}}-273.15\right),\hspace{1em}\mathrm{if} T_{\mathrm{st}}<312.15\hspace{0.5em}\mathrm{K}\end{array}\\ C_{{\mathrm{O}}_2}=\frac{P_{{\mathrm{O}}_2}}{5.08\times 10^6\exp \left(-\frac{498}{T_{\mathrm{st}}}\right)}\end{array}.$$(3)

Considering the restricted power rating of an individual cell, PEMFC stack generally incorporates multiple single cells connected in series, constituting a group aimed at augmenting the output voltage. The output voltage can be expressed as [26],$$V_{\mathrm{st}}=n{V}_{\mathrm{cell}}.$$(4)

2.2 Actuator models

2.2.1 Air compressor model

The primary function of the air compressor is to provide adequate air pressure and flow to the fuel cell stack, thereby ensuring its normal operation [14]. The mathematical modeling formula is presented as,$$\{\begin{array}{c}\frac{\mathrm{d}{n}_{\mathrm{motor}}}{\mathrm{d}t}=\frac{{\eta }_{\mathrm{cm}}{k}_t{k}_v}{J_{\mathrm{cp}}{R}_{\mathrm{cm}}}{n}_{\mathrm{motor}}-\frac{C_P{T}_{\mathrm{atm}}}{J_{\mathrm{cp}}{\eta }_{\mathrm{cp}}{n}_{\mathrm{motor}}}\left({\left(\frac{M_{a,\mathrm{atm}}{V}_{\mathrm{sm}}{P}_{\mathrm{sm}}}{R{T}_{\mathrm{atm}}\gamma }\right)}^{\frac{\gamma -1}{\gamma }}-1\right){\omega }_{\mathrm{cp}}+\frac{{\eta }_{\mathrm{cm}}{k}_t}{J_{\mathrm{cp}}{R}_{\mathrm{cm}}}{u}_{\mathrm{cp}}\\ {\omega }_{\mathrm{cp}}=\frac{{\omega }_{\mathrm{cp}}^{\max }{n}_{\mathrm{motor}}}{n_{\mathrm{motor}}^{\max }}\left(1-\exp \left(\frac{-r(s+\frac{n_{\mathrm{motor}}^2}{q}-P_{\mathrm{sm}})}{s+\frac{n_{\mathrm{motor}}^2}{q}-P_{\mathrm{sm}}^{\min }}\right)\right),\end{array}$$(5)where n_motor is the rotational speed of motor, u_cp represents the compressor input voltage, ω_cp denotes the airflow rate of compressor, P_sm indicates the pressure of supply manifold, and motor and compressor relevant parameters are summarized in Tables 3 and 4 [16].

2.2.2 Solenoid valve

Solenoid valve located at the cathode side outlet serves as an actuator to regulate and maintain cathode pressure, with the control signal input determining the valve opening. The mathematical modeling formula is presented as [20],$$\{\begin{array}{c}\frac{\mathrm{d}{P}_{\mathrm{rm}}}{\mathrm{d}t}=\frac{R_{\mathrm{air}}{T}_{\mathrm{rm}}}{V_{\mathrm{rm}}}\left(W_{\mathrm{ca},\mathrm{out}}-W_{\mathrm{rm}}\right)\\ W_{\mathrm{rm}}=\{\begin{array}{c}\frac{C_{\mathrm{D}}{A}_{\mathrm{T}}{u}_{\mathrm{\theta }}{P}_{\mathrm{rm}}}{\sqrt{{\mathrm{RT}}_{\mathrm{rm}}\left(\gamma -1\right)}}\sqrt{2\gamma \left({\left(\frac{P_{\mathrm{atm}}}{P_{\mathrm{rm}}}\right)}^{\frac{2}{\gamma }}-{\left(\frac{P_{\mathrm{atm}}}{P_{\mathrm{rm}}}\right)}^{\frac{1+\gamma }{\gamma }}\right)}\hspace{1em}\mathrm{if} \frac{P_{\mathrm{atm}}}{P_{\mathrm{rm}}}>0.528\\ \frac{C_{\mathrm{D}}{A}_{\mathrm{T}}{u}_{\mathrm{\theta }}{P}_{\mathrm{rm}}}{\sqrt{{\mathrm{RT}}_{\mathrm{rm}}}}{\left(\gamma {\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{\gamma -1}}\right)}^{\frac{1}{2}}\hspace{1em}\mathrm{if} \frac{P_{\mathrm{atm}}}{P_{\mathrm{rm}}}<0.528\end{array}\end{array},$$(6)where P_rm is the pressure of exhaust manifold, W_ca,out and W_rm represent the cathode outlet and exhaust airflow respectively, u_θ denotes the valve opening, and other parameters are collected in Table 4.

2.3 Establishment of system space model

The controller design is based on a fifth-order model of the PEMFC cathode supply system. Both the system state quantities and supply system structure are illustrated in Figure 2. In Figure 2, W_ca,in and W_ca,out are inlet and exhaust airflow of the cathode respectively [37]. P_sm, P_ca, and P_rm denote the pressures in the supply manifold, PEMFC cathode, and return manifold, respectively. u_cp and u_θ are two signal outputs of the control system: u_cp is used to control the voltage of the air compressor, and u_θ regulates the opening of the solenoid valve. This nonlinear spatial state equations can be described as [32, 34],$$\{\begin{array}{c}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}=b_1(x_4-x_1-x_2-b_2)-\frac{b_3}{b_4}(x_1+x_2-x_5+b_2)\frac{x_1}{x_1+x_2+b_2}\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=b_8(x_4-x_1-x_2-b_2)-\frac{b_3}{b_5}(x_1+x_2-x_5+b_2)\frac{x_2}{x_1+x_2+b_2}\\ \begin{array}{c}\frac{\mathrm{d}{x}_3}{\mathrm{d}t}=-b_9{x}_3-\frac{b_{10}}{x_3}\left({\left(\frac{x_4}{b_{14}}\right)}^{b_{12}}-1\right)\phi \left(x\right)+b_{13}{u}_{\mathrm{cp}}\\ \frac{\mathrm{d}{x}_4}{\mathrm{d}t}=b_{14}\left(1+b_{15}\left({\left(\frac{x_4}{b_{14}}\right)}^{b_{12}}-1\right)\right)\times (\phi \left(x\right)-b_{16}-x_1-x_2-b_2)\\ \frac{\mathrm{d}{x}_5}{\mathrm{d}t}=b_{25}(x_1+x_2-x_5+b_2)-b_{20}{x}_5{u}_{\theta }\end{array}\end{array}$$(7)where x_i ∈ R⁵ (i = 1, …, 5) represents the system state quantities (x₁ is the cathode oxygen pressure, x₂ represents the cathode nitrogen pressure, x₃ = n_motor, x₄ = P_sm, x₅ = P_rm), the control input vector is u = (u_cp, u_θ) ∈ R², φ(x) = ω_cp, and b₁ to b₂₅ are constant coefficients which are collected in Appendix.

Figure 2 Supply system structure.

2.4 Control objectives

An insufficient supply of oxygen not only diminishes the efficiency of the electrochemical reaction but also precipitates extremely high concentrations that facilitate nitrogen penetration, ultimately posing a risk of irreversible damage to the PEMFC stack.

The Oxygen Excess Ratio (OER) serves as a crucial indicator for evaluating the oxygen content within PEMFC. It represents the ratio of the oxygen flow to the stoichiometric oxygen requirement necessary for complete fuel oxidation [16].$${\lambda }_{{\mathrm{O}}_2}=\frac{Q_{m,\mathrm{air},\mathrm{provide}}}{Q_{m,\mathrm{air}}}=\frac{4{k^2}_{\mathrm{ca},\mathrm{in}}{y}_{{\mathrm{O}}_2,\mathrm{atm}}F(P_{\mathrm{sm}}-P_{O_2}-P_{N_2}-P_{\mathrm{sat}})}{n{M}_{\mathrm{a},\mathrm{atm}}(1+{\omega }_{\mathrm{atm}})}$$(8)where Q_{m,air,provide} is the incoming oxygen flow and Q_m,air is the consumed oxygen flow (Fig. 3).

Figure 3 The diagram of power variation with OER under selected currents.

By analyzing the change curve between the net power and OER under various current conditions in, the optimal OER for each current level is identified. Fitting these data points of D. Yang’s research in Figure 3, the best optimal OER can be described as [38],$${\lambda }_{\mathrm{opt}{\mathrm{O}}_2}=a_1{I}_{\mathrm{st}}^2+a_2{I}_{\mathrm{st}}+a_3{I}_{\mathrm{st}}^{-1}+a_4,$$(9)where a_x is the correlation coefficient (a₁ = 1.368e–5, a₂ = −9.784e–3, a₃ = −56.84, a₄ = 3.873).

Meanwhile, the other control value is the cathode pressure P_ca and it is described as [32],$$P_{\mathrm{ca}}=P_{O_2}+P_{N_2}+P_{\mathrm{sat}}.$$(10)

Based on these, the affine form for PEMFC cathode supply system is presented as,$$\{\begin{array}{c}\frac{\mathrm{d}{x}_i}{\mathrm{d}t}=f(x_i,D)+\mathrm{gu}\\ y_i=h(x_i,D)\end{array},$$(11)where D = I_st ∈ R is regarded as a measurable disturbance input vector, $h=({\lambda }_{{\mathrm{O}}_2},P_{\mathrm{ca}}{)}^T\in R^2$ is the vector of control objectives.

3 Hybrid controller design

In this section, a complex composite robust control strategy is proposed to achieve effective synergistic control of the cathode, considering various perturbations. Two Non-singular Fast Terminal Sliding Mode (NFTSM) structure controllers are combined as control system to optimize oxygen content and pressure control under dynamic loading conditions. The non-singular terminal sliding mode structure ensures finite-time convergence while mitigating potential singularity issues associated with traditional terminal sliding modes. The control block diagram is illustrated in Figure 4.

Figure 4 The control block diagram of the sliding mode.

In Terminal Sliding Mode (TSM) control, when the system state approaches the equilibrium point, certain terms in the control law may lead to a zero denominator, subsequently resulting in the singularity problem. By redefining the sliding mode surface and designing the corresponding control law, the NFTSM control effectively avoids the singularity problem while ensuring system convergence within a finite time.

According to above equations, the error vector can be denoted as [35],$$\mathit{e}=(e_1,e_2{)}^T=(y_{1\mathrm{opt}}-y_1,y_{2\mathrm{ref}}-y_2{)}^T\in R^2,$$(12)where $e_1=\frac{b_{23}}{b_{24}{I}_{\mathrm{st}}}(x_4-x_1-x_2-b_2)-y_{1\mathrm{opt}}$ is the error signal of compressor, and e₂ = (x₁ + x₂ + b₂) − y_2ref is the error signal of solenoid valve.

Then, the derivative of the error vector can be described as,$$\frac{\mathrm{d}\mathit{e}}{\mathrm{d}t}={\left(\frac{\mathrm{d}{e}_1}{\mathrm{d}\mathrm{t}},\frac{\mathrm{d}{e}_2}{\mathrm{d}\mathrm{t}}\right)}^T={\left[\frac{b_{23}}{b_{24}{I}_{\mathrm{st}}}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}\mathrm{t}}-\frac{\mathrm{d}{x}_1}{\mathrm{d}\mathrm{t}}-\frac{\mathrm{d}{x}_2}{\mathrm{d}\mathrm{t}}\right),\frac{\mathrm{d}{x}_4}{\mathrm{d}\mathrm{t}}-\frac{\mathrm{d}{x}_1}{\mathrm{d}\mathrm{t}}-\frac{\mathrm{d}{x}_2}{\mathrm{d}\mathrm{t}}\right]}^T.$$(13)

Further,$$\frac{{\mathrm{d}}^{2}e}{\mathrm{d}{t}^2}={\left(\frac{{\mathrm{d}}^2{e}_1}{\mathrm{d}{t}^2},\frac{{\mathrm{d}}^2{e}_2}{\mathrm{d}{t}^2}\right)}^T={\left[\frac{b_{23}}{b_{24}{I}_{\mathrm{st}}}\left(\frac{{\mathrm{d}}^2{x}_4}{\mathrm{d}{t}^2}-\frac{{\mathrm{d}}^2{x}_1}{\mathrm{d}{t}^2}-\frac{{\mathrm{d}}^2{x}_2}{\mathrm{d}{t}^2}\right),\left(\frac{{\mathrm{d}}^2{x}_1}{\mathrm{d}{t}^2}+\frac{{\mathrm{d}}^2{x}_2}{\mathrm{d}{t}^2}\right)\right]}^T,$$(14)where above second order derivatives of state quantities are represented as,$$\{\begin{array}{c}\frac{{\mathrm{d}}^2{x}_1}{\mathrm{d}{t}^2}=b_1\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}-\frac{\mathrm{d}{x}_1}{\mathrm{d}t}-\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)+\frac{b_3}{b_4}\left(\frac{\mathrm{d}{x}_5}{\mathrm{d}t}\left(\frac{x_5(x_1+x_2)\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}-\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)+x_1(x_1+x_2+b_5)}{(x_1+x_2+b_2{)}^2}-\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)\right)\\ \frac{{\mathrm{d}}^2{x}_2}{\mathrm{d}{t}^2}=b_8\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}-\frac{\mathrm{d}{x}_1}{\mathrm{d}t}-\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)+\frac{b_3}{b_5}\left(\left(\frac{x_5\left((x_1+b_2)\frac{\mathrm{d}{x}_2}{\mathrm{d}t}-x_2\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)+x_2(x_1+x_2+b_5)\frac{\mathrm{d}{x}_5}{\mathrm{d}t}}{(x_1+x_2+b_2{)}^2}\right)-\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)\\ \frac{{\mathrm{d}}^2{x}_4}{\mathrm{d}{t}^2}=b_{14}{b}_{16}\left(1+b_{15}\left({\left(\frac{x_4}{b_{11}}\right)}^{b_{12}}-1\right)\right)\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}+\frac{\mathrm{d}{x}_2}{\mathrm{d}t}+\frac{\mathrm{d}{x}_3}{\mathrm{d}t}\right)+\frac{\partial }{\partial x_4}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\end{array}$$(15)where $\frac{\partial }{\partial x_4}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)=b_{12}{b}_{14}{b}_{15}{b_{11}}^{-b_{12}}\left(\phi \right(x)-q_{16}(x_4-x_1-x_2-b_2\left)\right)+b_{14}(1+b_{15}\left({\left(\frac{x_4}{b_{11}}\right)}^{b_{12}}-1)\right)\left(\frac{\partial \phi \left(x\right)}{\partial x_4}-b_{16}\right)$.

3.1 Compressor control

Based on Section 2.2, compressor can be described as [16],$$\phi \left(x\right)=\frac{{\phi }_{\max }}{x_3^{\max }}\left(x_3-\exp \left(\frac{-r\left(s+\frac{x_3^2}{q}-x_4\right)}{s+\frac{x_3^2}{q}-x_{{ }_4}^{\min }}\right)\right).$$(16)

The sliding variable of compressor control is designed as,$$s_1=\frac{\mathrm{d}{e}_1}{\mathrm{d}t}+{\kappa }_1{e}_1+{\kappa }_2\vartheta \left(e_1\right),$$(17)where κ₁ and κ₂ are the NFTSM gains, and ϑ(e₁) is,$$\vartheta \left(e_1\right)=\{\begin{array}{c}{e}_1^{c1/c2}\hspace{1em}\mathrm{if} s_1^{\mathrm{mark}}=0\hspace{1em}\mathrm{or}\hspace{0.5em}{s}_1^{\mathrm{mark}}\ne 0, \left|e_1\right|\ge \frac{1}{\infty }\\ d_1{e}_1+d_2\mathrm{sign}\left(e_1\right)e_1^2\hspace{1em}\mathrm{if}\hspace{0.5em}{s}_1^{\mathrm{mark}}\ne 0, \left|e_1\right|<\frac{1}{\infty }\end{array},$$(18)where $s_1^{\mathrm{mark}}=\frac{\mathrm{d}{e}_1}{\mathrm{d}t}+{\kappa }_1{e}_1+{\kappa }_2{e_1}^{c_1/c_2}$, c₁ and c₂ are the positive odd numbers, and they satisfy 1/2 < c₁/c₂ < 1. d₁ and d₂ are selected as $d_1=(2-c_1/c_2)(\frac{1}{\infty }{)}^{{ }^{c_1/c_2-1}}$ and $d_2=(c_1/c_2-1)(\frac{1}{\infty }{)}^{{ }^{c_1/c_2-2}}$.

Let $\frac{\mathrm{d}{s}_1}{\mathrm{d}t}=0$, u_cmp can be inversely calculated as,$$u_{\mathrm{cmp}}=\frac{1}{q\left(x\right)}\left(-\sigma \mathrm{sign}\left(s_1\right)-i{s}_1-p\left(x\right)\right),$$(19)where p(x) and q(x) can be represented as,$$\{\begin{array}{c}\begin{array}{c}p\left(x\right)={\kappa }_1{b}_{16}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}-\frac{\mathrm{d}{x}_1}{\mathrm{d}t}-\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)-b_{16}\left(\frac{{\mathrm{d}}^2{x}_1}{\mathrm{d}{t}^2}+\frac{{\mathrm{d}}^2{x}_2}{\mathrm{d}{t}^2}\right)+b_{16}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\frac{\partial }{\partial x_4}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)+\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\frac{\partial }{\partial x_1}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)+\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\frac{\partial }{\partial x_2}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)\right)\\ +b_{16}\left(-b_{{ }_9}{x}_3-\frac{b_{10}}{x_3}\left({\left(\frac{x_4}{b_{14}}\right)}^{b_{12}}-1\right)\right)\phi \left(x\right)\frac{\partial }{\partial x_3}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)+{\kappa }_2\theta \left(e_1\right)\end{array}\\ q\left(x\right)=b_{13}{b}_{16}\frac{\partial }{\partial x_3}\left(\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\right)\end{array}.$$(20)

3.2 Solenoid valve control

To minimize computation and complexity, a Simplified NFSTM (SNFSTM) controller is employed in this section. The sliding variable of solenoid valve control is designed as,$$s_2=\frac{\mathrm{d}{e}_2}{\mathrm{d}t}+\frac{1}{\lambda }{e}_2^{\frac{j}{k}},$$(21)where λ > 0 is the sliding mode gain of SNFSTM, j and k are the positive odd numbers (j > k).

Same as the air compressor, the control law of the solenoid valve can be expressed as,$$u_{\vartheta }=\frac{1}{n\left(x\right)}\left(-\sigma \mathrm{sign}\left(s_2\right)-i{s}_2-m\left(x\right)\right),$$(22)where m(x) and n(x) can be represented as,$$\{\begin{array}{c}m\left(x\right)=b_{25}\left(x_1+x_2+c_2-x_5\right)(\frac{\delta }{\delta x_5}\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)+\left(\frac{\delta }{\delta x_1}\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)+\frac{\delta }{\delta x_1}\left(\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)\right)\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\\ +\left(\frac{\delta }{\delta x_2}\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)+\frac{\delta }{\delta x_2}\left(\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)\right)\frac{\mathrm{d}{x}_2}{\mathrm{d}t}+\left(\frac{\delta }{\delta x_4}\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)+\frac{\delta }{\delta x_4}\left(\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)\right)\frac{\mathrm{d}{x}_4}{\mathrm{d}t}\\ \begin{array}{c}+\frac{\delta }{\delta x_5}\left(\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right))+\frac{j}{\lambda k}\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}+\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)(x_1+x_2{)}^{\frac{j}{k}-1}\\ n\left(x\right)=b_{20}\left(\frac{\delta }{\delta x_5}\left(\frac{\mathrm{d}{x}_1}{\mathrm{d}t}\right)+\frac{\delta }{\delta x_5}\left(\frac{\mathrm{d}{x}_2}{\mathrm{d}t}\right)\right)\end{array}\end{array}.$$(23)

3.3 Control stability verification

Firstly, the Lyapunov function is designed as,$$v=\frac{1}{2}{s}_1^2+\frac{1}{2}{s}_2^2.$$(24)

Because v > 0 satisfies Lyapunov stability, further the derivative of v can be described as,$$\begin{array}{c}\mathrm{d}v=s_1\frac{\mathrm{d}{s}_1}{\mathrm{d}t}+s_2\frac{\mathrm{d}{s}_2}{\mathrm{d}t}\\ =s_1\left(p\right(x)+q(x\left)u_{\mathrm{cm}p}\right)+s_2\left(m\right(x)+m(x\left)u_{\vartheta }\right)\\ \begin{array}{c}=s_1\left(p\left(x\right)+q\left(x\right)\left(\frac{1}{q\left(x\right)}(-\sigma \mathrm{sign}(s_1)-i{s}_1-p(x)\right)\right)+s_2\left(m\left(x\right)+n\left(x\right)\frac{1}{n\left(x\right)}(-\sigma \mathrm{sign}(s_2)-i{s}_2-m(x\left)\right)\right)\\ =s_1(-\sigma \mathrm{sign}(s_1)-i{s}_1)+s_2(-\sigma \mathrm{sign}(s_2)-i{s}_2)\end{array}\end{array}.$$(25)

It’s obvious that dv ≤ 0 is always valid, so it fulfills the criteria of the Lyapunov stability theorem, thereby ensuring the stability of the system.

4 Result and discussion

In this research paper, two distinct types of variable loads are employed to validate the performance of the control system. Case 1 is the step current load and Case 2 is the step current load with noise disturbance. In Case 2, a white noise perturbation is introduced into the system’s signal. This case is specifically designed to evaluate the system’s capability to reject noise perturbations and its effectiveness in tracking the optimal trajectory under conditions that closely mimic the actual operating environment.

Case 1: Step current load

In Case 1, the current is varied at intervals of 20 s, and this disturbance is intended to simulate fluctuations in the external load of stack. As depicted in Figure 5a, the current exhibits a gradual increase from 0 to 40 s. Then, at the 60 s, it initially undergoes a decrease and subsequently experiences an increase. As shown in Figures 5b and 5c, the trend observed in the output voltage and power changes aligns with the current, demonstrating that the control system effectively maintains stable and efficient output of PEMFC stack. With the sudden change of input current, the output exhibits a response time of approximately 0.6 s. This delay can be attributed to the PEMFC electrochemical reaction.

Figure 5 Step current load and system outputs in step current load.

The variations in the system’s state parameters further confirm the satisfactory operation of the PEMFC stack as depicted in Figure 6. Owing to the controllers’ design as NTSM structures, the air compressor motor speed, gas supply pressure, and gas exhaust pressure all attain finite-time convergence with minimal delay. Meanwhile, the resultant overshoot by controller remains within acceptable boundaries. A slight delay in the nitrogen and oxygen pressures can also be attributed to the electrochemical reactions occurring within the PEMFC stack.

Figure 6 System state parameters.

Figures 7a and 7b illustrate the fluctuations in the system error and the OER signal, respectively. The red line represents the proposed NFTSMC, the blue line corresponds to TSMC, and the black line denotes SMC. Although NFTSMC and TSMC achieve finite time convergence, TSMC is significantly downperformed by NFTSMC in terms of overshoot and tracking pefformance from subplots. As for SMC, due to its traditional sliding mode structure, it has the lower performance. Relevant control performance indexes of Figure 7 are summarized in Table 5. It is evident that the NFTSMC exhibits superior tracking performance and robustness.

Figure 7 Performance comparisons in step current load.

Meanwhile, the cathode pressure error curve and cathode pressure curve are denoted in Figures 7c and 7d, respectively. The red dotted line is SNFTSMC and the blue is TSMC. It is evident that both controllers exhibit comparable tracking performance. Nonetheless, it is observable that the TSMC possesses a steady-state error which is approximately 100 Pa higher than that of the SNFTSMC during steady-state conditions. Although this discrepancy in error is deemed insignificant for the operation of the stack, the simplified NTSM structure controller is employed in this study due to its superior performance and the potential reduction in singularities of TSMC.

Case 2: Step current load with noise disturbance

A disturbance of white noise is incorporated into the system’s signal in this case, with the purpose of assessing the control system’s further potentials. The trend of current changes is consistent with Case 1, as shown in Figure 8a. It’s obvious that the controller exhibits a moderate level of suppression against current noise disturbances as depicted in Figures 8b and 8c.

Figure 8 Step current load and system outputs in step current load with disturbance.

NTSMC demonstrated superior control performance, particularly in managing noise within the locally amplified regions of the OER error and OER tracking curves. Compared to the traditional SMC and the TSMC, the NTSMC outperforms tracking accuracy, which is same as Case 1.

Notably, in the subgraph of Figure 9, NTSMC shows enhanced anti-interference capabilities, as it exhibits less fluctuation under noise interference compared to both TSMC and SMC, which display a more pronounced amplitude of change. Meanwhile, in the subgraph of Figure 10, the TSMC and SMC struggle to stabilize their errors at zero due to their slower convergence speed when subjected to continuous noise disturbances, whereas NTSMC still ensures quick convergence, driving the error towards zero.

Figure 9 OER error curve in step current load with disturbance.

Figure 10 OER curve in step current load with disturbance.

Based on the pressure error curve depicted in Figure 11, it is evident that the steady-state error in the step current is insignificant as a result of the inclusion of noise. Nevertheless, an examination of both the error and tracking curves reveals that the TSMC exhibits notably poorer suppression of noise perturbations compared to SNFTSMC. Furthermore, TSMC demonstrates a more pronounced overshoot when confronted with abrupt current changes. The data with disturbances is organized in the subsequent Table 6.

Figure 11 Cathode pressure and error curves in step current load with disturbance.

5 Conclusion

In this paper, a complex robust controller is employed to manage oxygen content and cathode pressure of PEMFC system. The primary conclusions derived are as follows,

• (i)

  The proposed control system is designed with a dual NFTSM controller for coordinated control, ensuring finite-time convergence of the system while effectively addressing potential singularity issues associated with TSM.

• (ii)

  The results demonstrate that NFTSMC attains a significant twofold improvement in tracking performance compared to TSMC, while also augmenting system stability. Additionally, NFTSMC exhibits superior performance in noise suppression when juxtaposed with other structural methodologies.

• (iii)

  The proposed controller successfully regulates the intake air with precision to achieve the optimal oxygen level in real-time, thereby ensuring the system consistently operates at peak output. Furthermore, it maintains uniform cathode and anode pressures, which is crucial for equilibrating internal stack pressures and mitigating nitrogen permeation.

• (iv)

  The proposed control method exhibits robustness and high performance in effectively managing load variations, external disturbances, and uncertainties associated with PEMFC operation, all of which are vital for both the operational efficiency and longevity of PEMFC.

Acknowledgments

We gratefully acknowledge the financial support from the National Natural Science Foundation of China (22179054), Ministry of Science and Technology of the People’s Republic of China (G2022014065L), and the Innovation Support Program of Science and Technology Program of Jiangsu Province (SBZ2023080107).

References

Appendix

The following are the internal and auxiliary system parameters, denoted as b_i where i ranges from 1 to 25, for the PEMFC system found in state equations.

$b_1=\frac{R{T}_{\mathrm{st}}{k}_{\mathrm{ca},\mathrm{in}}}{M_{\mathrm{O}2}{V}_{\mathrm{ca}}}\left(\frac{{\chi }_{{\mathrm{O}}_2,\mathrm{atm}}}{1+{\omega }_{\mathrm{atm}}}\right)$, $b_2=P_{\mathrm{sat}}$, $b_3=\frac{R{T}_{\mathrm{st}}{k}_{\mathrm{ca},\mathrm{out}}}{V_{\mathrm{ca}}}$, $b_4=M_{O_2}$, $b_5=M_{N_2}$, $b_6=M_V{P}_{\mathrm{sat}}$, $b_7=\frac{R{T}_{\mathrm{st}}n}{4F{V}_{\mathrm{ca}}}$, $b_8=\frac{R{T}_{\mathrm{st}}{k}_{\mathrm{ca},\mathrm{in}}}{M_{N2}{V}_{\mathrm{ca}}}\left(\frac{1-{\chi }_{O_2,\mathrm{atm}}}{1+{\omega }_{\mathrm{atm}}}\right)$, $b_9=\frac{k_t{k}_v{\eta }_{\mathrm{cm}}}{R_{\mathrm{cm}}{J}_{\mathrm{cp}}}$, $b_{10}=\frac{T_{\mathrm{atm}}{C}_p}{J_{\mathrm{cp}}{\eta }_{\mathrm{cp}}}$, $b_{11}=P_{\mathrm{atm}}$, $b_{12}=\frac{\gamma -1}{\gamma }$, $b_{13}=\frac{k_t{\eta }_{\mathrm{cm}}}{R_{\mathrm{cm}}{J}_{\mathrm{cp}}}$, $b_{14}=\frac{R{T}_{\mathrm{atm}}\gamma }{M_{\mathrm{atm}}{V}_{\mathrm{sm}}}$, $b_{15}=\frac{1}{{\eta }_{\mathrm{cp}}}$, $b_{16}=k_{\mathrm{ca},\mathrm{in}}$, $b_{17}=\frac{A_T{C}_D}{\sqrt{R{T}_{\mathrm{st}}}}\sqrt{\frac{2\gamma }{\gamma -1}}$, $b_{18}=\frac{1}{\gamma }$, $b_{19}=(\frac{2}{\gamma +1}{)}^{\frac{\gamma }{\gamma -1}}$, $b_{20}=\frac{R{T}_{\mathrm{st}}{A}_T{C}_D}{V_{\mathrm{rm}}\sqrt{R{T}_{\mathrm{st}}}}{\gamma }^{0.5}(\frac{2}{\gamma +1}{)}^{\frac{\gamma +1}{2\gamma -1}}$, $b_{21}=\frac{1}{R_{\mathrm{cm}}}$, $b_{22}=k_v$, $b_{23}=\frac{X_{O_2,\mathrm{atm}}}{1+{\omega }_{\mathrm{atm}}}$, $b_{24}=\frac{n{M}_{O_2}}{4F}$, $b_{25}=\frac{R{T}_{\mathrm{st}}{k}_{\mathrm{ca},\mathrm{out}}}{V_{\mathrm{rm}}}$.

All Tables

All Figures

	Figure 1 Sketch diagram of typical PEMFC stack.
In the text

	Figure 2 Supply system structure.
In the text

	Figure 3 The diagram of power variation with OER under selected currents.
In the text

	Figure 4 The control block diagram of the sliding mode.
In the text

	Figure 5 Step current load and system outputs in step current load.
In the text

	Figure 6 System state parameters.
In the text

	Figure 7 Performance comparisons in step current load.
In the text

	Figure 8 Step current load and system outputs in step current load with disturbance.
In the text

	Figure 9 OER error curve in step current load with disturbance.
In the text

	Figure 10 OER curve in step current load with disturbance.
In the text

	Figure 11 Cathode pressure and error curves in step current load with disturbance.
In the text