Robust observer-based-α-variable model-free supertwisting fractional order sliding mode control for nonlinear PEMFC system with uncertainties and disturbance

Abstract

The service life and efficiency of proton exchange membrane fuel cells (PEMFCs) are significantly related to the control performance of the air supply system. Therefore, this research develops a novel robust observer-based- $\alpha$ -variable model-free supertwisting fractional-order sliding mode control ( $\alpha (t)$ -MF-STFOSMC) for complex nonlinear PEMFC air supply system with unmeasurable state variables, model uncertainties, and external disturbance. First, the fifth-order nonlinear PEMFC air supply system model is presented, and the unmeasured state variables are estimated using a nonlinear disturbance observer (NDOB1). Second, the ultra-local model (ULM) algorithm is employed to reconstruct the complex nonlinear PEMFC air supply system, avoiding the need for precise modeling and reducing the controller design difficulty. To estimate and compensate for the uncertain dynamics in the ULM algorithm, another NDOB2 is proposed to realize zero estimation error and finite-time observation. Besides, to achieve optimal control performance with less input chattering, finite-time convergence, and fast response speed, the $\alpha (t)$ -MF-STFOSMC is designed using super-twisting fractional-order sliding mode control (STFOSMC) algorithm. Furthermore, the stability of $\alpha (t)$ -MF-STFOSMC via a closed-loop system is verified using the Lyapunov theorem. Finally, the nonlinear PEMFC air supply system model with the proposed controller is realized in MATLAB/Simulink environment, and the simulation results are given to show the superiority and effectiveness of the proposed technique.

Introduction

With the growing environmental challenges created by traditional fuels, fuel cells have recently received much attention. A proton exchange membrane fuel cell (PEMFC) is an emerging renewable energy source that show great potential in reducing emissions. Due to its promising characteristics, it is considered for further development.^1–3 Moreover, the PEMFC plant is capable of being one of the most popular innovative technologies replacing conventional internal combustion engines. It is ideal for ground vehicles because of its advantages, such as fast startup, low operating temperature, low operating noise, small volume, and zero carbon dioxide emission.^4,5

An automotive PEMFC plant mainly includes an air (oxidant) supply subsystem, hydrogen (fuel) supply subsystem, humidification/thermal management subsystem, fuel cell stack, and balance of plant (BOP) control unit.^6–8 To achieve optimal performance and safe operation, it is critical to keep good cooperation among subsystems via the control system of the BOP to provide suitable conditions such as pressure, flow rate, temperature, and humidity for the electrochemical reaction inside the PEMFC system. However, the reaction performance, output efficiency, and the PEMFC stack’s lifetime are significantly influenced by the air supply subsystem, which provides the necessary compressed air to the cathode channel of the PEMFC stack under operating conditions and time-varying load current. ⁹ In particular, when the current load changes abruptly during the PEMFC system operation, the compressed air in the cathode channel of the PEMFC stack will be rapidly consumed to fulfill the power requirements. On the one hand, if the air compressor is unable to react immediately to rapid current load changes via the control unit to provide enough oxygen, the stack cathode will suffer from the well-known phenomenon of so-called oxygen starvation, which leads to undesirable voltage drops and may cause physical and thermal harm to the PEMFC system. ¹⁰ On the other hand, a larger amount of oxygen may result in “oxygen saturation,” which will inevitably increase the air compressor’s parasitic power, even up to 30% of the PEMFC stack power.^8,11 Furthermore, the performance of the PEMFC stack will also deteriorate due to membrane dehydration and increased membrane impedance. ¹² It is well-known that the oxygen excess ratio (OER) assesses the sufficiency of air feeding. Therefore, the OER must be maintained within/to its desired range/value via effective control under various operating conditions.

In the recent two decades, various control schemes have been widely proposed for the PEMFC air-feeding system. The most significant framework in the early stage is common model-based control methods. Pukrushpan et al. ¹³ proposed a linear quadratic regulator (LGR) for a 75 kW PEMFC system model to control the desired OER value. a proportional-integral-derivative (PID) controller is presented based on a linearized model established at an obtained equilibrium point. ¹⁴ In addition, a feedback linearization method-based-optimal control theory framework is proposed by Chen et al. ¹⁵ Similarly, a feedforward combined with an LQR feedback control approach is developed.^16,17 Unfortunately, in practical applications, the precise identification of the structure and parameters of the PEMFC system mode is often challenging due to the time-varying disturbances, strong nonlinearity, and high uncertainties. Consequently, methods that heavily rely on high-accuracy plant models are unlikely to ensure the most efficient global regulation of OER. Considering the nonlinearities and dynamic characteristics of the air supply subsystem, nonlinear model predictive control is developed, ¹⁸ which can effectively address the nonlinearities and operating constraints; however, the fast solution of optimization issues adds some extra challenges in practice. Moreover, some researchers switched to intelligent modeling-based-control approaches to avoid the high dependence on precise system models. For example, fuzzy logic systems^19,20 and radial basis function neural networks^21,22 have been developed to online estimate the unmolded or uncertain dynamics of the PEMFC air supply system, laying the foundation for the effective control of OER regulation. Moreover, fuzzy rules are formulated to automatically tune the coefficients of the PID controller to improve the adaptability of intricate operating conditions. ⁴ Nevertheless, there is a lack of dependable guidelines for offline tuning of fuzzy logic systems/rules, as well as for determining the structure, learning parameters, and selections of neural networks. Consequently, these tasks often rely on the experience and expertise of the engineering developer. It is well known that sliding mode control (SMC) is a robust nonlinear control approach, which has very low sensitivity against uncertainties and disturbances of the controlled system while providing fast response speed and high control precision with a simple structure. Different structures of SMC have been extensively employed in addressing the control problem of OER; for example, a simple SMC of a fifth-order nonlinear PEMFC model ²³ and cascaded adaptive SMC^8,24 are developed. Although SMC can address nonlinear characteristics and uncertainties, chattering on control variables is usually generated near the equilibrium point under a high feedback gain, leading to system performance deterioration. In this regard, Rakhtala et al. ²⁵ and Pilloni et al. ²⁶ developed a high-order SMC based on a super-twisting algorithm to regulate OER at its pre-defined value. The findings exhibit that the designed controller performs a good effect on the desired OER tracking and disturbance rejection immunity. Similarly, a new variable gain second-order SMC using a super-twisting algorithm is designed to control the air-breathing system of the PEMFC system. ²⁷ The simulation results show that the presented method has the advantages of the fixed gain in terms of robustness against parametric uncertainties, disturbance with the unknown boundary, and chattering reduction. Still, this strategy did not cover all operating conditions of the PEMFC system, such as the cathode pressure, which is an unmeasurable variable in practice. Consequently, how to calculate the OER and design an effective controller to achieve given steady-state and transient responses without intense chattering has hardly been reported. Furthermore, most of the methods mentioned above rely heavily on highly accurate models and are unlikely to ensure the most efficient possible control of OER.

On the other hand, the total pressure inside the cathode channel of the PEMFC stack, which is called cathode pressure (CP), cannot be accurately measured by any advanced sensor, and CP will be used to calculate the OER value; therefore, in practice, it is not a directly available variable. Nevertheless, the most relevant references neglected this fact and assumed that both CP and OER are measurable variables. In order to address such a research problem, various state observers have been proposed to estimate/reconstruct the unmeasured CP and OER, such as Kalman Filters^13,28 and Luenberger observer. ²⁹ However, due to high nonlinearity, model uncertainties in fuel cell systems, and requirements of the availability of some dynamic state variables make these observers intricate and sometimes undesirable. Algebraic observers based on numerical differentiation have been proposed,^30,31 but these observers are sensitive to measurement noise and can not ensure globally consistent estimation accuracy. Moreover, a high-order sliding mode observer has been proposed for the PEMFC system,^25,26 but the solution of the observability matrix requires many calculations due to the complexity and high nonlinearities of PEMFC system as well as how to estimate the second and third order of the output estimation errors is added extra challenge in practical application. Although the existing research provided different degrees of success in estimating unmeasurable state variables and OER control, there are still some limitations and shortcomings.

Fortunately, model-free control (MFC), which only relies on knowledge of input and output information. Therefore, it has been developed gradually to avoid the requirement for obtaining model acknowledgment. Traditional MFCs are widely applied in different fields and offer satisfactory performance with appropriate coefficients. ³² Furthermore, other advanced model-free schemes have received significant attention, such as intelligent PID (iPID),^33–36 data-driven model-free adaptive control, ³⁷ iterative learning control, ³⁸ and active disturbance rejection control. ³⁹

Among the different MFCs, iPID based on an ultra-local model has been intensively studied in recent years.^1,33–36,40 In order to improve the control performance and ensure the stability of the controlled system, the iPID has been developed and integrated with other methods; for instance, SMC, ¹ fractional-order SMC, ⁴¹ adaptive fuzzy logic control, ⁴² and neural network control. ⁴³ However, the robustness of these methods depends on the designed observer/estimator to approximate the unknown lumped uncertainties, such as algebraic observer (AO) ⁴⁰ and time-delay estimation.^42,43 However, due to time delays and time windows, respectively, both the TDE and AO inevitably have approximation errors. Additionally, MFC based on extended state observers (ESOs) is developed.^1,41,44–46 The ESO can only achieve asymptotic observation when the time derivative of the unknown uncertain dynamics reaches to zero, the estimation error will converge to zero as time tends to infinity. Therefore, the zero estimation error and finite-time observation of the uncertain system dynamics are not considered in the aforementioned approaches. Thus, the NDOB is the most effective one ⁴⁷ to estimate the uncertain or unmodeled system dynamics, thanks to its advantages such as simplicity in design structure, a few coefficients, and low sensitivity against disturbances and measurement noise immunity.

Considering the valuable contributions introduced by the studies mentioned above, however, it is notable that most existing methods regarding the OER regulation in the PEMFC system focused on reducing the effect of external disturbance and parameter uncertainties, ignoring the observation of cathode pressure. Thus, this research develops a novel $\alpha (t)$ -MF-STFOSMC using NDOB for complex nonlinear PEMFC air supply system with unmeasurable state variables, model uncertainties, and external disturbance. On the one hand, the designed NDOB estimates unmeasurable state variables and uncertain system dynamics to achieve fast convergence and zero estimation error. On the other hand, designing the fractional-order SMC using a supertwisting approach offer optimal control performance and achieve finite-time convergence. In addition, the supertwisting algorithm is introduced due to its advantages of reducing the chattering phenomenon and accelerating the reaching speed.

The main contribution and novelties of the proposed $\alpha (t)$ -MF-STFOSMC approach are summarized as follows:

1. This study represents the pioneering effort to address the both problems of estimation and control of OER. Different from existing methods,^{21,22,25,26,30,31} a new NDOB1 is designed for OER reconstruction, aiming to provide high estimation precision and robustness against parameter variations and noise immunity.

2. To avoid the need for a precise model and simplify the controller design process, the ULM algorithm is utilized to reconstruct the complex nonlinear PEMFC system thanks to MFC theory, wherein the second NDOB2 is proposed to observe the unknown, uncertain system dynamics.

3. On the basis of the ULM algorithm, the designed NDOB2 and STFOSMC are utilized to calculate the control law of the proposed $\alpha (t)$ -MF-STFOSMC method, aiming to provide optimal control performance with fast response speed, less chattering, and finite-time convergence.

4. In the existing ULM algorithm-based MFCs,^{1,33,34,40,42,43} the value of $\alpha$ is fixed and pre-defined. In He et al., ³⁵ and Wang et al. ³⁶ an adaptive techniques using a least-squares algorithm and new nonlinear function are presented to tune $\alpha$ automatically. Unlike,^35,36 a novel $\alpha$ -variable law is proposed to adjust the value of control gain automatically, aiming to enhance the tracking accuracy of the MFC.

5. The stability of the proposed $\alpha (t)$ -MF-STFOSMC approach is analyzed using the Lyapunov approach, and the nonlinear PEMFC air supply system model with the proposed method is realized in the MATLAB/Simulink environment. Moreover, unmeasurable variables, external disturbances, and uncertainties are considered in the simulation to verify the robustness and efficiency of the proposed technique.

This reminder of this paper can be organized as follows: Problem formulation and dynamic modeling of air supply system are presented in the next section. Then, the design structure of the proposed $\alpha (t)$ -MF-STFOSMC approach and its stability analysis are introduced. Next, the simulation results and discussion are presented. Finally, the conclusion of this paper is elaborated.

Preliminaries and problem formulation

Preliminaries

The definition of the fractional-order can be represented by the general operator structures, ${{}_{t_0}D}_t^{\overline{\gamma }}$ and ${{}_{t_0}I}_t^{\overline{\gamma }}$ , which represent a generalization of the differential and integral operators, respectively.^48,49

Definition 1. The derivative of the function $z(t)$ with fractional order $\overline{\gamma }$ based on Riemann-Liouville definition can be given as:

$${{}_{t_0}D}_t^{\overline{\gamma }}z(t)=\frac{d^{\overline{\gamma }}z(t)}{d{t}^{\overline{\gamma }}}=\frac{1}{\mathrm{\Gamma }(n-\overline{\gamma })}\frac{d^n}{d{t}^n}{\int }_{t_0}^t(t-s{)}^{n-\overline{\gamma }-1}z(s)\mathrm{ds}$$

(1)

where $t_0$ is the initial time, and $n$ is the first integer larger than $\overline{\gamma }$ , that is, $n-1\le \overline{\gamma }<n$ .

The integration of function $z(t)$ with fractional order $\overline{\gamma }$ based on Riemann-Liouville definition can be expressed as:

$${{}_{t_0}I}_t^{\overline{\gamma }}z(t)=\frac{1}{\mathrm{\Gamma }(\overline{\gamma })}{\int }_{t_0}^t(t-s{)}^{\overline{\gamma }-1}z(s)\mathrm{ds}$$

(2)

where $\overline{\gamma }\in {\mathbf{R}}^{+}$ and $\mathrm{\Gamma }(·)$ is Euler’s Gamma function.

Property 1. If $n$ is an integer, there exists

$$\frac{d^n}{d{t}^n}[{{}_{t_0}D}_t^{\overline{\gamma }}z(t)]{=}_{t_0}{D}_t^{n+\overline{\gamma }}z(t)$$

(3)

Thus, the following notations can be utilized for convenience: (i) ${{}_{t_0}D}_t^{\overline{\gamma }}z(t)=D^{\overline{\gamma }}z(t)$ (ii) $\frac{dz(t)}{dt}=\stackrel{·}{z}(t)$ .

Lemma 1. ⁵⁰ Young’s inequality: Let $a>0$ , $b>0$ , $r>1$ , $g>1$ , $\overline{\epsilon }>1$ , and $\frac{1}{r}+\frac{1}{g}=1$ , then the following inequality holds:

$$\mathrm{ab}\le \frac{\overline{\epsilon }{a}^r}{r}+\frac{{\overline{\epsilon }}^{-r/g}{b}^g}{g}$$

(4)

In particular, when $r=g=2$ , the above (4) becomes

$$\mathrm{ab}\le \frac{\overline{\epsilon }}{2}{a}^2+\frac{1}{2\overline{\epsilon }}{b}^2$$

(5)

Lemma 2. ⁵¹ Consider the Lyapunov function $V(t)$ with initial value $V(0)$ such that the following inequality holds,

$$\stackrel{·}{V}(t)\le -m{V}^d(t),\forall t\ge 0,V(t_0)\ge 0$$

(6)

where $0<d<1$ and $m>0$ . Then, the finite time $t_f$ can be calculated as in the following equation:

$$t_{\mathrm{f}}\le \frac{V^{1-d}(t_0)}{m(1-d)}$$

(7)

Problem formulation

The structure of the integrated PEMFC plant, illustrated in Figure 1, is mainly composed of the hydrogen supply sub-system, the air supply sub-system, and the thermal and humidification management sub-system. ⁶ A high-pressure hydrogen tank stores and supplies hydrogen for the hydrogen supply subsystem, and an electronic valve can rapidly adjust the hydrogen flow into the anode channel of the PEMFC stack to match the airflow. Therefore, hydrogen flow is considered to be sufficient and under control to keep track of the cathode’s air pressure. Furthermore, because the behavior of thermodynamics is significantly slower than the behavior of aerodynamics, the humidity and temperature of the PEMFC system can be regulated independently by local control. In particular, it is assumed that the humidity and temperature of the PEMFC system can operate at the required level because the humidity and temperature values cannot increase or decrease fast. Thus, the main task is concentrated on the air supply system, and it is critical to construct a robust and efficient controller for the air supply system to keep the mass flow rate and partial pressure of oxygen at desired levels on the cathode.

Figure 1.

The main subsystems of the PEMFC system.

Dynamic modeling of PEMFC air supply system

In practice, auxiliary components such as an air compressor, an air supply manifold, a cooler, and a humidifier assist the air supply system operate reliably. It is reasonable to assume that other subsystems’ constraint requirements are most likely well-controlled. According to References 5 and 23, the model of the reduced fifth-order nonlinear PEMFC system is presented to describe the dynamic behaviors of airflow, and its state-space representation can be defined as follows:

$$\stackrel{·}{X}=F(X)+G_{U}U+G_{D}D$$

(8)

The state vector $X=[X_1,X_2,X_3,X_4,X_5{]}^T=[P_{O_2},P_{N_2},$ ${P_{v,\mathrm{ca}},{\omega }_{\mathrm{cp}},P_{\mathrm{sm}}]}^T$ respectively denote the oxygen partial pressure, the nitrogen partial pressure, vapor partial pressure, the angular speed of the compressor motor, and the supply manifold pressure. $U=v_{\mathrm{cm}}$ is the control input (the voltage of the air compressor motor), and $D=I_{\mathrm{st}}$ is the external disturbance (load current). The vector of nonlinear functions $F(X)=[F_1(X),F_2(X),F_3(X),F_4(X),F_5(X){]}^{\mathrm{T}}$ is given as follows

