Real-Time Estimation of PEMFC Parameters Using a Continuous-Discrete Extended Kalman Filter Derived from a Pseudo Two-Dimensional Model

Abstract

Proton Exchange Membrane Fuel Cells (PEMFCs) are clean energy conversion devices that are widely used in various energy applications. In most applications, the main challenge is accurately estimating the state of health (SoH) of the PEMFCs during dynamic operating conditions. Moreover, their behavior is affected by numerous physical phenomena such as heat and membrane flooding. This paper proposes the design of an observer to estimate the PEMFC parameters. A state-space model is first built from 2D physical equations solved by a finite difference in a discretized space domain. The discretized dynamic model is then used to design an observer based on the continuous-discrete extended Kalman filter. The observer has been validated experimentally and is used to estimate the parameters of a PEMFC under dynamic operating conditions. For several load variations, the results obtained using the proposed observer accurately characterize the dynamic responses of PEMFC in real-time.

1. Introduction

1.1. Background

The need for clean energy production across a number of application sectors is expected to increase significantly in the years to come. Fuel cells (FCs) offer a promising alternative clean-energy conversion technology and are widely used in multiple sectors such as transportation and microgrid systems. More particularly, FCs are also associated with optimizing renewable energy storage. Critically, FCs perform best when supplied with pure hydrogen, which results in zero carbon dioxide emissions [1,2]. However, pure hydrogen production, storage, and distribution represent major obstacles to commercializing fuel cell technology [3,4,5]. Another important challenge associated with FC systems is accurately estimating both their dynamic response in real-time and their state of health (SoH). This paper aims to build an observer to reliably predict these performance indicators.

Based on the kind of electrolytes they employ, there are five types of fuel cells [6]: Proton Exchange Membrane Fuel Cells (PEMFCs), Phosphoric Acid Fuel Cells (PAFCs), Alkaline Fuel Cells (AFCs), Molten Carbonate Fuel Cell (MCFCs), and Solid Oxide Fuel Cells (SOFCs). The PAFC, MCFC, and SOFC operate at high temperatures (220 °C, 650 °C, and up to 1000 °C, respectively), whereas the AFC operating temperature may vary from 50 to 200 °C, offering higher performance than the PEMFC. However, in the case of AFC, the hydrogen and oxygen must be pure for optimal performance [7]. Among the various types of fuel cells, Proton Exchange Membrane Fuel Cells (PEMFCs) offer considerable advantages due to their low operating temperature (generally under 100 °C), high power density, quick start-up capability, and long lifetime [8,9]. This overview of the different types of FC demonstrates the key parameters to consider when evaluating overall performance, notably operating temperature and fuel purity. This work focuses exclusively on PEMFCs.

Building an observer for PEMFCs requires understanding the different physical phenomena within this type of FC, such as voltage losses and voltage/current undershoot. Voltage drop in PEMFCs is caused by activation, mass transport, and ohmic losses, which can be characterized by the Tafel slope, mass transport, and ohmic resistances. These three parameters can be estimated experimentally using two dynamic tests: current steps and current sweeps [10]. The PEMFC voltage undershoot phenomenon happens when the current transiently increases to the next step value. This voltage overshoot/undershoot is mainly affected by certain operating conditions such as diverging water saturation levels appearing particularly after prolonged FC operation. The current does not show undershoot due to the effect of double-layer charge [11,12].

1.2. Literature Review

The main challenge when considering the observer design is the complexity of modeling PEMFC systems given that they are non-linear and multi-variable [13,14], and the fact that several phenomena such as heat and membrane flooding can impact PEMFC performance [15]. Numerous models have been presented in the literature, however, these models are either based on the steady-state representation of FCs or are limited to account for specific dynamic behaviors [16,17,18]. These models are either based on sophisticated mathematical equations [19] or data collected from real-life PEMFC [12]. The authors in [12] developed a semi-empirical, two-dimensional, and non-isothermal dynamic model to simulate the transient response in the step changes of the cell potential. The data used to build the model in [12] was collected from performing a specific number of tests.

Diagnosis methods have been also widely used to monitor the degradations of the fuel cell performance [20,21]. Electrochemical impedance spectroscopy (EIS) was used in [22,23] to identify the Warburg impedance and the constant phase element (CPE), which help to observe the hydration of the electrode-membrane state assembly and SoH. However, EIS is not a very convenient solution for onboard integration because the measurement time is long and the EIS device is expensive. Other papers offered advanced diagnostic techniques, such as particle filter-based prognostics [24,25], but these techniques are based on complicated calculation algorithms which limit their application for real-time estimation of PEMFC parameters. Other observers proposed in [26,27] are based on the measurements of non-electrical quantities such as pressure or oxygen concentration. Spatial discretization was applied on the 2D or 3D models of FCs to solve the models using these measurements.

1.3. Study Aim

Developing a PEMFC observer which is implementable in real-time is extremely important in industrial applications. The model on which this observer is based must take into account the most important phenomena and parameters affecting the FC behavior, such as the double-layer and voltage potentials (concentration, activation, and ohmic polarization). In this paper, a simplified 2D physical model that describes these phenomena with physical equations is derived from the first and second Fick’s law and Tafel’s law. To solve these equations, a spatial discretization method is suggested. The space-discretized dynamic model for gas channels, the Gas Diffusion Layer (GDL), and Catalyst Layer (CL) is then obtained. After that, an observer based on the continuous-discrete extended Kalman filter (EKF) is built using the obtained discretized dynamic model. This observer is based on measuring electrical quantities of voltage and current.

Based on the presented literature review, the observer developed in this paper offers certain advantages over previously proposed observers in that it only uses current and voltage measurements, it does not require complicated measurements, it is able to estimate the dynamic response of the PEMFCs, and it can be implemented in real-time.

1.4. Nomenclature and Constants

All the cell geometry parameters and variables are described in

Table 1. Parameters and variables of PEMFC.