$$\{\begin{array}{c}{F}_1(X)=b_1(X_4-\chi )-\frac{b_2{X}_1\varphi (\chi )}{b_3{X}_1+b_4{X}_2+b_5{X}_3}\\ F_2(X)=b_7(X_4-\chi )-\frac{b_2{X}_2\varphi (\chi )}{b_3{X}_1+b_4{X}_2+b_5{X}_3}\\ F_3(X)=b_8(X_4-\chi )-\frac{b_2{X}_2\varphi (\chi )}{b_3{X}_1+b_4{X}_2+b_5{X}_3}+b_{10}\\ F_4(X)=-b_{11}{X}_4-\frac{b_{12}}{X_4}\left({\left(\frac{X_5}{b_{16}}\right)}^{b_{14}}-1\right)\times W_{cp}(X_4,X_5)\\ F_5(X)=b_{16}\left(1+b_{17}\left({\left(\frac{X_5}{b_{13}}\right)}^{b_{14}}-1\right)\right)\times (W_{cp}(X_4,X_5)-b_{18}(X_5-\chi ))\end{array}$$

(9)

where $\chi =X_1+X_2+X_3$ stands for the cathode pressure, $\varphi (\chi )$ represents the total flow rate at the outlet of cathode channel and $W_{\mathrm{cp}}(X_4,X_5)$ denotes the air flow rate inside the air compressor.

$$\varphi (\chi )=\{\begin{array}{ccc}{b}_{19}\chi {\left(\frac{b_{13}}{\chi }\right)}^{b_{20}}\sqrt{1-{\left(\frac{b_{13}}{\chi }\right)}^{b_{14}}},\hfill & \mathrm{if}\hfill & \frac{b_{13}}{\chi }>b_{21}\hfill \\ \multicolumn{1}{c}{b_{22}\chi ,}& \mathrm{if}\hfill & \frac{b_{13}}{\chi }\le b_{21}\hfill \end{array}$$

(10)

$$W_{\mathrm{cp}}(X_4,X_5)=\frac{{W_{\mathrm{cp}}}^{\max }{X}_4}{x_4^{\max }}(1-e^{\left(\frac{-c(\kappa +\frac{X_4^2}{q}-X_5)}{\kappa +\frac{X_4^2}{q}-X_5^{\min }}\right)})$$

(11)

and

$$\begin{array}{c}{G}_U=[000b_{14}{0]}^{\mathrm{T}},\hfill \\ \multicolumn{1}{c}{G_D=[-b_{6}02b_6{00]}^{\mathrm{T}}.}\end{array}$$

(12)

The constants $b_i$ , for $\{i=1,...,25\}$ are functions of the physical parameters of the PEMFC system, which are listed in Table 1.^13,23,52

Table 1.

Constants of the PEMFC air supply system.

$b_1=\frac{R{T}_{\mathrm{st}}{k}_{\mathrm{ca},\mathrm{in}}}{M_{O_2}{V}_{\mathrm{ca}}}\left(\frac{x_{O_2,\mathrm{atm}}}{1+w_{\mathrm{atm}}}\right)$	$b_{14}=\frac{{\eta }_{\mathrm{cm}}{k}_t}{J_{\mathrm{cp}}{R}_{\mathrm{cm}}}$
$b_2=\frac{R{T}_{\mathrm{st}}}{V_{\mathrm{ca}}}$	$b_{15}=\frac{R{T}_{\mathrm{atm}}\gamma }{M_{a,\mathrm{atm}}{V}_{\mathrm{sm}}}$
$b_3=M_{O_2}$	$b_{16}=\frac{1}{{\eta }_{\mathrm{cp}}}$
$b_4=M_{N_2}$	$b_{17}=k_{\mathrm{ca},\mathrm{in}}$
$b_5=M_{v,\mathrm{ca}}$	$b_{18}=\frac{C_D{A}_T}{\sqrt{R{T}_{\mathrm{st}}}}\sqrt{\frac{2\gamma }{\gamma -1}}$
$b_6=\frac{R{T}_{\mathrm{st}}N}{4F{V}_{\mathrm{ca}}}$	$b_{19}=\frac{1}{\gamma }$
$b_7=\frac{R{T}_{\mathrm{st}}{k}_{\mathrm{ca},\mathrm{in}}}{M_{N_2}{V}_{\mathrm{ca}}}\left(\frac{1-x_{O_2,\mathrm{atm}}}{1+{\omega }_{\mathrm{atm}}}\right)$	$b_{20}={\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma }{\gamma -1}}$
$b_8=k_{\mathrm{ca},\mathrm{in}}\frac{R_v{T}_{\mathrm{st}}{\omega }_{\mathrm{atm}}}{1+{\omega }_{\mathrm{atm}}}$	$b_{21}=\frac{C_D{A}_T}{\sqrt{R{T}_{\mathrm{st}}}}\sqrt{\gamma }{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{2\gamma -1}}$
$b_9=N_{v,\mathrm{membr}}{A}_{\mathrm{FC}}N\frac{R{T}_{\mathrm{st}}}{V_{\mathrm{ca}}}$	$b_{22}=k_{\mathrm{ca},\mathrm{in}}\left(\frac{x_{O_2,\mathrm{atm}}}{1+w_{\mathrm{atm}}}\right)$
$b_{10}=\frac{{\eta }_{\mathrm{cm}}{k}_t{k}_v}{J_{\mathrm{cp}}{R}_{\mathrm{cm}}}$	$b_{23}=\frac{N{M}_{O_2}}{4F}$
$b_{11}=\frac{C_p{T}_{\mathrm{atm}}}{J_{\mathrm{cp}}{\eta }_{\mathrm{cp}}}$	$b_{24}=\frac{1}{R_{\mathrm{cm}}}$
$b_{12}=P_{\mathrm{atm}}$	$b_{25}=k_v$
$b_{13}=\frac{\gamma -1}{\gamma }$	$x_{O_2,\mathrm{atm}}=\frac{y_{O_2,\mathrm{atm}}{M}_{O_2}}{M_{a,\mathrm{atm}}}$
	${\omega }_{\mathrm{atm}}=\frac{M_v}{M_{a,\mathrm{atm}}}\frac{{\varphi }_{\mathrm{atm}}{p}_{\mathrm{sat}}}{p_{\mathrm{atm}}-{\varphi }_{\mathrm{atm}}{p}_{\mathrm{sat}}}$

Control objectives

Oxygen excess ratio (OER) ${\lambda }_{O_2}$ is considered a key variable that impacts the operating performance of the PEMFC system, which is defined as the ratio of the inlet flow of oxygen to oxygen reaction. It can be expressed as follows:

$${\lambda }_{O_2}=\frac{W_{O_2\mathrm{in}}}{W_{O_2\mathrm{react}}}=\frac{b_{22}(X_5-\chi )}{b_{23}D}$$

(13)

where $W_{O_2\mathrm{in}}$ is the oxygen inlet flow, $W_{O_2\mathrm{react}}$ is the oxygen reaction consumption flow.

OER is a crucial performance variable that indicates the current status of airflow in the PEMFC system. Oxygen starvation in the PEMFC system is typically caused by low OER, whereas high OER can avoid oxygen starvation and enhance the PEMFC system’s efficiency. Unfortunately, high OER will cause the high-speed operation of the air compressor and lead to a reduction in the net power. Consequently, to prevent oxygen starvation and maximize net power, the OER value must be adjusted to its optimal point.

The output net power $P_{\mathrm{net}}$ can be defined as the difference between the power output of the PEMFC stack and the consumed power by the air compressor.

$$P_{\mathrm{net}}=P_{\mathrm{st}}-P_{\mathrm{cp}}$$

(14)

where $P_{\mathrm{st}}$ is the total output power generated by the PEMFC stack and $P_{\mathrm{cp}}$ is the power consumed by the air compressor.

$$P_{\mathrm{st}}=V_{\mathrm{st}}{I}_{\mathrm{st}}=N{V}_{\mathrm{fc}}{I}_{\mathrm{st}}$$

(15)

$$P_{\mathrm{cp}}=b_{24}U(U-b_{25}{X}_4)$$

(16)

where $N$ is the number of cells and $V_{\mathrm{fc}}$ the single cell’s voltage; more details can be found in Pukrushpan et al. ¹³

According to equations (13)–(16), it can be concluded that the OER value and the net power are significantly influenced by the magnitude of the load current. The relationship between the output net power and OER value ${\lambda }_{O_2}$ under the different stack currents is displayed in Figure 2(a). It can be observed that from Figure 2(a), the maximum value of output net power $P_{\mathrm{net}}$ is obtained when the OER ${\lambda }_{O_2}$ within the range between 2 and 2.5 under different stack currents from 100 to 300 A. According to References 13, 14, and 52, the ${\lambda }_{O_2}=2$ is the most studied at present, which has been considered as the desired value. The other opinion is to take the reference value of OER ${\lambda }_{O_2}$ as a variable and keep OER track of its optimal value ${\lambda }_{O_2}^{\mathrm{opt}}$ to obtain maximum output net power. Based on Pukrushpan et al., ¹³ a lookup table can be generated using OER, which maximizes the net power, as illustrated in Figure 2(b). By forcing the actual OER to track its optimal OER, the maximum output net power can be obtained. Hence, the main control goal for the air supply system is to adjust the OER value to avoid oxygen starvation and achieve maximum power point tracking under various load conditions. The ${\lambda }_{{\mathrm{O}}_2}^{\mathrm{Opt}}$ is determined based on the stack current $I_{\mathrm{st}}$ as follows:

$${\lambda }_{O_2}^{\mathrm{Opt}}=c_3\times I_{\mathrm{st}}^3+c_2\times I_{\mathrm{st}}^2+c_1\times I_{\mathrm{st}}+c_0$$

(17)

where the coefficients $c_0=3.3263554$ , $c_1=-1.211507\times {10}^{-2}$ , $c_2=4.33996\times {10}^{-5}$ , and $c_3=-5.7967\times {10}^{-8}$ are obtained using the least squares method.

Figure 2.

(a) Curves for the relationship between output net power and OER at different load currents and (b) fitting curve for the optimal OER and load currents.

Control design

In view of (8), the lack of the cathode pressure, which is difficult to quantify, is the greatest barrier to determining the oxygen excess ratio (OER) ${\lambda }_{O_2}$ . Moreover, sensors cannot always be used for measurements because of the high costs associated with the detecting technology or because the quantity is not always easily quantifiable, mainly when humidified gas streams are available inside a fuel cell stack. State observers are of tremendous importance in these situations because they can estimate the unavailable values in place of physical sensors.^30,31 Thus, we propose an NDOB1 that provides a high level of precision and rapid convergence for reconstructing the OER value ${\lambda }_{O_2}$ . Then, the estimated value is employed in the design of $\alpha$ -variable model-free super-twisting fractional-order sliding mode control ( $\alpha (t)$ -MF-STFOSMC; see Figure 3). The proposed control method is composed of two parts: the first part is NDOB1 for estimating the unmeasurable state variables, and the second part is $\alpha (t)$ -MF-STFOSMC for calculating the control law, and their details will be given in the following subsections.

Figure 3.

Schematic diagram of proposed $\alpha (t)$ -MF-STFOSMC controller.

Nonlinear disturbance observer-1 (NDOB1) for estimation of CP and OER

From the view of (8), one has

$${\stackrel{·}{X}}_5=h(X_5)(W_{\mathrm{cp}}(X_4,X_5)-b_{17}{X}_5)+\psi (\chi ,X_5)$$

(18)

where $h(X_5)=b_{15}(1+b_{16}({\left(\frac{X_5}{b_{12}}\right)}^{b_{13}}-1))$ and $\psi (\chi ,X_5)=b_{17}h(X_5)\chi$ is considered as unknown function due to including unmeasurable variables $\chi$ .

Remark 1. Given that $X_5$ is a physical variable and its corresponding parameters $b_{12},b_{13},b_{15},b_{16}$ , and $b_{17}$ are all positive, the part ${\left(\frac{X_5}{b_{12}}\right)}^{b_{13}}-1>0$ is ensured in the operation region, it means that the function $h(X_5)$ is also positive. Therefore, it is reasonable to consider that the reciprocal value of $h(X_5)$ exists and is bounded in the operation range $X_5\in {\mathcal{X}}_5$ .

Assumption 1. There exists a constant $\eta >0$ such that the derivative $\stackrel{·}{\psi }(\chi ,X_5)$ of the unknown uncertain function $\psi (\chi ,X_5)$ satisfies $|\stackrel{·}{\psi }(\chi ,X_5)|\le \eta$ .

Note that the Assumption 1 is a standard condition widely used for the convergence analysis of NDOB1.

Let’s assume that the estimated value of the unknown function $\psi (t)$ , and its estimation error can be defined as $\hat{\psi }(t)$ and $\tilde{\psi }(t)=\psi (t)-\hat{\psi }(t)$ , respectively. We can denote $\mathrm{\Phi }(t)=\psi (t)-\epsilon X_5$ as a new vertical variable, where $\epsilon >0$ is a positive constant.

Next, when differentiating $\mathrm{\Phi }(t)$ , we obtain

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathrm{\Phi }}(t)=\stackrel{·}{\psi }(t)-\epsilon {\stackrel{·}{X}}_5(t)}\\ \multicolumn{1}{c}{=\stackrel{·}{\psi }(t)-\epsilon [h(X_5)(W_{\mathrm{cp}}(X_4,X_5)-b_{17}{X}_5)+\psi (t)]}\\ \multicolumn{1}{c}{=\stackrel{·}{\psi }(t)-\epsilon [\epsilon X_5+h(X_5)(W_{\mathrm{cp}}(X_4,X_5)-b_{17}{X}_5)]}\\ \multicolumn{1}{c}{-\epsilon \mathrm{\Phi }(t)}\end{array}$$

(19)

Therefore, the proposed NDOB1 can be formulated as follows ⁴⁷ :

$$\{\begin{array}{c}\stackrel{·}{\hat{\mathrm{\Phi }}}(t)=-\epsilon \hat{\mathrm{\Phi }}(t)-\epsilon [\epsilon X_5+h(X_5)(W_{\mathrm{cp}}(X_4,X_5)-b_{17}{X}_5)]\hfill \\ \multicolumn{1}{c}{\hat{\psi }(t)=\hat{\mathrm{\Phi }}(t)+\epsilon X_5}\end{array}$$

(20)

Remark 2. According to NDOB formula ⁴⁷ the defined variable $\mathrm{\Phi }(t)$ in (19) is used to transfer (18) into general form of NDOB1 (19).

Then, by using the proposed NDBO1 (20), the estimated value of $\psi (\chi ,X_5)$ can be calculated. According to Remark 1, and thus through the reciprocal function $h(X_5)$ , we can get the estimated value of $\chi (t)$ as follows:

$$\hat{\chi }=b_{17}^{-1}h(X_5{)}^{-1}\hat{\psi }(\chi ,X_5)$$

(21)

where $\hat{\chi }(t)$ is the estimation of CP. Accordingly, the OER can be exactly reconstructed as follows:

$${\hat{\lambda }}_{O_2}=\frac{b_{22}(X_5-\hat{\chi })}{b_{23}D}$$

(22)

Theorem 1. Considering the suggested NDOB1 (20) for the system dynamics (8), the estimating error $\tilde{\psi }(t)$ can converge to a bounded region such that $|\tilde{\psi }(t)|\le \frac{\eta }{\epsilon }$ , where $\epsilon >0$ .

Proof. From (19) and (20), we have $\tilde{\mathrm{\Phi }}(t)=\mathrm{\Phi }(t)-\hat{\mathrm{\Phi }}(t)$ , $\stackrel{·}{\stackrel{~}{\mathrm{\Phi }}}(t)=\stackrel{·}{\psi }(t)-\epsilon \tilde{\mathrm{\Phi }}(t)$ , and $\tilde{\psi }(t)=\tilde{\mathrm{\Phi }}(t)$ . Then, the Lyapunov function is defined as follows:

$$V_1(t)=\frac{1}{2}{\tilde{\psi }}^2(t)$$

(23)

Differentiating (23), one can get

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_1(t)=\tilde{\psi }(t)\stackrel{·}{\stackrel{~}{\psi }}(t)=\tilde{\mathrm{\Phi }}(t)\stackrel{·}{\stackrel{~}{\mathrm{\Phi }}}(t)}\\ \multicolumn{1}{c}{=\tilde{\mathrm{\Phi }}(t)[\stackrel{·}{\psi }(t)-\epsilon \tilde{\mathrm{\Phi }}(t)]}\\ \multicolumn{1}{c}{=-\epsilon {\tilde{\mathrm{\Phi }}}^2(t)+\tilde{\mathrm{\Phi }}\stackrel{·}{\psi }(t)}\end{array}$$

(24)

According to Lemma 1, the part $\tilde{\mathrm{\Phi }}(t)\stackrel{·}{\psi }(t)$ in (24) is given as $|\tilde{\mathrm{\Phi }}(t)\stackrel{·}{\psi }(t)|\le \frac{1}{2}\epsilon {\overline{\mathrm{\Phi }}}^2(t)+\frac{1}{2\epsilon }{\stackrel{·}{\psi }}^2(t)$ . Therefore, (24) can be rewritten as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_1(t)\le -\epsilon {\tilde{\mathrm{\Phi }}}^2(t)+\frac{1}{2}\epsilon {\tilde{\mathrm{\Phi }}}^2(t)+\frac{1}{2\epsilon }{\stackrel{·}{\psi }}^2(t)}\\ \multicolumn{1}{c}{\le -\frac{1}{2}\epsilon {\tilde{\mathrm{\Phi }}}^2(t)+\frac{1}{2\epsilon }{\stackrel{·}{\psi }}^2(t)}\\ \multicolumn{1}{c}{\le -\frac{1}{2}\epsilon {\tilde{\mathrm{\Phi }}}^2(t)+\frac{1}{2\epsilon }{\eta }^2}\\ \multicolumn{1}{c}{\le -\epsilon V_1(t)+\frac{1}{2\epsilon }{\eta }^2}\end{array}$$

(25)

Next, by integrating both sides of (25), we obtain:

$$\ln \{-\epsilon V_1(t)+\frac{1}{2\epsilon }{\eta }^2\}\le -\epsilon t$$

(26)

By simplifying (26), one has

$$V_1(t)\le \frac{e^{-\epsilon t}-\frac{1}{2\epsilon }{\eta }^2}{-\epsilon }$$

(27)

Substituting (23) into (27) yields to:

$$|\tilde{\psi }(t)|\le \sqrt{\frac{2e^{-\epsilon t}-\frac{1}{\epsilon }{\eta }^2}{-\epsilon }}$$

(28)

According to the result that $e^{-\epsilon t}\to 0$ with time increasing $t$ , one can conclude that the convergence of estimating error $|\tilde{\psi }(t)|\le \frac{\eta }{\epsilon }$ for $\epsilon >0$ is ensured.

Observer-based- $\alpha$ -variable model-free super-twisting fractional-order sliding mode control

This part introduces an $\alpha (t)$ -MF-STFOSMC design, whose control structure is directed in Figure 3. The proposed controller is composed of NDOB2-based- $\alpha$ -variable ULM algorithm and STFOSMC. The NDOB2 is used to estimate uncertain system dynamics, and ST-FOSMC is constructed to achieve optimal performance with less input chattering, rapid convergence speed, and finite-time stabilization. First, the design principle of NDOB-based- $\alpha$ -fixed model-free intelligent proportional integral controller (NDOB-MF-iPI) is presented and followed by $\alpha (t)$ -MF-STFOSMC design. Moreover, by using the Lyapunov theorem, the stability of the proposed control strategy via a closed-loop system is investigated.

NDOB-MF-iPI design

For simplicity’s sake to system with a single control variable ${\nu }_{\mathrm{cm}}(t)$ and a single output variable ${\lambda }_{O_2}(t)$ , according to the basic principle of the ultra-local model (ULM), the complex nonlinear PEMFC air supply system (8) can be re-formulated and decoupled by the ULM as follows ¹ :

$${\lambda }_{O_2}^{(\iota )}(t)=\alpha {\nu }_{\mathrm{cm}}(t)+\mathrm{\Upsilon }(t)$$

(29)

where ${\lambda }_{O_2}^{(\iota )}(t)$ is the $i$ th-order derivative of the system output (oxygen excess ratio), $\mathrm{\Upsilon }(t)$ denotes the continuously updated lumped uncertainties which is estimated via the measurements of the control input ${\nu }_{\mathrm{cm}}(t)$ and the controlled output ${\lambda }_{O_2}(t)$ . The term $\mathrm{\Upsilon }(t)$ does not distinguish between system nonlinearities, unknown dynamics, and external disturbances. $\alpha$ is a positive parameter which only indicates ${\nu }_{\mathrm{cm}}(t)$ and ${\lambda }_{O_2}^{(\iota )}(t)$ are of the same magnitude. This approach avoids the dependence of controller design on accurate model information and greatly reduces the difficulty of controller design.

Remark 3. According to (8), the degree between the control input ${\nu }_{\mathrm{cm}}(t)$ and controlled output ${\lambda }_{O_2}(t)$ is 2, the functions $F_i(X)$ ( $i$  = 1,2, …, 5.) will be highly coupled and complex. This makes the design of model-based controllers by using the knowledge of $F_i(X)$ difficult to design and implement. Moreover, the high order derivatives for the functions $F_i(X)$ will lead to high and unknown nonlinearities. In the existing research, these high and unknown nonlinearities are always approximated and eliminated by using intelligent algorithm, such as fuzzy logic and neural networks. However, the use of intelligent algorithm increases the computational burden on the controller. Therefore, due to advantages of the ULM that only depends on the control input and output variables, the ULM algorithm in (29) is used to design the controller. Hence, $\iota$ always be chosen 1 or 2. In this article, the case where $\iota$  = 1 is considered.

Thus, by selecting $\iota$  = 1, the NDBO-MF-iPI can be designed as follows ¹ :

$${\nu }_{\mathrm{cm}}(t)=\frac{1}{\alpha }[-\hat{\mathrm{\Upsilon }}(t)+{\stackrel{·}{\lambda }}_{O_2}^{*}+K_{\mathrm{p}}e(t)+K_{\mathrm{i}}\int e(t)\mathrm{dt}]$$

(30)

where ${\lambda }_{O_2}^{*}$ represents the reference trajectory, $K_{\mathrm{p}}$ and $K_{\mathrm{i}}$ are the coefficients of proportional and integral, respectively. $e(t)={\lambda }_{O_2}^{*}-{\lambda }_{O_2}(t)$ represents the tracking error and $\hat{\mathrm{\Upsilon }}(t)$ denotes the estimated value of unknown system dynamics $\mathrm{\Upsilon }(t)$ . Substituting (30) into (31), the error equation can be defined as

$$\stackrel{·}{e}(t)+K_{\mathrm{p}}e(t)+K_{\mathrm{i}}\int e(t)\mathrm{dt}-\tilde{\mathrm{\Upsilon }}(t)=0$$

(31)

where $\tilde{\mathrm{\Upsilon }}(t)=\hat{\mathrm{\Upsilon }}(t)-\mathrm{\Upsilon }(t)$ denotes the estimation error.

The steady state error can be guaranteed, if the values of $K_{\mathrm{p}}$ and $K_{\mathrm{i}}$ are appropriately selected according to the Hurwitz criterion. ¹

Nonlinear disturbance observer-2 (NDOB2) for estimation of uncertain system dynamics

In this paper, NDOB2 is used to estimate the lumped uncertainties $\mathrm{\Upsilon }(t)$ via the control input and output data.

Assumption 2. There exists a constant $\varrho >0$ such as the derivative $\stackrel{·}{\mathrm{\Upsilon }}(t)$ of the uncertain system dynamics term $\mathrm{\Upsilon }(t)$ satisfies $|\stackrel{·}{\mathrm{\Upsilon }}(t)|\le \varrho$ , $\varrho >0$ is positive constant.

The Assumption 2 is considered to theoretically ensure the existence of the constant $\varrho$ for the derivative $\stackrel{·}{\mathrm{\Upsilon }}(t)$ of the uncertain system dynamics term $\mathrm{\Upsilon }(t)$ .