Geometry and Model Parameters
Parameter	Unit	Signification	${\mathrm{V}}_{\mathrm{d}}$	m/s	Constant
$\mathrm{b}$	V	Tafel slope	$\mathrm{x}$	M	Spatial coordinate
${\mathrm{C}}_{\mathrm{c}}$	mol/m3	Oxygen concentration in the channel	$\mathrm{y}$	M	Spatial coordinate
${\mathrm{C}}_{\mathrm{g}}$	mol/m3	Oxygen concentration in the GDL	${\mathrm{x}}_{\mathrm{n}},{ \mathrm{y}}_{\mathrm{n}}$	-	Normalized spatial coordinates
${\mathrm{C}}^{\mathrm{ref}}$	mol/m3	Reference concentration on the channel inlet	$\mathsf{\eta }$	V	Cathode potential
${\mathrm{C}}_{\mathrm{cn}}$ ${\mathrm{C}}_{\mathrm{gn}}$	-	Normalized oxygen concentrations	${\mathsf{\eta }}_{\mathrm{n}}$	-	Normalized potential
${\mathrm{C}}_{\mathrm{DL}}$	F/m2	Double-layer capacity	Discretization parameters
${\mathrm{D}}^{\mathrm{eff}}$	m2/s	Effective diffusivity of the GDL	$\mathrm{d}$	-	$\mathrm{Variable} \mathrm{related} \mathrm{to} \mathrm{x}=\frac{d}{D}{h}_d$
$\mathrm{e}$	V	Cell voltage	$\mathrm{D}$	-	Number of grid points in the x-axis
${\mathrm{e}}^{\mathrm{ref}}$	V	Fuel cell potential reference	$\mathrm{k}$	-	$\mathrm{Variable} \mathrm{related} \mathrm{to} \mathrm{y}=\frac{k}{K}{L}_c$
$\mathrm{F}$	C/mol	Faraday constant	$\mathrm{K}$	-	Number of grid points in the y-axis
${\mathrm{h}}_{\mathrm{c}}^{\mathrm{eff}}$	m	Channel depth	n	-	Number of state variables
${\mathrm{h}}_{\mathrm{d}}$	m	GDL thickness	Observer parameters
${\mathrm{i}}_{\mathrm{c}}$	A/m2	Exchange current density	$\mathrm{a}$	-	Constant
$\mathrm{i}$	A	Fuel cell current	$\mathrm{H}$	-	Observation matrix
$\mathrm{j}$	A/m2	Local current density	${\mathrm{K}}_{\mathrm{N}}$	-	Kalman gain
${\mathrm{j}}_{\mathrm{n}}$	-	Normalized Local current density	$\mathrm{F}$	-	Dynamic matrix
${\mathrm{L}}_{\mathrm{c}}$	m	Channel length	$\mathrm{N}$	-	Time index
${\mathrm{l}}_{\mathrm{c}}$	m	Channel width	$\mathrm{q}$	-	White noise covariance
${\mathrm{j}}_{\lim }$	A/m2	Limiting current density	$\mathrm{P}$	-	Covariance matrix
${\mathrm{r}}_{\mathsf{\Omega }}$	Ω m2	Ohmic resistance	${\mathrm{T}}_{\mathrm{s}}$	s	Sampling period
$\mathrm{t}$	s	Temporal coordinate	$\mathrm{R}$	-	Measurement noise covariance
$\mathrm{T}$	°C	Temperature	$\mathrm{v}$	-	Process noise
${\mathrm{u}}_{\mathrm{c}}$	m/s	Air velocity	$\mathrm{x}$	-	State variables

.

2. Physical Two-Dimensional Model

2.1. Geometry and Assumptions

A half-cell including a straight channel, a porous gas diffusion layer (GDL), and a catalyst layer (CL) is considered. A cross-section of this cell is given in Figure 1a. Oxygen transport is considered through the electrode and GDL. Only cathode electrochemical transfers are modeled. The anode is not considered in the model due to much faster kinetic reaction compared to those occurring on the cathode side. Liquid water has not been considered in the model. Figure 1b shows a simple fuel cell electrical circuit.

The model is based on the following assumptions:

• Oxygen transport in the cathode catalyst layer (CCL) is fast;
• The oxygen transport in the GDL is purely 1D in the x-direction;
• The oxygen transport in the channel is assumed to be a plug-flow with an averaged velocity and the choice of plug flow in fuel cell microchannels is justified by the small dimensions of the channels. It was shown in previous studies [29,30] that the laminar velocity profile in rectangular microchannels is flat at the center. Thus, at the first order and given the high mass diffusivity of air in the cathode channel, one can assume that the average velocity of a plug flow is able to give a good description of the mass transfer in the fuel cell microchannels;
• The model is isothermal;
• The CL is supposed to be infinitely small;
• The voltage potential is equal along the CL in the y-dimension;
• All channels of the cell are supposed to be operated in the same condition such as temperature, humidity, etc.; and
• The concentrations in the canal and the GDL are supposed to be equal ${\mathrm{C}}^{\mathrm{ref}}$ (concentration on the channel inlet) before any current is applied to the FC.

2.2. General Equations

The first equation is the one describing the physical model in the channel and the oxygen transport. This equation is represented by the local continuity equation [31]:

$$\frac{\partial C_c}{\partial t}\left(y,t\right)+u_c \frac{\partial C_c}{\partial y}\left(y,t\right)=\frac{D^{eff}}{{h_c}^{eff}} \frac{\partial C_g}{\partial x}\left(0,y,t\right)$$

(1a)

with

$$C_c\left(y=0,t\right)=C^{ref}\left(t\right)$$

(1b)

Equation (1b) represents the boundary conditions at the channel entry. In the GDL, the oxygen transport is described by the first and second Fick’s law as in [31]:

$$\frac{\partial C_g}{\partial t}\left(x,y,t\right)=D^{eff} \frac{{\partial }^2{C}_g}{\partial x^2}\left(x,y,t\right)$$

(2a)

with

$$C_g\left(x=0,y,t\right)=C_c\left(y,t\right)$$

(2b)

The instantaneous current density $\mathrm{j}$ produced by a fuel cell can be written as the sum of a faradaic ${\mathrm{j}}_{\mathrm{f}}$(derived from Tafel law) and a capacitive current density ${\mathrm{j}}_{\mathrm{DL}}$ [31]:

$$j\left(y,t\right)=j_f\left(y,t\right)+j_{DL}\left(y,t\right)$$

(3a)

with

$$j_f\left(y,t\right)=i_c e^{\frac{\eta \left(y,t\right)}{b}} \frac{C_g\left(h_d,y,t\right)}{C^{ref}}, j_{DL}\left(y,t\right)=C_{DL }\frac{\partial \eta }{\partial t}\left(y,t\right)$$

The boundary conditions at the GDL/CL interface are described by

$$\frac{\partial C_g}{\partial t}\left(h_d,y,t\right)=-\frac{1}{4 F D^{eff}} j\left(y,t\right)$$

(3b)

The fuel cell voltage satisfies [31]:

$$e\left(t\right)=e^{ref}-\eta \left(y,t\right)-r_{\Omega }j\left(y,t\right)$$

(4)

The fuel cell current can be computed using

$$i\left(t\right)=l_c{{\displaystyle \int }}_0^{L_c}j\left(y,t\right) dy$$