Let’s assume that the estimated value of the unknown system dynamics and its estimation error are defined as $\hat{\mathrm{\Upsilon }}(t)$ and $\tilde{\mathrm{\Upsilon }}(t)=\mathrm{\Upsilon }(t)-\hat{\mathrm{\Upsilon }}(t)$ , respectively. Then, define a new vertical variable as $\theta (t)=\mathrm{\Upsilon }(t)-\mu {\lambda }_{O_2}(t)$ , where $\mu >0$ is positive constant.

Next, by differentiating $\theta (t)$ , one can get

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\theta }(t)=\stackrel{·}{\mathrm{\Upsilon }}(t)-\mu {\stackrel{·}{\lambda }}_{O_2}(t)}\\ \multicolumn{1}{c}{=\stackrel{·}{\mathrm{\Upsilon }}(t)-\mu [\alpha u(t)+\mathrm{\Upsilon }(t)]}\\ \multicolumn{1}{c}{=\stackrel{·}{\mathrm{\Upsilon }}(t)-\mu [\mu {\lambda }_{O_2}(t)+\alpha {\nu }_{\mathrm{cm}}(t)]-\mu \theta (t)}\end{array}$$

(32)

Hence, the suggested NDOB2 can be constructed as bellow ⁴⁷ :

$$\{\begin{array}{c}\stackrel{·}{\hat{\theta }}(t)=-\mu \hat{\theta }(t)-\mu [\mu {\lambda }_{O_2}(t)+\alpha {\nu }_{\mathrm{cm}}(t)]\hfill \\ \multicolumn{1}{c}{\hat{\mathrm{\Upsilon }}(t)=\hat{\theta }(t)+\mu {\lambda }_{O_2}(t)}\end{array}$$

(33)

Theorem 2. Considering the suggested NDOB2 (33) for the system dynamics (29), the estimating error $\tilde{\mathrm{\Upsilon }}(t)$ can converge to bounded region as $|\tilde{\mathrm{\Upsilon }}(t)|\le \frac{\varrho }{\mu }$ for $\mu >0$ .

Proof. From (32) and (33), we have $\tilde{\theta }(t)=\theta (t)-\hat{\theta }(t)$ , $\stackrel{·}{\stackrel{~}{\theta }}(t)=\stackrel{·}{\mathrm{\Upsilon }}(t)-\mu \tilde{\theta }(t)$ , and $\tilde{\mathrm{\Upsilon }}(t)=\tilde{\theta }(t)$ .

Then, let’s define the Lyapunov function as:

$$V_2(t)=\frac{1}{2}{\tilde{\mathrm{\Upsilon }}}^2(t)$$

(34)

By differentiating (34), we can get

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_2(t)=\tilde{\mathrm{\Upsilon }}(t)\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)=\tilde{\theta }(t)\stackrel{·}{\stackrel{~}{\theta }}(t)}\\ \multicolumn{1}{c}{=\tilde{\theta }(t)[\stackrel{·}{\mathrm{\Upsilon }}(t)-\mu \tilde{\theta }(t)]}\\ \multicolumn{1}{c}{=-\mu {\tilde{\theta }}^2(t)+\tilde{\theta }(t)\stackrel{·}{\mathrm{\Upsilon }}(t)}\end{array}$$

(35)

According to Lemma 1, the part $\tilde{\theta }(t)\stackrel{·}{\mathrm{\Upsilon }}(t)$ in (35) is given as $|\tilde{\theta }(t)\stackrel{·}{\mathrm{\Upsilon }}(t)|\le \frac{1}{2}\mu {\tilde{\theta }}^2(t)+\frac{1}{2\mu }{\stackrel{·}{\mathrm{\Upsilon }}}^2(t)$ . Thus, (35) can be rewritten as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_2(t)\le -\mu {\tilde{\theta }}^2(t)+\frac{1}{2}\mu {\tilde{\theta }}^2(t)+\frac{1}{2\mu }{\stackrel{·}{\mathrm{\Upsilon }}}^2(t)}\\ \multicolumn{1}{c}{\le -\frac{1}{2}\mu {\tilde{\theta }}^2(t)+\frac{1}{2\mu }{\stackrel{·}{\mathrm{\Upsilon }}}^2(t)}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{\le -\frac{1}{2}\mu {\tilde{\theta }}^2(t)+\frac{1}{2\mu }{\varrho }^2}\\ \multicolumn{1}{c}{\le -\mu V_2(t)+\frac{1}{2\mu }{\varrho }^2}\end{array}$$

(36)

Next, by integrating both sides of (36), we obtain:

$$\ln \{-\mu V_2(t)+\frac{1}{2\mu }{\varrho }^2\}\le -\mu t$$

(37)

By simplifying (37), we can get

$$V_2(t)\le \frac{e^{-\mu t}-\frac{1}{2\mu }{\varrho }^2}{-\mu }$$

(38)

By substituting (34) into (38), one has

$$|\tilde{\mathrm{\Upsilon }}(t)|\le \sqrt{\frac{2e^{-\mu t}-\frac{1}{\mu }{\varrho }^2}{-\mu }}$$

(39)

According to the result that $e^{-\mu t}\to 0$ with time increasing of $t$ , one can conclude that the convergence of estimating error $|\tilde{\mathrm{\Upsilon }}(t)|\le \frac{\varrho }{\mu }$ for $\mu >0$ is ensured.

$\alpha (t)$ -MF-STFOSMC design

According to the ULM principle, (29) can be rewritten as follows

$${\stackrel{·}{\lambda }}_{O_2}(t)=\alpha (t){\nu }_{\mathrm{cm}}(t)+\mathrm{\Upsilon }(t)$$

(40)

where $\alpha (t)$ is updated parameter and its discerption will be given in (51) later.

Then, the control law of $\alpha (t)$ -MF-STFOSMC can be constructed as follows:

$${\nu }_{\mathrm{cm}}(t)=\frac{{\stackrel{·}{\lambda }}_{O_2}^{*}(t)-\hat{\mathrm{\Upsilon }}(t)+{\nu }_{\mathrm{cm};\mathrm{TEC}}(t)}{\alpha (t)}$$

(41)

where ${\nu }_{\mathrm{cm};\mathrm{TEC}}(t)$ denotes the sub-control law for tracking error convergence (TEC) to be constructed and $\hat{\mathrm{\Upsilon }}(t)$ is estimated value of uncertain system dynamics using the designed NDOB2 (33). Note that the $\alpha (t)$ variable can be updated in (33) accordingly.

By substituting (41) into (40), a new error equation can be obtained as follows:

$$\stackrel{·}{e}(t)-\tilde{\mathrm{\Upsilon }}(t)+{\nu }_{\mathrm{cm};\mathrm{TEC}}(t)=0$$

(42)

Now, new state variables are defined as $e_1(t)=\int e(t)\mathrm{dt}$ and $e_2(t)=e(t)$ ; then, (42) can be expressed in a state space form as follows:

$$\{\begin{array}{c}{\stackrel{·}{e}}_1(t)=e_2(t)\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{e}}_2(t)=\tilde{\mathrm{\Upsilon }}(t)-{\nu }_{\mathrm{cm};\mathrm{TEC}}(t)}\end{array}$$

(43)

In order to achieve finite-time convergence for tracking error with less input chattering, a fractional-order non-singular fast terminal sliding surface is designed as follows:

$$\sigma (t)={\stackrel{·}{e}}_1(t)+m_1{e}_1^{p/q}(t)+m_2{D}^{ϵ-1}\mathrm{sgn}(e_1(t){)}^{\kappa }$$

(44)

where $1<p/q<2$ , $ϵ>0$ , and $m_1,m_2>0$ are positive constants, the notation $\mathrm{sgn}(e_1(t){)}^{\kappa }=|e_1(t){|}^{\kappa }\mathrm{sgn}(e_1(t))$ with $0<\kappa <1$ . By taking the first time derivative for (44) and using (43), one can obtain

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\sigma }(t)=\tilde{\mathrm{\Upsilon }}(t)-{\nu }_{\mathrm{cm};\mathrm{TEC}}(t)+m_1\frac{p}{q}{e}_1^{(p-q)/q}{\stackrel{·}{e}}_1(t)}\\ \multicolumn{1}{c}{+m_2{D}^{ϵ}\mathrm{sgn}(e_1(t){)}^{\kappa }}\end{array}$$

(45)

The above proposed switching surface has fast convergence speed and high tracking accuracy. Then, by using ST-FOSMC, the sub-control law ${\nu }_{\mathrm{cm};\mathrm{TEC}}(t)$ can be formulated as bellow:

$${\nu }_{\mathrm{cm};\mathrm{TEC}}(t)={\nu }_{\mathrm{cm};\mathrm{TEC}}^{\mathrm{eq}}(t)+{\nu }_{\mathrm{cm};\mathrm{TEC}}^{\mathrm{re}}(t)$$

(46)

where ${\nu }_{\mathrm{cm};\mathrm{TEC}}^{\mathrm{eq}}(t)$ denotes the equivalent controller law and its corresponding equation is defined as follows:

$${\nu }_{\mathrm{cm};\mathrm{TEC}}^{\mathrm{eq}}(t)=m_1\frac{p}{q}{e}_1^{(p-q)/q}(t){\stackrel{·}{e}}_1(t)+m_2{D}^{ϵ}\mathrm{sgn}(e_1(t){)}^{\kappa }$$

(47)

Simultaneously, the reaching controller law ${\nu }^{\mathrm{re}}\mathrm{cm};\mathrm{TEC}(t)$ is formulated using the ST algorithm to ensure the reaching of the switching surface with reduced chattering and improved control performance. Consequently, the establishment of ${\nu }^{\mathrm{re}}\mathrm{cm};\mathrm{TEC}(t)$ is as follows:

$${\nu }_{\mathrm{cm};\mathrm{TEC}}^{\mathrm{re}}(t)=k_1\mathrm{sgn}(\sigma (t){)}^{1/2}+k_2\int \mathrm{sgn}(\sigma (t))\mathrm{dt}$$

(48)

where $k_1$ and $k_2$ are positive constants.

Correspondingly, from (46) to (48), the tracking error convergence law ${\nu }_{\mathrm{cm};\mathrm{TEC}}(t)$ is given as:

$$\begin{array}{c}\multicolumn{1}{c}{{\nu }_{\mathrm{cm};\mathrm{TEC}}(t)=m_1\frac{p}{q}{e}_1^{(p-q)/q}(t){\stackrel{·}{e}}_1(t)+m_2{D}^{ϵ}\mathrm{sgn}(e_1(t){)}^{\kappa }}\\ +k_1\mathrm{sgn}(\sigma (t){)}^{1/2}+k_2\int \mathrm{sgn}(\sigma (t))\mathrm{dt}\hfill \end{array}$$

(49)

Finally, based on NDOB2 (33) and ST-FOSMC law (49), the actual control law of the proposed $\alpha (t)$ -MF-STFOSMC controller is formulated as follows:

$$\begin{array}{c}\multicolumn{1}{c}{{\nu }_{\mathrm{cm}}(t)=\frac{1}{\alpha (t)}[m_1\frac{p}{q}{e}_1^{\frac{(p-q)}{q}}{\stackrel{·}{e}}_1(t)+m_2{D}^{ϵ}\mathrm{sgn}{(e_1(t))}^{\kappa }}\\ +k_1\mathrm{sgn}(\sigma (t){)}^{1/2}+k_2\int \mathrm{sgn}(\sigma (t))\mathrm{dt}-\hat{\mathrm{\Upsilon }}(t)+{\stackrel{·}{\lambda }}_{O_2}^{*}]\hfill \end{array}$$

(50)

In order to further improve the tracking performance of the proposed $\alpha (t)$ -MF-STFOSMC method, an $\alpha (t)$ -variable approach is proposed to automatically update the value of $\alpha$ and its correspondingly formula is expressed as follows:

$$\stackrel{·}{\alpha }(t)=-\mathrm{\Omega }|\sigma (t){|}^M\mathrm{sgn}(\alpha (t)-{\alpha }_{\min })$$

(51)

where ${\alpha }_{\min }$ represents the positive lower bound of $\alpha (t)$ , $\mathrm{\Omega }$ and $M$ are two positive designed parameters. Due to that $\alpha (t)=0$ will cause the singularity in control input ${\nu }_{\mathrm{cm}}(t)$ ; therefore, the part of $\mathrm{sgn}(\alpha (t)-{\alpha }_{\min })$ is formulated to ensure that $\alpha (t)$ is not less than ${\alpha }_{\min }$ .

Theorem 3. Considering the nonlinear PEMFC air supply system (8) reformulated by ultra-local model (40), the control law proposed in (50) with the NDOB2 (33) and the ST algorithm (48) can ensure that the sliding mode surface $\sigma (t)$ converge to 0 in finite-time.

Proof. Substituting (49) into (45), one gets

$$\begin{array}{c}\stackrel{·}{\sigma }(t)=-k_1\mathrm{sgn}{(\sigma (t))}^{\frac{1}{2}}+\tau (t),\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\tau }(t)=-k_2\mathrm{sgn}(\sigma (t))-\omega (t)}\end{array}$$

(52)

where $\omega (t)=\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)$ .

Then, let’s define a new vector as:

$$\zeta (t)=[{\zeta }_1(t){\zeta }_2(t){]}^{\mathrm{T}}=[|\sigma (t){|}^{1/2}\mathrm{sgn}(\sigma (t))\tau (t){]}^{\mathrm{T}}$$

(53)

Next, the following equation (54) describe the Lyapunov candidate function.

$$V(t)=V_0(t)+\frac{1}{2}{\tilde{\alpha }}^2(t)$$

(54)

where $\tilde{\alpha }(t)=\alpha (t)-{\alpha }_{\min }$ and $V_0(t)$ is given as:

$$V_0(t)={\zeta }^{\mathrm{T}}(t)\mathrm{\Xi }\zeta (t)$$

(55)

where $\mathrm{\Xi }=\left[\begin{array}{cc}{\mathrm{\Xi }}_{11}\hfill & {\mathrm{\Xi }}_{12}\hfill \\ \multicolumn{1}{c}{*}& {\mathrm{\Xi }}_{22}\hfill \end{array}\right]$ is positive definite, namely,

$${\lambda }_{\min }(\mathrm{\Xi })\parallel \zeta (t){\parallel }^2\le V_0(t)\le {\lambda }_{\max }(\mathrm{\Xi })\parallel \zeta (t){\parallel }^2$$

(56)

where ${\lambda }_{\min }(\mathrm{\Xi })$ and ${\lambda }_{\max }(\mathrm{\Xi })$ are minimum and maximum eigenvalues of $\mathrm{\Xi }$ and $\parallel \zeta (t){\parallel }^2=|\sigma (t)|+{\tau }^2(t)$ denotes the Euclidean norm of $\zeta (t)$ .

By differentiating the $\zeta (t)$ , we can get

$$\stackrel{·}{\zeta }(t)=[{\stackrel{·}{\zeta }}_1(t){\stackrel{·}{\zeta }}_2(t){]}^{\mathrm{T}}={[\frac{\stackrel{·}{\sigma }(t)}{2|\sigma (t{)|}^{1/2}}\stackrel{·}{\tau }(t)]}^{\mathrm{T}}$$

(57)

Based on (52) and (57), one has

$$\stackrel{·}{\zeta }(t)=\frac{A\zeta (t)}{|\sigma (t{)|}^{1/2}}+B\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)$$

(58)

with $A=\left[\begin{array}{cc}-\frac{k_1}{2}\hfill & \frac{1}{2}\hfill \\ \multicolumn{1}{c}{-k_2}& 0\hfill \end{array}\right]$ and $B=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{-1}\end{array}\right]$ .

Then, taking the time derivative of $V(t)$ yields to:

$${\stackrel{·}{V}}_0(t)={\stackrel{·}{\zeta }}^{\mathrm{T}}(t)\mathrm{\Xi }\zeta (t)+{\zeta }^{\mathrm{T}}(t)\mathrm{\Xi }\stackrel{·}{\zeta }(t)$$

(59)

Substituting (58) into (59), one obtains

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_0(t)=\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}(A^{\mathrm{T}}\mathrm{\Xi }+\mathrm{\Xi }A)\zeta (t)}\\ \multicolumn{1}{c}{+2[0-\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)]P\zeta (t)}\\ \multicolumn{1}{c}{=\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}(A^{\mathrm{T}}\mathrm{\Xi }+\mathrm{\Xi }A)\zeta (t)}\\ \multicolumn{1}{c}{+2[0-|\sigma (t{)|}^{1/2}\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)]P\zeta (t)}\\ \multicolumn{1}{c}{=\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}(A^{\mathrm{T}}\mathrm{\Xi }+\mathrm{\Xi }A)\zeta (t)}\\ \multicolumn{1}{c}{+\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}({\mathrm{\Pi }}^{\mathrm{T}}\mathrm{\Xi }+\mathrm{\Xi }\mathrm{\Pi })\zeta (t)}\\ \multicolumn{1}{c}{=\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}\{{(A+\mathrm{\Pi })}^{\mathrm{T}}\mathrm{\Xi }+\mathrm{\Xi }(A+\mathrm{\Pi })\}\zeta (t)}\end{array}$$

(60)

with $\mathrm{\Pi }=\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{-\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)\mathrm{sgn}(\sigma (t))}& 0\hfill \end{array}\right]$ . Then, one can define

$$\begin{array}{c}\multicolumn{1}{c}{Q=-(A+\mathrm{\Pi }{)}^{\mathrm{T}}\mathrm{\Xi }-\mathrm{\Xi }(A+\mathrm{\Pi })}\\ \multicolumn{1}{c}{=\left[\begin{array}{cc}\frac{k_1}{2}\hfill & \varpi (t)\hfill \\ \multicolumn{1}{c}{-\frac{1}{2}}& 0\hfill \end{array}\right]\left[\begin{array}{cc}{\mathrm{\Xi }}_{11}\hfill & {\mathrm{\Xi }}_{12}\hfill \\ \multicolumn{1}{c}{*}& {\mathrm{\Xi }}_{22}\hfill \end{array}\right]}\\ \multicolumn{1}{c}{+\left[\begin{array}{cc}{\mathrm{\Xi }}_{11}\hfill & {\mathrm{\Xi }}_{12}\hfill \\ \multicolumn{1}{c}{*}& {\mathrm{\Xi }}_{22}\hfill \end{array}\right]\left[\begin{array}{cc}\frac{k_1}{2}\hfill & -\frac{1}{2}\hfill \\ \multicolumn{1}{c}{\varpi (t)}& 0\hfill \end{array}\right]}\\ \multicolumn{1}{c}{=\left[\begin{array}{cc}{k}_1{\mathrm{\Xi }}_{11}+2\varpi (t){\mathrm{\Xi }}_{12}\hfill & -\frac{1}{2}{\mathrm{\Xi }}_{11}+\frac{1}{2}{k}_1{\mathrm{\Xi }}_{12}+\varpi (t){\mathrm{\Xi }}_{22}\hfill \\ \multicolumn{1}{c}{*}& -{\mathrm{\Xi }}_{12}\hfill \end{array}\right]}\end{array}$$

(61)

with $\varpi (t)=k_2+\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)\mathrm{sgn}(\sigma (t))$ . Since $\stackrel{·}{\stackrel{~}{\mathrm{\Upsilon }}}(t)$ is limited, and there exists $\delta \ge 0$ such that $k_2-\delta \le \varpi (t)\le k_2+\delta$ . Furthermore, if parameter $k_2>0$ is selected appropriately, we can ensure that $\varpi (t)>0$ . For convenience, let ${\mathrm{\Xi }}_{11}$ be defined as 1, which leads to the following expressions for $\mathrm{\Xi }$ and $Q$ :

$$\begin{array}{l}\mathrm{\Xi }=\left[\begin{array}{cc}1& {\mathrm{\Xi }}_{12}\\ \ast & {\mathrm{\Xi }}_{22}\end{array}\right]\\ Q=\left[\begin{array}{cc}{k}_1+2\varpi (t){\mathrm{\Xi }}_{12}& -\frac{1}{2}+\frac{1}{2}{k}_1{\mathrm{\Xi }}_{12}+\varpi (t){\mathrm{\Xi }}_{22}\\ \ast & -{\mathrm{\Xi }}_{12}\end{array}\right]\end{array}$$

(62)

Therefore, the following inequalities (63) and (64) should hold to make $\mathrm{\Xi }$ and $Q$ both positive definite,

$${\mathrm{\Xi }}_{22}>{\mathrm{\Xi }}_{12}^2,{\mathrm{\Xi }}_{12}<0$$

(63)

$$\begin{array}{c}\multicolumn{1}{c}{k_1{\mathrm{\Xi }}_{12}+2\varpi (t){\mathrm{\Xi }}_{12}^2+\frac{1}{4}(k_1{\mathrm{\Xi }}_{12}-1{)}^2-(1-k_1{\mathrm{\Xi }}_{12})\varpi (t){\mathrm{\Xi }}_{22}}\\ \multicolumn{1}{c}{+{\varpi }^2(t){\mathrm{\Xi }}_{22}^2<0}\end{array}$$