(5)

In this section, a 2D physical model of the fuel cell has been presented. In the next section, the system of equations is transformed into dimensionless coordinates. The model is then spatially discretized in Section 3 using a finite difference technique. To validate the proposed model, simulations are performed using the Matlab/Simulink environment and compared to experimental data, as detailed in Section 5.

2.3. Normalized Adimensional Equations

For the transformation to dimensionless coordinates, the following normalized variables are used [28]:

$$\begin{gathered} j_n\left(y,t\right)=\frac{j\left(y,t\right)}{j_{lim}}, {\eta }_n\left(y,t\right)=\frac{\eta \left(y,t\right)}{b}, C_{cn}\left(y,t\right)=\frac{C_c\left(y,t\right)}{C^{ref}}, \\ {C}_{gn}\left(x,y,t\right)=\frac{C_g\left(x,y,t\right)}{C^{ref}}, x_n=\frac{x}{h_d}, y_n=\frac{y}{L_c} \end{gathered}$$

(6)

Appendix A shows the value of ${\mathrm{j}}_{\lim }$ and ${\mathrm{V}}_{\mathrm{d}}$.

The physical model of Equations (1)–(5) are rewritten in dimensionless coordinates. In the channel, the oxygen transport equation and the boundary conditions become

$$\frac{\partial C_{cn}}{\partial t}\left(y_n,t\right)+\frac{u_c}{L_c} \frac{\partial C_{cn}}{\partial y_n}\left(y_n,t\right)=\frac{D^{eff}}{h_d {h_c}^{eff}} \frac{\partial C_{gn}}{\partial x_n}\left(0,y_n,t\right)$$

(7a)

$$C_{cn}\left(y_n=0,t\right)=1$$

(7b)

In the GDL, the oxygen transport equation is rewritten as

$$\frac{\partial C_{gn}}{\partial t}\left(x_n,y_n,t\right)=\frac{D^{eff}}{{h_d}^2}\frac{{\partial }^2{C}_{gn}}{\partial {x_n}^2}\left(x_n,y_n,t\right)$$

(8a)

$$C_{gn}\left(0,y_n,t\right)=C_{cn}\left(y_n,t\right)$$

(8b)

In the CL, the normalized current density represented by Equation (3a) becomes

$$j_n\left(y_n,t\right)=\frac{1}{4 F C^{ref}{V}_d}\left(i_c e^{{\eta }_n\left(y_n,t\right)} C_{gn}\left(1,y_n,t\right)+C_{DL} b \frac{\partial {\eta }_n}{\partial t}\left(y_n,t\right)\right)$$

(9)

The fuel cell voltage can be represented as

$$e\left(t\right)=e^{ref}-b {\eta }_n\left(y_n,t\right)-4 F C^{ref }{V}_{d }{r}_{\Omega } j_n\left(y_n,t\right)$$

(10)

The boundary conditions equation at CL/GDL can be rewritten as

$$\frac{\partial C_{gn}}{\partial t}\left(1,y_n,t\right)=\frac{h_d{V}_d}{D^{eff}} j_n\left(y_n,t\right)$$

(11)

The cell current density is calculated by

$$j_n\left(y_n,t\right)=\frac{e^{ref}-e\left(t\right)-b\cdot {\eta }_n\left(y_n,t\right)}{4 F C^{ref} V_d r_{\Omega }}$$

(12a)

$$\frac{\partial {\eta }_n}{\partial t}\left(y_n,\mathrm{t}\right)=\frac{1}{b C_{Dl}} \left(\begin{array}{c}\frac{1}{r_{\Omega }} \left(e^{ref}-e\left(t\right)-b {\eta }_n\left(y_n,\mathrm{t}\right)\right)\\ {-\mathrm{i}}_c e^{{\eta }_n\left(y_n,\mathrm{t}\right)} C_{gn}\left(1,y_n,t\right)\end{array}\right)$$

(12b)

The current density in Equation (12a) is replaced by that in Equation (9) to give the potential Equation (12b).

3. Finite-Difference Discretization

To numerically solve the mathematical model represented by Equations (7)–(12) with known boundary conditions, a finite-difference discretization technique is used in this paper. For the channel and the GDL, a spatial meshing${\mathrm{C}}_{\mathrm{cn}}\left(\mathrm{k}\mathsf{\Delta }\mathrm{y},\mathrm{t}\right),$ ${\mathrm{C}}_{\mathrm{gn}}\left(\mathrm{d}\mathsf{\Delta }\mathrm{x},\mathrm{k}\mathsf{\Delta }\mathrm{y},\mathrm{t}\right)$with discretization steps $\mathsf{\Delta }\mathrm{y}=1/\mathrm{K}$, $\mathsf{\Delta }\mathrm{x}=1/\mathrm{D}$ is considered. The discretization uses central differences when possible and forward or backward differences when not. It is formulated as:

$$\frac{\partial C_{cn}}{\partial y_n}\left(\frac{k}{K},t\right)\approx \left\{\begin{array}{c}\frac{K}{2}\left(C_{cn}\left(\frac{k+1}{K},t\right)-C_{cn}\left(\frac{k-1}{K},t\right)\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}for\hspace{0.33em}k=1\dots K-1\\ K\left(C_{cn}\left(\frac{K}{K},t\right)-C_{cn}\left(\frac{K-1}{K},t\right)\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em} for\hspace{0.33em}k=K\end{array}\right.$$

(13a)

$$\frac{\partial C_{gn}}{\partial x_n}\left(\frac{0}{D},\frac{k}{K},t\right)\approx D\left(C_{gn}\left(\frac{1}{D},\frac{k}{K},t\right)-C_{gn}\left(\frac{0}{D},\frac{k}{K},t\right)\right)$$

(13b)

$$\begin{gathered} \frac{{\partial }^2{C}_{gn}}{{\partial x_n}^2}\left(\frac{d}{D},\frac{k}{K},t\right)\approx D^2\left(\begin{array}{c}{C}_{gn}\left(\frac{d+1}{D},\frac{k}{K},t\right)\\ -2 C_{gn}\left(\frac{d}{D},\frac{k}{K},t\right)+C_{gn}\left(\frac{d-1}{D},\frac{k}{K},t\right)\end{array}\right)\hspace{0.33em}\hspace{0.33em} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for\hspace{0.33em}d=1\dots D-1\hspace{0.33em}and\hspace{0.33em}k=0\dots K \end{gathered}$$

(13c)

$$C_{gn}\left(\frac{D}{D},\frac{k}{K},t\right)\approx C_{gn}\left(\frac{D-1}{D},\frac{k}{K},t\right)-\frac{1}{D}{j}_n\left(\frac{k}{K},t\right)$$

(13d)

with

$$\frac{\partial C_{gn}}{\partial x_n}\left(\frac{0}{D},\frac{k}{K},t\right)\approx C_{cn}\left(\frac{k}{K},t\right)$$

(13e)

$$j_n\left(\frac{k}{K},\mathrm{t}\right)\approx \frac{e^{ref}-e\left(t\right)-b {\eta }_n\left(\frac{k}{K},\mathrm{t}\right)}{4 F C^{ref} V_d r_{\Omega }}$$

(13f)

Equations (7) and (8) can be rewritten using Equation (13a–d).

In the canal, the equations can be written as:

$$\begin{gathered} for\hspace{0.33em}k=1\dots K-1 \\ \frac{\partial C_{cn}}{\partial t}\left(\frac{k}{K},t\right)=-\frac{u_c K}{2 L_c} \left(C_{cn}\left(\frac{k+1}{K},t\right)-C_{cn}\left(\frac{k-1}{K},t\right)\right) \\ +\frac{D^{eff} D}{h_d {h_c}^{eff}} \left(C_{gn}\left(\frac{1}{D},\frac{k}{K},t\right)-C_{cn}\left(\frac{k}{K},t\right)\right) \end{gathered}$$

(14a)

$$\begin{gathered} for\hspace{0.33em}k=K \\ \frac{\partial C_{cn}}{\partial t}\left(\frac{k}{K},t\right)=-\frac{u_c K}{L_c} \left(C_{cn}\left(\frac{K}{K},t\right)-C_{cn}\left(\frac{K-1}{K},t\right)\right)+\frac{D^{eff} D}{h_d {h_c}^{eff}} \left(C_{gn}\left(\frac{1}{D},\frac{K}{K},t\right)-C_{cn}\left(\frac{K}{K},t\right)\right) \end{gathered}$$

(14b)

$\mathrm{For}\hspace{0.33em}\mathrm{k}$ = 0, Equation (1b) is used.

In the GDL, we can write:

$$\begin{gathered} for\hspace{0.33em}d=1\dots D-1\hspace{0.33em}and\hspace{0.33em}k=0\dots K \\ \frac{\partial C_{gn}}{\partial t}\left(\frac{\mathrm{d}}{D},\frac{K}{K},t\right)\approx \frac{D^{eff}\cdot D^2}{{h_d}^2} \left(\begin{array}{c}{C}_{gn}\left(\frac{d+1}{D},\frac{k}{K},t\right)-2 C_{gn}\left(\frac{d}{D},\frac{k}{K},t\right)\\ +C_{gn}\left(\frac{d-1}{D},\frac{k}{K},t\right)\end{array}\right) \end{gathered}$$

(14c)

$\mathrm{For}\hspace{0.33em}\mathrm{d}=\mathrm{D}\hspace{0.33em}\mathrm{and}\hspace{0.33em}\mathrm{k}=0\dots \mathrm{K},$ Equation (13d) is used.

In the Catalyst Layer (CL), the potential can be computed as

$$\frac{\partial {\eta }_n}{\partial t}\left(\frac{k}{K},\mathrm{t}\right)=\frac{1}{b C_{Dl}} \left(\begin{array}{c}\frac{1}{r_{\Omega }} \left(e^{ref}-e\left(t\right)-b {\eta }_n\left(\frac{k}{K},\mathrm{t}\right)\right)\\ {-\mathrm{i}}_c e^{{\eta }_n\left(\frac{k}{K},\mathrm{t}\right)} C_{gn}\left(\frac{D}{D},\frac{k}{K},t\right)\end{array}\right)$$

(14d)

Finally, for $\mathrm{k}=0\dots \mathrm{K}$, the current density is calculated from Equation (13f). The fuel cell current represented in Equation (5) can be calculated by

$$i\left(t\right) \approx \frac{j_{lim} L_c l_c}{K}\left(\frac{1}{2} j_n\left(\frac{0}{K},\mathrm{t}\right)+{\displaystyle \sum }_{k=1}^{K-1}{j}_n\left(\frac{k}{K},\mathrm{t}\right)+\frac{1}{2} j_n\left(\frac{K}{K},\mathrm{t}\right)\right)$$

(15)

4. Description of the Experimental Test Bench and Considered Fuel Cell

A 12 cm² single cell was specially designed to validate the proposed model. It includes two parallel channels of 1.5$\mathrm{mm}$ width and $200 \mathrm{mm}$ in length. The flow field was machined directly in the copper current collector to enable an optical access inside the channel via a transparent Plexiglas plate on the top of the channel (see Figure 2). The current collectors were also gold-plated to ensure minimal ohmic resistance. Between the current collector, a membrane electrode assembly (MEA) was inserted; it was made of a $15 \text{µ}\mathrm{m}$ thick Gore PEM with $0.5 \mathrm{mg}/{\mathrm{cm}}^2$ platinum loading in the CL. The thin membrane and high platinum loading were specially chosen to ensure low membrane resistance and high current density. The MEA was sandwiched between two 10BC GDL from SGL^® compressed to a width of $250 \mu \mathrm{m}$ using two rigid spacers. Finally, two aluminum end-plates were used to assemble all the fuel cell components. A water-cooling circuit connected to a temperature control system was drilled inside the endplates to keep the fuel cell temperature constant.

The fuel cell was connected to a FCT 50 test station from BioLogic^® to control the air and hydrogen flow rates, temperatures, humidity, and pressure. The current and voltage were also recorded through the FCT 50 test station.

5. Experimental Validation of PEMFC Model

Reliable voltage/current curves require a stable environment where the temperature, pressure, humidity, and flow rates are kept constant while the tests are being conducted. If the conditions are fluctuating, the voltage/current characteristics may change. In addition to keeping the testing environment stable, a second consideration is the condition of the fuel cell itself which can take time to stabilize. Depending upon the fuel cell design and size, this stabilization phase may take up to 30 min following a change in current or voltage. All measurements were conducted using the parameters listed in

Table 2. PEMFC model parameter values.

Parameter	Value	Parameter	Value
$\mathrm{b}$	$0.03$ V	${\mathrm{h}}_{\mathrm{c}}^{\mathrm{eff}}$	$1.5\times {10}^{-3}$ m
${\mathrm{C}}_{\mathrm{DL}}$	$200$ F/m2	${\mathrm{h}}_{\mathrm{d}}$	$250\times {10}^{-6}$ m
${\mathrm{C}}^{\mathrm{ref}}$	$7.42$ mol/m3	${\mathrm{i}}_{\mathrm{c}}$	$300$ A/m2
${\mathrm{D}}^{\mathrm{eff}}$	${10}^{-5}$ m2/s	${\mathrm{r}}_{\mathsf{\Omega }}$	$2.2\times {10}^{-5}$ Ωm2
${\mathrm{e}}^{\mathrm{ref}}$	$0.95$ V	$\mathrm{T}$	$70$ °C
$\mathrm{F}$	96,487 C/mol	${\mathrm{u}}_{\mathrm{c}}$	$0.42$ m/s

and under a constant stoichiometric ratio of H₂ (1.1) and air (5), with a relative humidity of 15%.

Using Equations (14) and (15), a PEMFC model was built with the MATLAB/SIMULINK environment. To compare the experimental and simulated results, several tests were carried out. The specifications of the PC used for these calculations are Intel Core i7-3770 CPU and 3.4 GHz with 16 GB RAM.

5.1. Current Profile with Step-Up

To validate our PEMFC model around the nominal operating point ($0.7 \mathrm{V}$), a piecewise constant current was applied. Figure 3 shows the dynamic response of voltage with different current levels of $9$, $10,$ and $11 \mathrm{A}$.

The fuel cell measurements presented in Figure 3 were obtained after an extended period of operation. The voltage variation shown here indicates that when the current steps up to a higher value, the voltage can undergo three changes: an instantaneous drop due to the set of ohmic resistances (electrodes, membrane, bipolar plates, etc.) [33], a quick recovery because the cell takes several seconds for the back-diffusion water to re-wet the anode side so as to ease water redistribution and recover the cell performance [11], and a slow increase due to changes in membrane resistance. Figure 4 compares the simulation and experimental results for the fuel cell current. The relative current error oscillates around zero with a mean value of 1.7%, as shown in Figure 4b.

5.2. Current Profile with Forward/Backward Sweeps

The fuel cell was swept with a current rate of $9.9 \mathrm{mA}/\mathrm{s}$, from $0$ to $20 \mathrm{A}$ (forward sweep) and back to $0 \mathrm{A}$ (back sweep). The measurements of the current and voltage are shown in Figure 5. The voltage gradually decreased from the open-circuit voltage to a low value.

The polarization curve is useful for troubleshooting issues in the fuel cell stack. If the polarization curve is recorded while increasing and decreasing the current, it may show hysteresis. This typically indicates a change in fuel cell conditions, such as drying or flooding of the membrane. Figure 6a represents the voltage of the fuel cell as a function of the current. It was observed that the voltage under the forward sweep is lower than that under the backward sweep. Figure 6b shows the output cell power vs. the current, where the power increases with the increase in current. The power reaches the maximum power point for a current of 18 A and then it decreases [34].

This improved performance during the backward sweep compared with the forward sweep can be explained by the variation of hydration levels in the membrane because the humidity level in the membrane directly influences its conductivity. When the current increases, more water vapor is produced, causing an increase of the internal conductivity. During backward swapping, the membrane was still saturated or at high humidity levels, thus it exhibited an improved performance. In addition, the difference between the forward and backward voltage at a lower current range is more important than at the higher current range. This is due to the increased water content observed in our electrolyte results caused by the increased electrode kinetics.

The PEMFC model was implemented in the Matlab/Simulink environment and was designed with voltage input and current output. Equation (14a–d) were solved with the ODE45 Dormand-Prince variable-step differential equation solver, which is considered to be one of the best methods for treating non-stiff general-form initial value problems for ODE. The current density of the fuel cell was then computed using Equation (13f). The total fuel cell current was obtained from Equation (15). The simulated and experimental currents are shown in Figure 7a. The simulation was conducted with different discretization steps, where various grid sizes (K, D) were used as this was expected to have a negligible impact on the accuracy. Sweeping the fuel cells at different current levels highlighted several phenomena. When the current exceeds 15 A, the voltage is distorted (Figure 5) and the error between the simulation and experimental currents becomes more important (Figure 7b). For the backward sweep, the relative current error increases linearly.

6. Estimating the PEMFC Parameters Using an Extended Kalman Filter

The voltage and current are considered as the input–output model. To design the observer, we used the continuous-discrete extended Kalman filter (CD-EKF) [35]. This technique uses a numerical integration of the ordinary differential equations to predict the state of the continuous-time dynamic system and the corresponding prediction error covariance matrix. A standard discrete-time correction step of the extended Kalman filter (EKF) is then applied to arrive at a more accurate state estimation. The state variables vector was chosen as:

$$x={\left(\dots C_{cn}\left(\frac{k}{K},t\right)\dots C_{gn}\left(\frac{d}{D},\frac{k}{K},t\right)\dots {\eta }_n\left(\frac{k}{K},t\right)\dots C_{DL}\hspace{0.33em}\hspace{0.33em}{r}_{\Omega }\hspace{0.33em}\hspace{0.33em}b\hspace{0.33em}\hspace{0.33em}{i}_c\hspace{0.33em}\hspace{0.33em}{D}^{eff}\right)}^T$$

(16)

where $\mathrm{k}$ goes from 1 to $\mathrm{K}$ for ${\mathrm{C}}_{\mathrm{cn}}$ and from 0 to $\mathrm{K}$ for ${\mathrm{C}}_{\mathrm{gn}}$, and${ \mathsf{\eta }}_{\mathrm{n}}$. $\mathrm{d}$ goes from $1$ to $\mathrm{D}-1$ (for ${\mathrm{C}}_{\mathrm{gn}})$. For $\mathrm{k}$=$0,$ ${\mathrm{C}}_{\mathrm{cn}}$ is determined by the boundary condition (Equation (1b)). For $\mathrm{d}$ = $\mathrm{D},$ ${\mathrm{C}}_{\mathrm{gn}}$ is calculated from the boundary condition (Equation (13d)). Thus, we have $\left(\mathrm{n}=\mathrm{KD}+\mathrm{D}+\mathrm{K}+5\right)$ state variables. The state equations for the concentration in the channel (${\mathrm{C}}_{\mathrm{cn}}$) and GDL (${\mathrm{C}}_{\mathrm{gn}}$), and potential (${ \mathsf{\eta }}_{\mathrm{n}}$) were described previously (Equation (14)). The time derivatives of ${\mathrm{C}}_{\mathrm{DL}}$, ${\mathrm{r}}_{\mathsf{\Omega }}$, $\mathrm{b}$, ${\mathrm{i}}_{\mathrm{c}}$, and ${\mathrm{D}}_{\mathrm{eff}}$ are considered equal to zero. The system of interest is a continuous-time dynamic system with a discrete-time measurement given by [36,37,38]:

$$\dot{x}=f\left(x\left(t\right),\hspace{0.33em}e\left(t\right)\right)+\nu \left(t\right)$$

(17a)

$$i\left(N T_s\right)=h\left(x_N,e\left(N T_s\right)\right)+w_N$$

(17b)

where $\mathrm{v}\left(\mathrm{t}\right)$ is the unmeasured “process noise” which is assumed to be a continuous-time zero-mean white noise of covariance $\mathrm{q}\left(\mathrm{t}\right)$ and ${\mathrm{R}}_{\mathrm{N}}$ is measurement noise covariance. ${\mathrm{w}}_{\mathrm{N}}$ is the measurement noise which is assumed to be discrete-time white noise with zero mean, ${\mathrm{x}}_{\mathrm{N}}=\mathrm{x}\left({\mathrm{N} \mathrm{T}}_{\mathrm{s}}\right)$ with ${\mathrm{t}}_{\mathrm{N}}={\mathrm{N} \mathrm{T}}_{\mathrm{s}}.$ The Kalman Filter mainly involves two steps: prediction and measurement update. The predicted state ${\hat{x}}_{\mathrm{p}}$ and its covariance $\mathrm{P}$ are calculated by solving Ordinary Differential Equations (ODE) [39]:

$${\dot{\hat{x}}}_p\left(t\right)=f\left({\hat{x}}_p\left(t\right),\hspace{0.33em}e\left(t\right)\right)$$

(18a)

$$\dot{P}\left({\hat{x}}_p\left(t\right)\right)=F\left({\hat{x}}_p\left(t\right)\right) P\left(t\right)+P\left(t\right) F^T\left({\hat{x}}_p\left(t\right)\right)+q\left(t\right)$$

(18b)

The dynamic matrix $\mathrm{F}$ (the Jacobian matrix of the partial derivatives of function $\mathrm{f}$) is computed as follows:

$$F\left({\hat{x}}_p\left(t\right)\right)=\frac{\partial f}{\partial x}\left({\hat{x}}_p\left(t\right),e\left(t\right) \right)=\left(\begin{array}{c}Par{t}_I\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}Par{t}_{II}\\ Par{t}_{III}\hspace{0.33em}Par{t}_{IV}\end{array}\right)$$