(64)

Therefore, the inequality (64) will be satisfied if

$$\begin{array}{l}{k}_1{\mathrm{\Xi }}_{12}+2(k_2+\delta ){\mathrm{\Xi }}_{12}^2+\frac{1}{4}{(k_1{\mathrm{\Xi }}_{12}-1)}^2\\ -(1-k_1{\mathrm{\Xi }}_{12})(k_2-\delta ){\mathrm{\Xi }}_{22}+{(k_2+\delta )}^2{\mathrm{\Xi }}_{22}^2<0\end{array}$$

(65)

It can follow that (63) and (64) will be satisfied if the following inequalities (66) and (67) are satisfied:

$$\begin{array}{c}\mathcal{M}>0,\mathcal{N}>0,\mathcal{C}>0,0<\mathfrak{J}<1\hfill \end{array}$$

(66)

$$\mathcal{M}-\frac{2}{\mathcal{N}}\mathcal{C}>\frac{1}{4}{(1+\mathcal{M})}^2-\mathfrak{J}(1+\mathcal{M})\mathcal{C}+{\mathcal{C}}^2$$

(67)

where

$$\mathcal{M}\stackrel{\mathrm{\Delta }}{=}-k_1{\mathrm{\Xi }}_{12},\mathcal{N}\stackrel{\mathrm{\Delta }}{=}\frac{{\mathrm{\Xi }}_{22}}{{\mathrm{\Xi }}_{12}^2},\mathcal{C}\stackrel{\mathrm{\Delta }}{=}{\mathrm{\Xi }}_{22}(k_2+\delta ),\mathfrak{J}\stackrel{\mathrm{\Delta }}{=}\frac{k_2-\delta }{k_2+\delta }$$

(68)

The inequality (66) denotes the interior of an ellipse within the system $(\mathcal{M},\mathcal{C})$ , which is parameterized by $(\mathcal{M},\mathcal{C})$ and its center located at the point as follows:

$${\mathcal{M}}_c=\frac{(\mathfrak{J}-\frac{2}{\mathcal{N}})\mathfrak{J}+1}{1-{\mathfrak{J}}^2},{\mathcal{C}}_c=\frac{\mathfrak{J}-\frac{1}{\mathcal{N}}}{1-{\mathfrak{J}}^2}$$

(69)

If $\mathcal{N}\mathfrak{J}>1$ , it is always possible to find positive values of $(\mathcal{M},\mathcal{C})$ that satisfy (67). In this case, the center $({\mathcal{M}}_c,{\mathcal{C}}_c)$ of the ellipse is located in the first quadrant, that is, ${\mathcal{M}}_c>0$ , ${\mathcal{C}}_c>0$ .

Therefore, it can be concluded that by selecting appropriate values of $\mathcal{N}>1$ and $0<\mathfrak{J}<1$ such that $\mathcal{N}\mathfrak{J}>1$ , a positive definite matrix $\mathrm{\Xi }$ is obtained as bellow:

$$\mathrm{\Xi }=\left[\begin{array}{cc}1\hfill & -\sqrt{\frac{(1-\mathfrak{J})\mathcal{C}}{2\mathcal{N}\delta }}\hfill \\ \multicolumn{1}{c}{*}& \frac{(1-\mathfrak{J})\mathcal{C}}{2\delta }\hfill \end{array}\right]$$

(70)

and the positive definite matrix $Q$ is ensured.

Next, from (60) and (61), we have

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_0(t)=-\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}Q\zeta (t)}\\ \multicolumn{1}{c}{\le -\frac{1}{|\sigma (t{)|}^{1/2}}{\zeta }^{\mathrm{T}}{\lambda }_{\min }\left\{Q\right\}\parallel \zeta (t)\parallel }\end{array}$$

(71)

where ${\lambda }_{\min }\left\{Q\right\}$ represents the minimum eigenvalue of matrix $Q$ .

Further, using (56) and (71) and from the fact that $|\sigma (t){|}^{1/2}\le \parallel \zeta (t)\parallel \le \frac{V_0^{1/2}(t)}{{\lambda }_{\min }^{1/2}\left\{\mathrm{\Xi }\right\}}$ , we can get

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_0(t)\le -\frac{{\lambda }_{\min }^{1/2}\left\{\mathrm{\Xi }\right\}{\lambda }_{\min }\left\{Q\right\}}{{\lambda }_{\max }^{1/2}\left\{\mathrm{\Xi }\right\}}{V}_0^{1/2}(t)}\\ \multicolumn{1}{c}{\le -{\vartheta }_1{V}_0^{1/2}(t)}\end{array}$$

(72)

with ${\vartheta }_1=\frac{{\lambda }_{\min }^{1/2}\left\{\mathrm{\Xi }\right\}{\lambda }_{\min }\left\{Q\right\}}{{\lambda }_{\max }^{1/2}\left\{\mathrm{\Xi }\right\}}$ is positive.

According to the widely accepted inequality $\sqrt{x^2+y^2}\le |x|+|y|$ , $\stackrel{·}{V}(t)$ in (54) can be derived as:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{V}(t)\le -{\vartheta }_1{V}_0^{1/2}(t)+\tilde{\alpha }(t)\stackrel{·}{\alpha }(t)}\\ \multicolumn{1}{c}{\le -{\vartheta }_1{V}_0^{1/2}(t)-\frac{{\vartheta }_2}{\sqrt{2}}\tilde{\alpha }(t)}\\ \multicolumn{1}{c}{\le -{\vartheta }_0{V}^{1/2}(t)}\end{array}$$

(73)

where ${\vartheta }_0=\min \{{\vartheta }_1,{\vartheta }_2\}$ , and ${\vartheta }_2=\sqrt{2}\mathrm{\Omega }|\sigma (t){|}^M$ . Therefore, it is oblivious from (72) that the origin along $\sigma (t)=0$ can be reached in a finite time.

According to Lemma 2, the finite time can be calculated as follows

$$t_{\mathrm{f}}\le \frac{2V^{1/2}(t_0)}{{\vartheta }_0}$$

(74)

This completes the proof.□

Remark 4. In order to obtain the continuous control input signal and avoid the chattering phenomenon, the $\mathrm{sgn}$ function in (50) and (51) is replaced with $\tanh (.)$ function in this work.

Remark 5. The parameters of the proposed $\alpha (t)$ -MF-STFOSMC method can be selected in accordance with specified range which is defined as $0<p/q<2$ , $m_1>0$ , $m_2>0$ , $k_1>0$ , $k_2>0$ , $0<\epsilon <1$ , $0<\kappa <1$ , and $\mu >0$ . Hence, if the controller’s parameters are not selected within the specified range, then the closed-loop system stability cannot be achieved. Thus, to obtain the stability of the closed-loop system desired trajectory tracking performance simultaneously, the appropriate value can be selected accordingly.

Simulation results and discussion

In this part, the robustness, efficiency, and superiority of the suggested controller are investigated, and a comparative study between the proposed $\alpha (t)$ -MF-STFOSMC and other controllers such as PID, NODOB-MF-iPI, and self-adaptive fuzzy PID (SFPID) ⁵³ is carried out in MATLAB/Simulink. Based on the different operating conditions, four cases are used to evaluate the performance of applied methods.

Operating condition for validation

• System model and controller parameters: The fifth-order PEMFC air supply system model (as described by equation (8)) is considered a real system and used in all simulation cases. Note that the designed controller based on the ultra-local model depends on the input and output signal of the controlled system (as described by equation (50) and shown in Figure 3). By try and error method, the parameters of applied controllers are selected as follows: PID as $K_{\mathrm{p}}=100$ , $K_{\mathrm{i}}=400$ , and $K_{\mathrm{d}}=0.9$ ; NDOB-MF-iPI as $K_{\mathrm{p}}=120,K_{\mathrm{i}}=650,\alpha =70.5,$ and $\mu =10.3$ ; the proposed $\alpha (t)$ -MF-STFOSMC as $q=0.7$ , $p=1.2$ , $m_1=42$ , $m_2=77.8$ , $k_1=25$ , $k_2=14$ , $ϵ=0.6$ , $\kappa =0.88$ , $\alpha (0)=70$ , $\epsilon =13.2$ , and $\mu =12$ .

• Performance indexes: The integral absolute error (IAE), the integral time-absolute error (ITAE), and the integral square error (ISE) are used to quantitatively assess the OER tracking error.

Results and analysis

Note that the main goal of the proposed $\alpha (t)$ -MF-STFOSMC method is to adjust the oxygen excess ratio ${\lambda }_{O_2}$ at its optimal value of reference trajectory ${\lambda }_{{\mathrm{O}}_2}^{\mathrm{Opt}}$ . In this paper, the simulation results are assessed in the following cases:

First case: Tracking performance with constant reference trajectory

In this case, the value of reference trajectory for ${\lambda }_{O_2}$ is kept constant at ${\lambda }_{O_2}^{*}=2$ , while different values of the disturbing current load $I_{\mathrm{st}}$ is applied (see Figure 4(a)), which cover the whole operational range of the nonlinear PEMFC system. The observation results of two NDOBs and the tracking performance of applied controllers are displayed in Figure 4. It can be observed that from Figure 4(b) and (c), the unmeasured CP and OER are effectively estimated by the proposed NDOB1 in the existence of the rapid time-varying load current, and the proposed $\alpha (t)$ -MF-STFOSMC significantly enhances the tracking accuracy and transient response speed. In addition, the proposed method has improved the values of rising time, settling time, overshoot, and undershoot compared with PID, NDOB-MF-iPI, and SFPID controllers. As shown in Figure 4(c) at $t=20$  s, a sudden increase in the current load will cause a sudden reduction of OER ${\lambda }_{O_2}$ and lead to a considerable drop in the stack voltage, as shown in Figure 4(d). The suggested $\alpha (t)$ -MF-STFOSMC method recovered from the impact of the perturbation faster, obtaining the least rising time and settling time with minimum overshoot and undershoot around the reference trajectory. Correspondingly, at $t=26$  s, a sudden decrease in the current load caused ${\lambda }_{O_2}$ to increase suddenly. This increase in ${\lambda }_{O_2}$ results in a decrease in the net power of the PEMFC system because of the increase in power consuming by the air compressor motor, as shown in Figure 4(e). Figure 4(f) illustrates the control inputs of applied controllers, and Figure 4(g) shows the estimated value of uncertain system dynamics which is obtained by proposed NDOB2. The sliding mode surface and variable $\alpha (t)$ of the proposed method are shown in Figure 4(h) and (i), respectively. The comparison of the total performance indexes is furnished in Table 2. The performance indexes in terms of IAE with the proposed $\alpha (t)$ -MF-STFOSMC is improved about 94% over PID controller, 79.5% over NDOB-MF-iPI controller, and 76% over SFPID controller. The enhancement of ISE value with the $\alpha (t)$ -MF-STFOSMC controller is 98% (PID controller), 82% (NDOB-MF-iPI controller), and 74% (SFPID controller), while the enhancement of ITAE value with the $\alpha (t)$ -MF-STFOSMC controller is 91% (PID controller), 80% (NDOB-MF-iPI controller), and 74% (SFPID controller). As it is obvious, the proposed $\alpha (t)$ -MF-STFOSMC method has a better dynamic response compared to other methods.