(19)

The matrix parts (Equation (19)) are given in Appendix B. The equations used to calculate the different terms in the Jacobian matrix are also given in Appendix C.

The two parts of Equation (18) can be solved simultaneously with an ODE solver [40]. Equation (18b) is vectorized. Here, the filter (CD-EKF) is simulated in Matlab/Simulink with ODE45. The integration is done in the interval $\left[{\mathrm{t}}_{\mathrm{N}-1},{ \mathrm{t}}_{\mathrm{N}}\right]$ with the initial conditions $\hat{\mathrm{x}}({\mathrm{t}}_{\mathrm{N}-1})={\hat{\mathrm{x}}}_{\mathrm{N}-1/\mathrm{N}-1 }$ and error covariance matrix $\mathrm{P}({\mathrm{t}}_{\mathrm{N}-1})={\mathrm{P}}_{\mathrm{N}-1/\mathrm{N}-1}$. The predicted state ${\hat{\mathrm{x}}}_{\mathrm{N}/\mathrm{N}-1 }={\hat{\mathrm{x}}}_p({\mathrm{t}}_{\mathrm{N}})$ and the error covariance matrix ${\mathrm{P}}_{\mathrm{N}/\mathrm{N}-1 }=\mathrm{P}({\mathrm{t}}_{\mathrm{N}})$ at the time point ${\mathrm{t}}_{\mathrm{N}}$ are computed. Then, the standard correction step of the EKF is used. Additionally, only the upper triangular part of the error covariance matrix, including the main diagonal because of its symmetry, is calculated. The standard measurement update is applied at the time (${\mathrm{t}}_{\mathrm{N}}={\mathrm{N} \mathrm{T}}_{\mathrm{s}}$). The observation function $\mathrm{h}$ can be deduced from Equations (15) and (13f) as:

$$i\left(t\right)=\frac{L_c\cdot l_c}{r_{\Omega }}\cdot \left(e^{ref}-e\left(t\right)\right)-b\cdot \frac{L_c\cdot l_c}{r_{\Omega }\cdot K}\cdot \left(\begin{array}{c}\frac{1}{2}{\eta }_n\left(\frac{0}{K},\mathrm{t}\right)\\ +{\displaystyle {\displaystyle \sum }_1^{K-1}}{\eta }_n\left(\frac{k}{K},\mathrm{t}\right)+\frac{1}{2}{\eta }_n\left(\frac{K}{K},\mathrm{t}\right)\end{array}\right)$$

(20a)

Supposing the constant $\mathrm{a}=\mathrm{KD}+\mathrm{D}-1$, we can rewrite Equation (20a) as a function of the state variables as

$$i\left(t\right)=\frac{L_c\cdot l_c}{r_{\Omega }\cdot K}\cdot \left(e^{ref}-e\left(t\right)\right)-b\cdot \frac{L_c\cdot l_c}{r_{\Omega }\cdot K}\cdot \left(\begin{array}{c}\frac{1}{2}{x}_a\left(\frac{0}{K},\mathrm{t}\right)\\ +{\displaystyle {\displaystyle \sum }_1^{K-1}}{x}_{\left(a+k\right)}\left(\frac{k}{K},\mathrm{t}\right)+\frac{1}{2}{x}_{\left(a+K\right)}\left(\frac{K}{K},\mathrm{t}\right)\end{array}\right)$$

(20b)

The observation matrix $\mathrm{H}$ can be calculated by

$$H\left(x\right)=\frac{\partial h}{\partial x}\left(x,e\left(N T_s\right)\right)=\left[\mathrm{0..0}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\frac{h_{\eta }}{2}\dots h_{\eta }\dots \frac{h_{\eta }}{2}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{h}_r{}_{\Omega }\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{h}_b\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0\hspace{0.33em}\right]$$

(21)

where

$$\begin{gathered} h_{\eta }\left(x\right)=-\frac{x_{\left(a+K+4\right)}\cdot L_c\cdot l_c}{x_{\left(a+K+2\right)}\cdot K} \\ {h}_r{}_{\Omega }\left(x\right)\hspace{0.33em}=-\frac{1}{x_{\left(a+K+2\right)}}\cdot i\left(t\right) \\ {h}_b\left(x\right)=-\frac{L_c\cdot l_c}{x_{\left(a+K+2\right)}\cdot K}\cdot \left(\begin{array}{c}\frac{1}{2}{x}_a\left(\frac{0}{K},t\right)\\ +{\displaystyle {\displaystyle \sum }_1^{K-1}}{x}_{\left(a+k\right)}\left(\frac{k}{K},t\right)+\frac{1}{2}{x}_{\left(a+K\right)}\left(\frac{K}{K},t\right)\end{array}\right) \end{gathered}$$

The Kalman gain can be computed as:

$$K_N=P_{N/N-1} H^T\left({\hat{x}}_{N/N-1}\right) S_N{}^{-1}$$

(22a)

$$S_N=H\left({\hat{x}}_{N/N-1}\right) P_{N/N-1} H^T\left({\hat{x}}_{N/N-1}\right)+R_N$$

(22b)

As presented in [41], Auger et al. proposed that the value of covariance $\mathrm{R}$ can be set to one (${\mathrm{R}}_{\mathrm{N}}=1$). The state estimate covariance is updated by the equation:

$$P_{N/N}=\left(1-K_N H\left({\hat{x}}_{N/N-1}\right)\right) P_{N/N-1}$$

(23)

We use the Kalman gain to update the state estimate:

$${\hat{x}}_{N/N}={\hat{x}}_{N/N-1}+K_N (i\left(N T_s\right)-h\left({\hat{x}}_{N/N-1},e\left(N T_s\right)\right)$$

(24)

Lastly, the calculation process repeats itself. The proposed EKF requires a high computation effort and is not appropriate for online estimation of PEMFC parameters. Hence, some simplifications have been implemented to calculate the covariance matrix $\mathrm{P}$ of the predicted state $\mathrm{x}$ (Equation (18b)) using Mazzoni’s formulae [42]. Mazzoni’s formulae, which are based on the implicit midpoint method, are rearranged as follows:

$$P_{N/N-1}=M_{1/2} P_{N-1/N-1} M_{1/2}^t+T_s N_{1/2} q N_{1/2}^t$$

(25a)

$$M_{1/2}=N_{1/2} \left(I_n+\frac{T_s}{2} F\left(\hat{x}\left(t_{N-1}+\frac{T_s}{2}\right)\right)\right)$$

(25b)

$$N_{1/2}=\left(I_n-\frac{T_s}{2} F\left(\hat{x}\left(t_{N-1}+\frac{T_s}{2}\right)\right)\right)$$

(25c)

where ${\mathrm{I}}_{\mathrm{n}}$ is the identity matrix of size $\mathrm{n}$. Due to the application of the simplification method, it is expected that the model accuracy will be decreased, albeit the model-solving is faster [43].

7. Results of the EKF Observer and Discussion

7.1. Current Profile with Forward/Backward Sweeps

To ensure that all results obtained in this study are not overly sensitive to the fineness of the spatial discretization, a mesh sensitivity analysis was performed. The values of$\mathrm{K}$ and $\mathrm{D}$ were optimized by simulation. To validate the EKF observer designed for $\mathrm{K}=20$ and $\mathrm{D}=15$, the same voltage/current measurements shown previously in Figure 5 were used. The PEMFC current measured experimentally was compared to the current estimated by the EKF ${\hat{\mathrm{i}}}_{\mathrm{N}/\mathrm{N}}$($\mathrm{h}\left({\hat{\mathrm{x}}}_{\mathrm{N}/\mathrm{N}},\mathrm{e}\left({\mathrm{N} \mathrm{T}}_{\mathrm{s}}\right)\right)$. The current estimated by EKF converged well with the experimental data, as shown in Figure 8a. Figure 8b shows a zoomed-in image of both curves between 2177 s and 2187 s; it confirms the estimated current’s convergence with the experimental current. Except at the beginning and end of the curve (where the PEMFC current is of the order of a few mA), the relative error is quite low, as shown in Figure 8c, with a mean value of $0.002\%.$ As shown in Figure 7a, there is a significant error between the experimental data and the output of the model simulated using time-invariant parameters. This error becomes negligible when using time-varying parameters estimated by a continuous-discrete EKF, as illustrated in Figure 8a.

Ambient air was employed as the oxygen source at the cathode, oxygen was consumed by the reaction at CL, and the concentration of oxygen was reduced between the inlet and outlet of the gas flow channel. When the FC current is higher, the electrochemical reaction becomes stronger and more oxygen is consumed at CL, leading to lower oxygen concentration inside the GDL. The variation of the oxygen concentrations ${\mathrm{C}}_{\mathrm{c}}$ in the channel is presented in Figure 9a as a function of time and $\mathrm{k}$. This concentration strongly depended on the current value (see Figure 8a). As an example, for $\mathrm{k}=20$, this concentration decreased with time until the minimum value was reached at time $\mathrm{t}=2181$ s for a current of 20 A. Diffusion across the cathode GDL is caused by a concentration differential between the gas flow channel and the CL. The variation of concentration ${\mathrm{C}}_{\mathrm{g}}$as a function of time and $\mathrm{d}$ for two values of $\mathrm{k}$ ($0$,$\mathrm{K}=20$) is presented in Figure 9b. Following this graph, it can be seen that the ${\mathrm{C}}_{\mathrm{g}}$ concentration values decrease with time as a result of the current variation, especially for the highest value of $\mathrm{d}$. For $\mathrm{d}=15$, for example, the concentration decreased over time until it reached its minimum value at time $\mathrm{t}=2181$ s for a current of 20 A.

Figure 10 shows the variation of the estimated parameters with time. Figure 10a shows the variation of potential $\mathsf{\eta }$ as a function of time for different values of $\mathrm{k}$. The value of potential $\mathsf{\eta }$ follows the PEMFC current values. The hysteresis phenomena can be observed on the potential curve for the forward/backward sweep. As illustrated in Figure 10b, the capacitances ${\mathrm{C}}_{\mathrm{DL}}$vary slightly. As presented in Figure 10c, the estimated resistance ${\mathrm{r}}_{\Omega }$increased hardly for the forward sweep and for the backward sweep, it was still constant. This could be attributed to the FC’s temperature as well as to hydration. A proper water content in a fuel cell will reduce the membrane resistance and promote the electrochemical reaction in the cell [44]. Tafel slope $\mathrm{b}$ represents the voltage losses due to activation polarization. It depends on current [31], humidity, and operating temperature [45]. The Tafel slope $\mathrm{b}$ decreased lightly, as depicted in Figure 10d. As shown in [31], this variation is fairly small because $\mathrm{b}$ can fluctuate from $20 \mathrm{mV}$ to $350 \mathrm{mV}$ as a function of various operating conditions [31]. The diffusion coefficient ${\mathrm{D}}^{\mathrm{eff}}$ was slightly variable within the backward sweep (Figure 10f). The results explained in [46] show that the through-plane gas diffusion coefficient decreases as the number of hydrophilic pores in the GDL increase due to blocking of pores by liquid water. The exchange current density depends on the roughness of the electrode surface; its catalytic properties; the presence of oxides and water; the concentration of the reactive species; the composition of the electrode and electrolyte; and the temperature [47]. As can be seen in Figure 10e, the exchange current density ${\mathrm{i}}_{\mathrm{c}}$ increased greatly. The exchange current density increases with the increase of temperature [48]. The exchange current density has three phases: a constant between 0 s and 1200 s because the current is low (the temperature is close to the starting temperature); an exponential increase between 1200 s and 2100 s, where the current rises and the internal temperature rises with it; and a second exponential increase between 2100 s and 4000 s.

7.2. Current Profile with Step-Up/Down

The EKF observer was applied with the same initial parameters as in the section above and for the experimental results presented in Figure 3. The Tafel slope $\mathrm{b}$, capacitance ${\mathrm{C}}_{\mathrm{DL}}$, and current ${\mathrm{i}}_{\mathrm{c}}$ were slightly variable for this short simulation period. The potential value $\mathsf{\eta }$, as well as the ${\mathrm{C}}_{\mathrm{c}}$ and ${\mathrm{C}}_{\mathrm{g}}$ concentrations, were heavily influenced by the PEMFC current levels. The resistance ${\mathrm{r}}_{\Omega }$ estimated by EKF is represented in Figure 11.

Figure 11 shows how the sharp voltage drops at t = 300 s and 600 s (see Figure 3) because of inevitable jumps in the membrane resistance. Then, the resistance decreases or increases slowly according to voltage variation. The flooding effect may cause transient fluctuations in membrane resistance, as seen in Figure 11 after a long course of running. Reference [49] explains the relationship between current density, membrane resistance, and water content. A positive feedback loop exists between current density, which produces water, and membrane resistancewhich increases with water contentas the water content increases. However, the membrane resistance is typically a small component of the overall cell resistance. For the feedback mechanism to destabilize the cell and generate hysteresis, the membrane resistance must be particularly sensitive to water content.

7.3. Modified PEMFC Cell with Step-Up/Down Profile

We used the same fuel cell mentioned previously but its channel width was changed from $1.5 \mathrm{mm}$to $0.5 \mathrm{mm}$. The current and voltage profiles with step-up/down were conducted as shown in Figure 12a,b. The current consisted of two step-up cycles (0 to 25 A) interspersed with a hard drop from 25 A to 5 A and then to 1 A, while the voltage varied in a range of 0.3 to 0.9 V (Figure 12b). This current/voltage profile can span the whole operating area of the investigated FC. It enables large changes in the parameters of the FC model. The EKF estimated and experimental currents are compared in Figure 12a. The relative current error is illustrated in Figure 12c. Its value is lower than 0.035% but takes a higher value when the current steps up or steps down.

Figure 13 shows the parameters estimated with EKF vs. time. Figure 13a shows the variation of potential $\mathsf{\eta }$ with time. The current value substantially influences the potential $\mathsf{\eta }$. As shown in Figure 13a, its variation was similar to the current one. As explained in [11], when the current drops from a high value to zero or a low value, the double-layer charge effect can appear. Figure 13b illustrates that despite the large current variation, the value of the double-layer capacitor dod not vary. Due to the low hydration and temperature at a low current, the resistance ${\mathrm{r}}_{\Omega }$ began with a high value and then gradually declined, as seen in Figure 13c. The resistance ${\mathrm{r}}_{\Omega }$decreases with the increase in humidity, current, and temperature [50,51,52]. The Tafel slope $\mathrm{b}$ increased with time, as indicated in Figure 13d. The rise in $\mathrm{b}$ could be attributable to an increase in temperatures and decreasing pore porosity, which decreases as humidity rises. The exchange current density vs. time is depicted in Figure 13e. The exchange current took a higher value compared with previous cases (Figure 10e). As depicted in Figure 13f, the diffusion coefficient ${\mathrm{D}}^{\mathrm{eff}}$increased from $1\times {10}^{-5}$ to $1.35\times {10}^{-5}{ \mathrm{m}}^2/\mathrm{s}$. The increase in ${\mathrm{D}}^{\mathrm{eff}}$ could be due to an increase in temperature. When the current dropped from $25 \mathrm{A}$ to $5 \mathrm{A}$, the diffusion coefficient jumped downward, as illustrated in Figure 13f. The effective diffusion coefficient decreases as the number of hydrophilic pores increases because liquid water blocks the pores in the GDL [46].

Table 3. EKF results.

Period of Real-Time Measurement (s)	Number of Samples	${\mathbf{T}}_{\mathbf{s}} \left(\mathbf{s}\right)$	Period of Simulation Time (s)
3979	23,471	0.1695	2340
6000	81,455	0.073	5188

summarizes the observer performances in terms of processing for the current profiles presented in Section 7.1 and Section 7.3. The EKF simulation with the data considered in Section 7.1 lasted 2340 s, which is substantially less than the duration of the real-time measurement (3779 s). The EKF simulation with data in Section 7.3 lasted 5188 s against 6000 s for the period of the real-time measurement. This means that the EKF can estimate the PEMFC parameters before the time period between two samples elapsed. In addition, these simulations demonstrate that despite the enormous number of state variables (n = 340), the observer was able to estimate the PEMFC parameters faster than the real-time measurement.

8. Conclusions

In this paper, the design of a 2D physics-based diagnostic model for the PEMFC is presented. The model is extracted from Fick’s laws and Tafel’s equation, which relate to the oxygen transport in fuel cell channels and the gas-diffusion layer (GDL). The proposed PEMFC model is spatially discretized and using this discretized PEMFC model, a continuous-discrete extended Kalman filter was designed. Mazzoni’s formulae are later applied on the EKF to accelerate the calculation time of the PEMFC state observer.

For various current profiles, the EKF was used to estimate the PEMFC parameters. This observer accurately estimates the PEMFC’s dynamic ohmic resistance considering the load variation. The most important characteristic of the observer is its capability to estimate the PEMFC parameters in real-time.

Other operating conditions affect the PEMFC, such as aging and the membrane state (dry, wet, and flooded). In addition, the ohmic resistance of the PEMFC depends on the aging of the PEMFC and, consequently, the SoH of the PEMFC affects its ohmic resistance. In future work, the PEMFC observer can be extended to estimate both lifetime and aging.