Figure 4.

(a) Current load profile (A), (b) estimation of cathode pressure (CP) by using NDOB1, (c) tracking performance of OER, (d) stack voltage (V), (e) net power (kW), (f) control input (voltage of air compressor motor (V)), (g) estimation of uncertain system dynamics by using NDOB2, (h) sliding mode surface of proposed method, and (i) the variable $\alpha (t)$ of proposed method.

Table 2.

Performance indexes: IAE, ISE, and ITAE for all cases.

Cases	Performance indexes	PID	NDOB-MF-iPI	SFPID	$\alpha (t)$ -MF-STFOSMC
First case	IAE	1.745	0.5073	0.4372	0.1037
	ISE	0.7186	0.0695	0.0485	0.0125
	ITAE	38.30	17.01	13.18	3.414
Second case	IAE	1.661	0.588	0.4901	0.0803
	ISE	0.7857	0.1668	0.1419	0.0092
	ITAE	32.07	12.45	11.74	2.511
Third case	IAE	2.754	1.305	1.224	0.5088
	ISE	0.8294	0.1858	0.1494	0.0239
	ITAE	66.31	35.46	32.49	16.54
Fourth case	IAE	15.69	5.054	3.798	1.32
	ISE	2.242	1.213	1.114	0.682
	ITAE	7660	2631	1879	626

Second case: Tracking performance with optimal reference trajectory ${\lambda }_{{\mathrm{O}}_2}^{\mathrm{Opt}}$

In this case, the proposed $\alpha (t)$ -MF-STFOSMC method is tested for tracking optimal reference trajectory ${\lambda }_{O_2}^{\mathrm{Opt}}$ to obtain the maximum power point operation of the PEMFC stack; wherein, the optimal value ${\lambda }_{{\mathrm{O}}_2}^{\mathrm{Opt}}$ (17) is determined based on the level of the load current $I_{\mathrm{st}}$ for achieving the maximum net power $P_{\mathrm{net}}$ from the PEMFC system. The same profile of the current load perturbation used in the first case is also used in the second case. However, the oxygen starvation risk due to sudden high power demands is increased because of the shortage of oxygen in the PEMFC stack. Therefore, the controller effort for regulating the optimal reference trajectory ${\lambda }_{{\mathrm{O}}_2}^{\mathrm{Opt}}$ is very crucial.

The observation results of two NDOBs and the tracking performance of applied controllers for regulating ${\lambda }_{O_2}$ are depicted in Figure 5. The simulation results illustrate that the proposed $\alpha (t)$ -MF-STFOSMC method has been successfully performed for reference trajectory tracking and disturbance rejection. The comparison of the total performance indexes is listed in Table 2. It can be seen from Table 2, that the proposed $\alpha (t)$ -MF-STFOSMC method has minimum values of performance indexes in terms of IAE, ISE, and ITAE compared to other methods. The performance indexes in terms of IAE with the proposed $\alpha (t)$ -MF-STFOSMC is improved about 95% over PID controller, 86% over NDOB-MF-iPI controller, and 83.6% over SFPID controller. The enhancement of ISE value with the $\alpha (t)$ -MF-STFOSMC controller is 98% (PID controller), 94% (NDOB-MF-iPI controller), and 93% (SFPID controller), while the enhancement of ITAE value with the $\alpha (t)$ -MF-STFOSMC controller is 92% (PID controller),79.8% (NDOB-MF-iPI controller), and 78% (SFPID controller). In comparison with the first case, the output voltage and output net power are shown in Figure 5(i) and (j), respectively. The obtained results exhibit that the second case (maximum power tracking) results in a higher net power output by 1% compared with the first case. According to the obtained results in the first and second cases, the suggested $\alpha (t)$ -MF-STFOSMC method provides reasonable tracking performance with higher efficiency compared to other approaches.

Figure 5.

(a) Estimation of cathode pressure (CP) by using NDOB1, (b) tracking performance of OER, (c) stack voltage (V), (d) net power (kW), (e) control input (voltage of air compressor motor (V)), (f) estimation of uncertain system dynamics by using NDOB2, (g) sliding mode surface of proposed method, (h) the variable $\alpha (t)$ of proposed method, (i) the comparison of the stack voltage between first case and second case, and (j) the comparison of the net power output between first case and second case.

Third case: Robustness analysis with parameter uncertainties and measurement noise

The main purpose of this subsection is to verify the robustness of the designed controller. In practical applications, some PEMFC parameters may vary with the operating conditions and external environment. Because we only used the MATLAB/Simulink environment to verify the performance of the proposed controller, parameter uncertainties and measurement noise are necessary to be considered for robustness testing of the proposed method. Therefore, we added slow time-varying perturbations for the system’s parameters, such as atmospheric temperature $T_{\mathrm{atm}}$ , stack temperature $T_{\mathrm{st}}$ , and relative humidity ${\varphi }_{\mathrm{atm}}$ , as shown in Figure 6. As recommended by Deng et al., ²⁴ a Gaussian noise is added to the measured variables $X_4$ , $X_5$ , and $I_{\mathrm{st}}$ . The simulation results of two NDOBs and applied controllers under the same operating conditions in the second case are shown in Figure 7. Naturally, when dealing with more complicated operating conditions, such as disturbed system parameters and measurement noise, different levels of performance degradation relative to the second case appear in all the examined methods. Compared to PID, NDOB-MF-iPI, and SFPID, the proposed $\alpha (t)$ -MF-STFOSMC method still provides reliable and robust tracking performance, and the resulting degradation is acceptable and slighter. Furthermore, a comparison of the total performance indexes such as IAE, ISE, and ITAE is given in Table 2. The performance indexes in terms of IAE with the proposed $\alpha (t)$ -MF-STFOSMC is improved about 81.5% over PID controller, 61% over NDOB-MF-iPI controller, and 58% over SFPID controller. The enhancement of ISE value with the $\alpha (t)$ -MF-STFOSMC controller is 97% (PID controller), 87% (NDOB-MF-iPI controller), and 84% (SFPID controller), while the enhancement of ITAE value with the $\alpha (t)$ -MF-STFOSMC controller is 75% (PID controller), 53% (NDOB-MF-iPI controller), and 49% (SFPID controller). Furthermore, it could be concluded that the proposed $\alpha (t)$ -MF-STFOSMC achieves stable and robust control, satisfying the controlled system performance, and there is no need for a reset even if the system is tested under a wide variation in the system’s parameters and loading conditions.

Figure 6.

Parameter uncertainties: (a) the atmospheric temperature $T_{\mathrm{atm}}$ (K), (b) stack temperature $T_{\mathrm{st}}$ , and (c) relative humidity ${\varphi }_{\mathrm{atm}}$ .

Figure 7.

(a) Estimation of cathode pressure (CP) by using NDOB1, (b) tracking performance of OER, (c) stack voltage (V), (d) net power (kW), (e) control input (voltage of air compressor motor (V)), (f) estimation of uncertain system dynamics by using NDOB2, (g) sliding mode surface of proposed method, and (h) the variable $\alpha (t)$ of proposed method.

Fourth case: Realistic drive cycle

It is necessary to investigate the tracking performance of the proposed method under a realistic drive cycle and time-varying reference trajectory. Thus, the effectiveness of the designed NDOB and controller for estimating and regulating OER is further tested by employing a real driving cycle ArtMw130. ⁵⁴ Since this paper does not focus on energy management, a simple filtering-based-energy management strategy for fuel cell vehicle hybrid power systems is used to calculate the current load profile. ⁵⁵ The vehicle speed and its corresponding current load profile are illustrated in Figure 8(a) and (b), respectively. Figure 8(c) to (j) illustrate the estimation of cathode pressure $P_{\mathrm{cp}}$ , tracking performance of OER, tracking error, stack voltage $V_{\mathrm{st}}$ , control input (voltage of air compressor motor ${\nu }_{\mathrm{cm}}$ ), estimation of uncertain system dynamics, sliding mode surface of proposed method and variable $\alpha (t)$ of proposed method, respectively. Table 2 shows the comparison of the total performance indexes for applied controllers. The performance indexes in terms of IAE with the proposed $\alpha (t)$ -MF-STFOSMC is improved about 91.5% over PID controller, 73.8% over NDOB-MF-iPI controller, and 65% over SFPID controller. The enhancement of ISE value with the $\alpha (t)$ -MF-STFOSMC controller is 69.5% (PID controller), 43% (NDOB-MF-iPI controller), and 38% (SFPID controller), while the enhancement of ITAE value with the $\alpha (t)$ -MF-STFOSMC controller 91% (PID controller), 76% (NDOB-MF-iPI controller), and 66% (SFPID controller). However, it can be concluded that the designed observer-based-controller strategy can achieve advantageous tracking performance and stabilize the PEMFC system under a complex and variable real driving cycle.

Figure 8.

(a) Vehicle speed (km/h), (b) current load profile (A), (c) estimation of cathode pressure (CP) by using NDOB1, (d) tracking performance of OER, (e) tracking error, (f) stack voltage (V), (g) control input (voltage of air compressor motor (V)), (h) estimation of uncertain system dynamics by using NDOB2, (i) sliding mode surface of proposed method, and (j) the variable $\alpha (t)$ of proposed method.

Conclusion

In this work, a novel $\alpha (t)$ -MF-STFOSMC controller based on NDOB is developed for the PEMFC system to accurately regulate the oxygen excess ratio in the presence of unmeasurable variables, external disturbance, uncertainties, and measurement noise. The designed NDOBs have a simple structure and can effectively estimate the unmeasurable variables and unmodeled PEMFC system dynamics. Moreover, a robust super-twisting fractional-order sliding mode control is formulated as a control law to achieve tracking performance. Compared to PID, NDOB-MF-iPI, and SFPID controllers, simulation results demonstrate high performances of the proposed $\alpha (t)$ -MF-STFOSMC controller in application for the PEMFC system under different operating conditions.