Nonlinear Model Predictive Control of an Automotive PEMFC System Considering Anode Recirculation Loop Bleeding

1. Introduction

Fuel cell electric vehicles have attracted increasing attention in recent years owing to their low pollution and extended driving range [1]. The primary source of power in this type of electric vehicle is the proton exchange membrane fuel cell (PEMFC). This type of fuel cell offers several benefits when compared with other fuel cell type, including zero emissions, low operating temperatures, high power output, minimal noise, excellent startup performance, and a long lifespan [2,3]. However, there are still some challenges that impede its massive commercialization. Lifetime still remaining one of the important parameters to measure the performance of a PEMFC [4]. It is affected by the degradation of the fuel cell components as the primary cause of short service time. The harsh operating condition of temperature, pressure, stoichiometry, and humidity are the main factors that affect the degradation and the lifetime. Indeed, in automotive applications, the operational conditions are subject to constant changes due to dynamic load cycling. As a result, the operating conditions undergo fluctuations [5] and consequently, there is a possibility for the operating conditions to take harsh values, leading to increased degradation and a deterioration in the system performance. For instance, during fast load increases, there is the possibility of oxygen and hydrogen starvation that mainly causes severe corrosion of the carbon support, uneven distribution, and cell voltage reversal [6].

Water is required to humidify the membrane but overhydrated or flooding causes reactants blockage and gas starvation [7]. Moreover, water and thermal management has to be coordinated. Indeed, non-uniform temperature distribution may cause condensation of water in the reactant gases [8]. If the temperature is maintained at excessively high levels, it causes dehydration and shrinkage of the membrane. As a consequence, the dry membrane exhibits higher Ohmic resistance, leading to a reduction in output voltage [9,10]. Finally, the proper stack temperature is needed to maintain optimum electrochemical reactions and keep the PEMFC stack material’s integrity.

The control of these variables is accomplished through three subsystems: the anode subsystem, the cathode subsystem, and the thermal subsystem. Each subsystem is equipped with suitable local controllers that receive the information of desired set points. Thus, the optimal operating conditions that keep the fuel cell system in a safe working region, avoid starvation and concurrently minimize fuel consumption are needed. These optimal operating conditions also determine the performance, lifetime, fuel utilization, and response times.

A supervisory controller that calculates the optimal operating conditions and contributes to the long-term operation of the fuel cell system is required. The supervisory controller is an algorithm that considers the system constraints satisfying the operating restrictions and determine the desired optimal operating conditions. The optimal operating conditions are converted into setpoints, which are provided to the local controllers of each subsystem.

According to the literature review, there are many studies in operating condition optimization and evaluation [11]–[14]. However, there are still shortcomings in obtaining and controlling the optimal operating conditions of the fuel cell stack in the case of dynamic load applications and specifically, in the automotive application.

Referring the operating conditions determination for the automotive applications, the following works can be mentioned. In [15] a multi-parametric model predictive control strategy was proposed to find out and control the air and hydrogen mas flow rate, temperature and stack current that provide the requested power. However, the air humidification was disregarded. A neural network model predictive control strategy based on nonlinear autoregressive exogenous (NARX) model was proposed in [16]. Despite their proposal effectively determine the optimal operating air humidity ratio, pressure, and stack temperature, it is worth mentioning that they overlooked the consideration of the dead-end anode and the potential consequences of gas species accumulation within the anode loop. A multi-objective optimization analysis was done in [17], but the controller was designed only to calculate the optimal oxygen stoichiometry. A data driven model predictive control strategy was designed for trajectory tracking of the oxygen stoichiometry. The strategies of using the experimental data and particle swarm optimization were used to obtain the optimal operating temperature in [18], however, the optimization was done just for temperature.

Apart from the crucial factors of air pressure, oxygen stoichiometry, humidity ratio, stack temperature, and hydrogen stoichiometry that significantly affect performance and efficiency, it is essential to also take into consideration the influence of anode recirculation loop purge on the overall efficiency. During the operation of a fuel cell system, although the electrolyte is considered to be a gas barrier, air can diffuse from the cathode to the anode through the membrane [19]. While the diffused oxygen reacts immediately and is converted into water, inert nitrogen accumulates in the anode channel. The concentration of nitrogen in the anode volume increases over time, leading to hydrogen starvation [20]. To reduce the nitrogen concentration in the anode module, the closed anode loop needs to be opened.

The purge valve is the solution that is used to purge the anode loop. Opening the purge valve could be periodic or continuous, the latter also called bleeding [21]. The challenge of periodic purging lies in determining the purge time duration and the purge cycling. Researchers have proposed some strategies to address this issue. Voltage-based [22] and nitrogen content-based [23] are two basic strategies used to control the purge.

The use of a constant voltage drop or a constant nitrogen concentration accumulation value as an index for determining purge specifications has certain drawbacks. This is attributed to the dependence of nitrogen crossover on the membrane’s humidification conditions, the stack temperature and the load current. In different operating conditions, the nitrogen crossover is different and consequently, the voltage drops value or nitrogen accumulation value is different. Therefore, different load current needs different voltage drop value as an index. MPC was proposed in [24] to address this limitation and determine the purge duration corresponding to the load current.

The anode recirculation loop equipped with bleed valve for continuous purging, despite experiencing continuous fuel loss through the bleed valve, has an optimal bleed rate. At steady state, the nitrogen that crossover reaches equilibrium with the bleed nitrogen. In the optimal bleed rate, the effect of bleed on the partial pressure of hydrogen and nitrogen and on the efficiency will be constant. Compared with the voltage-based purge strategy, the energy efficiency is higher in the bleed strategy [22]. A purge strategy based on the hydrogen content at the stack outlet was proposed by [21], however, this strategy does not improve efficiency.

Additionally, comparing the periodic purge and the continue purge or bleed, the later one is more stable in terms of voltage ripple. However, it is crucial to achieve optimal bleeding that serves two purposes. On one hand, it ensures proper nitrogen bleeding, preventing nitrogen accumulation in the anode channel, which could lead to fuel starvation. On the other hand, it aims to limit excessive purging, which would result in hydrogen loss and reduced efficiency.

The current research has made improvements compared to the previous work [25]. Specifically, the stack model has been enhanced with more accurate nitrogen crossover modeling and the anode recirculation modeling. Additionally, this study takes into consideration the impact of bleed and water separator in the anode recirculation loop and optimizes the bleed rate accordingly.

The purpose of this study is to analyze and consider the effect of all operating variables that influence the overall efficiency of the fuel cell system. The stack efficiency and corresponding system efficiency, simultaneously depend on the operating parameters such as inlet air pressure, inlet air mass flow rate and humidity ratio, stack temperature, the fuel pressure, and bleed valve opening.

Nonlinear Model Predictive Control (NMPC) as a multivariable control strategy takes into account the actuator limitations, states and outputs constraints to calculate the optimal control variables. In NMPC, a multi-objective optimization problem is tackled within the prediction horizon. This optimization procedure accounts for both the physical and manufacturing constraints imposed on the operating variables such as the compressor map. In order to leverage the capabilities of Model Predictive Control (MPC) in addressing constraints, predicting system behavior within the prediction horizon, and adapting to changing load conditions in automotive applications, this strategy is employed in our work.

The main contributions and sections of this paper are listed as follows:

• NMPC design strategy for the real-time computation of the optimal control variables for the fuel cell system, taking into account the stack current. The considered controlled variables are the inlet air pressure, inlet air mass flow rate, inlet air humidity ratio, stack temperature and bleed rate in the anode recirculation loop.
• Comprehensive consideration and modeling of various components in the anode recirculation loop. Specifically, the dynamic behavior of water separator, bleed valve, ejector, return manifold and supply manifold. This holistic modeling enables a more accurate representation of the system dynamics, allowing for a better understanding and analysis of the anode recirculation loop’s overall performance.
• Use as simulation model of an improved model corresponding to the Inn-balance project [26] fuel cell system. The model includes all sensors and actuators of the real system and it is experimentally validated. Additionally, the control-oriented model implemented in the MPC considers the nitrogen crossover through the membrane and the influence of liquid water in the cathode catalyst layer. This liquid water affects the electrochemically active surface area (ECSA) and, consequently, the stack voltage.

The paper is organized into the following sections. Section 2 describes the structure of the fuel cell system used. The fuel cell system modeling is presented in section 3. In section 4, an offline setpoint generator used as baseline for comparison analysis is explained. In section 5, the proposed nonlinear MPC structure, nonlinear constraints and cost function is presented. The simulation results of the proposed control structure are presented in section 6. The main conclusions and key results are summarized in section 7.

2. Fuel Cell System Structure

The fuel cell system considered in this work is shown in Figure 1. It corresponds to the Inn-Balance project [26]. It has four balance of plant subsystems including: anode, cathode, thermal, and power electronic subsystem. The cathode subsystem supplies the air with the desired mass flow rate and pressure. A membrane humidifier is used to humidify the inlet air to control the humidity of the stack. The thermal subsystem provides the coolant liquid at the necessary temperature to cool down the stack. Furthermore, the thermal subsystem is equipped with a loop including a heater and pump to warm the stack in case of freeze startup. The anode subsystem supplies the hydrogen for the stack with proper stoichiometry. To improve the efficiency of hydrogen utilization and promote a more uniform concentration within the anode channel, the anode subsystem has a recirculation pathway equipped with an ejector-injector. Additionally, a water separator and bleeding valve are used to remove the accumulated water and prevent anode flooding.

3. Modeling the Fuel Cell System

NMPC is a model-based control strategy that utilizes the system model to predict the behavior of the system during the prediction horizon. Based on the system’s behavior, it determines the corresponding control signals. In this section, the model of the fuel cell system which is used by NMPC to predict the system behavior, called the control oriented model (COM) is described.

3.1. Hydrogen Recirculation System

The fuel cell system in this paper includes a hydrogen recirculation comprising the anode supply manifold, anode return manifold, hydrogen tank, pressure control valve, water separator, bleed valve, and a passive ejector-injector for the flow recirculation. Indeed, in order to enhance the hydrogen utilization efficiency and maintain a consistent concentration within the stack, the anode subsystem incorporates an ejector-injector in the recirculation pathway. This passive component allows unreacted hydrogen to recirculate back to the supply manifold. Additionally, the anode subsystem includes a water separator and purge valve to remove accumulated liquid water and nitrogen diffused through the membrane from the cathode to the anode side.

3.1.1. Anode Supply Manifold

The pipeline that connects the ejector output with the anode channel inlet is the supply manifold. The pressure dynamics of the different gas species in the supply manifold can be written as Eq. (1)

$${\displaystyle \frac{{dp}_{sm,i}}{dt}}={\displaystyle \frac{R{T}_{sm}^a}{V_{sm}^a{M}_i}}(W_{i,in}^{sm}-W_{i,out}^{sm})$$

(1)

where $W_{i,in}^{sm}$, $W_{i,out}^{sm}$ $(i:N_2,H_2,H_{2}O)$ are respectively the inlet and outlet mass flow rates of the different species. $T_{sm}^a\mathrm{a}\mathrm{n}\mathrm{d}{V}_{sm}^a$are the supply manifold temperature and volume. The pressure difference between the anode supply manifold and the anode channel is relatively small and the supply manifold outlet flow can be calculated as Eq.(2), where flow is proportional to the pressure difference between the supply manifold and the anode channel [27].

$$W_{out,sm}=K_{sm,out}(p_{sm}-p_{an})$$

(2)

where the $K_{sm,out}$ constant depends on the physical characteristics of the orifice.

3.1.2. Anode Return Manifold

The anode return manifold is the pathway that connects the anode channel outlet to the water separator. It is modelled same as the supply manifold. The dynamic equations of the pressure inside the return manifold are given by Eq.(3)

$${\displaystyle \frac{{dp}_{rm,i}}{dt}}={\displaystyle \frac{R{T}_{rm}^a}{V_{rm}^a{M}_i}}(W_{i,in}^{rm}-W_{i,out}^{rm})$$

(3)

where $W_{i,in}^{rm}$, $W_{i,out}^{rm}$ $(i:N_2,H_2,H_{2}O)$ are respectively the inlet and outlet flow rates of the different species. $T_{rm}^a\mathrm{a}\mathrm{n}\mathrm{d}{V}_{rm}^a$are the return manifold temperature and volume [27]. The total outlet flow is given by the nozzle dynamic Eq.(4)

$$W_{out,rm}=K_{rm,out}(p_{rm}-p_{sep})$$

(4)

where the $K_{rm,out}$ constant depends on the physical characteristics of the orifice.

3.1.3. Ejector

In order to avoid the decrease in hydrogen partial pressure caused by the accumulation of inert nitrogen gas and water in the anode during prolonged operation, a purged recirculation loop is a good solution. This loop allows for continuous gas flow, preventing localized gas shortages and maintaining sufficient hydrogen partial pressure by removing accumulated inert gas and water. By adopting this approach, the overall performance and longevity of the fuel cell stack can be improved. In the recirculation loop, either a pump or an ejector can be used, each having its own advantages and disadvantages [28]. In this study, an ejector acting as an alternative to the hydrogen recirculation pump is used to recirculates the unconsumed hydrogen. The structure of the ejector is shown in Figure 2.

The ejector has two inlet ports. The primary and the secondary inlet ports, that are connected to two sources or manifolds. The primary inlet is connected to a high-pressure source at pressure $p_p$. The flow entering this port is referred to as the primary flow. On the other hand, the secondary inlet is connected to a low-pressure source at a pressure$p_s$. The outlet flow from the cell stack is sucked into the secondary inlet. The primary flow is obtained using the equation for a chocked flow nozzle, as

$$W_p={\displaystyle \frac{P_p{A}_t}{\sqrt{R{T}_P}}}\sqrt{\gamma {\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{\gamma -1}}}}\sqrt{{\eta }_P}$$

(5)

where ${\eta }_P$ is the efficiency of the primary nozzle (typical value 0.95)$. \gamma$ and $T_p$ are the ratio of specific heats of the recirculation gases and the temperature of the primary gas. To calculate the secondary flow, it is assumed that the secondary flow is chocked at section$y-y$, i.e., the Mach number of the secondary flow at section$y-y$, $M_{sy}$, is 1. If the effective area occupied by the secondary flow at section$y-y$, $A_{sy}$ is known, the secondary flow can be determined as Eq.(6) [29]

$$W_s={\displaystyle \frac{P_{sep}{A}_{sy}}{\sqrt{R{T}_s}}}\sqrt{\gamma {\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{\gamma -1}}}}\sqrt{{\eta }_s}$$

(6)

where ${\eta }_s$ is the efficiency for the secondary flow path (typical value 0.85). $P_{set}$ and $T_s$ are the separator pressure and separator temperature. However, $A_{sy}$ is a hypothetical area and cannot be directly measured. Since the total area occupied by the primary and the secondary flows equals the mixing tube area, the effective area occupied by the secondary flow can be obtained using

$$A_{sy}=A_{mix}-A_{py}$$

(7)

where $A_{py}$ is the area occupied by the primary flow. The parameters $A_{py}$and Mach number of primary flow at section $y-y$, $M_{py}$, are obtained using the compressible flow equations Eqs.(8), (9)

$$A_{py}=A_t{\displaystyle \frac{1}{M_{py}}}{\left({\displaystyle \frac{2}{\gamma +1}}\left(1+{\displaystyle \frac{\gamma -1}{\gamma }}{M}_{py}^2\right)\right)}^{{\displaystyle \frac{\gamma +1}{2(\gamma -1)}}}{\varphi }_p$$

(8)

$$M_{py}^2=\left({\left({\displaystyle \frac{P_p}{P_{py}}}\right)}^{{\displaystyle \frac{\gamma -1}{\gamma }}}-1\right){\displaystyle \frac{2}{\gamma -1}}$$

(9)

where $P_{py}$ is the pressure of the primary flow at the section $y-y$ and ${\varphi }_p$ is a factor to account for losses due to viscose effect (typical value is 0.88). The mixing is assumed to take place at section $y-y$, and the pressure of the two streams are equal at that section, i.e., $P_{py}=P_{sy}$, where $P_{sy}$ is the pressure of the secondary flow at section $y-y.$ Using $M_{sy}=1$ in the following equation, $P_{sy}$ is given by

$$P_{sy}=P_a{\left(1+{\displaystyle \frac{\gamma -1}{\gamma }}{M}_{sy}^2\right)}^{{\displaystyle \frac{-\gamma }{\gamma -1}}}=P_a{\left({\displaystyle \frac{2}{\gamma -1}}\right)}^{{\displaystyle \frac{-\gamma }{\gamma -1}}}$$

(10)

The outlet flow of the ejector is the sum of the primary and the secondary flows, hence, the rate and the temperature of the flow out of the ejector, $W_e$and $T_e,$are given by

$$W_e=W_p+W_s$$

(11)

$$T_e={\displaystyle \frac{W_p{T}_p+W_s{T}_s}{W_e}}$$

(12)

3.1.4. Separator

Even though the membrane is considered to be a gas barrier, air can diffuse from the cathode to the anode through the thin membrane. Therefore, in addition to the water crossed due to gradient effect, the oxygen diffused to the anode side reacts with hydrogen and causes water accumulation in the anode side. The presence of accumulated water disrupts the homogeneous distribution of hydrogen throughout the channel, and it can cause blockage of the gas channels. To prevent blockage, it is essential to eliminate any liquid water droplets that have the potential to cause such obstructions. This can be achieved by employing a water separator, which removes the water droplets from the recirculation stream. The dynamic equations of pressure for the different species of gases inside the water separator are detailed in Eq.(13) to (16)

$${\displaystyle \frac{{dp}_{sp,i}}{dt}}={\displaystyle \frac{R{T}_{sp}}{V_{sp}{M}_i}}(W_{i,in}^{sep}-W_{i,bleed}^{sep}-W_{i,s}^{sep})$$

(13)

where $W_{i,in}^{sep},W_{i,bleed}^{sep}, W_{i,s}^{sep}$ $(i:N_2,H_2)$ are respectively the inlet and outlet flow rates of the different species.

The vapor pressure inside the separator is the saturation pressure. It is assumed that the liquid water is in equilibrium with its vapor. The mass of vapor water inside the separator can be estimated using the ideal gas law and the saturation pressure of water at the given ambient temperature as in Eq.(14)

$$m_{H2O}={\displaystyle \frac{P_{sat}\left(T_{sep}\right)V_{sep}}{{RT}_{sep}}}$$

(14)

The individual mass flows of the separator outlet flow can be calculated based on the mass fraction of the different species as Eq.(15)

$$W_{i,j}^{sep}={\displaystyle \frac{m_{i,sep}}{m_{H_2,sep}+m_{N_2,sep}+m_{H_{2}O,sep}}}{W}_j$$

(15)

where $\left(i:H_2,H_{2}O,N_2\right)$, and the parameter $j,$stand for bleed or Ejector secondary. The mass of a species gas can be calculated using the separator pressure, temperature and volume as Eq.(16)

$$m_i={\displaystyle \frac{P_{i, sep}{V}_{sep,}{M}_i}{R{T}_{sep}}}, i\in [H_2,N_2]$$

(16)

3.1.5. Bleed Valve

During the recirculation process of the exhaust gas in the anode module, the concentration of nitrogen in the anode gas volume progressively builds up over time. This accumulation of nitrogen can cause a shortage of hydrogen supply to the fuel cell stack, leading to hydrogen starvation. To address this issue and reduce the nitrogen concentration in the anode module, a bleed valve is used, which allows for the opening of the closed anode loop. The valve opening is dependent on the diffusion coefficients of the materials used in the fuel cell stack and the operating conditions. At higher loads, higher inlet mass flows are required in the stack, thus more air diffuses through the membrane, which in turn leads to a faster accumulation of nitrogen. The nozzle flow equation is used to calculate the outlet flow of the manifold${ (W}_{out,bleed})$. The flow rate exiting a nozzle is influenced by the pressure differential between the upstream and downstream sides of the nozzle. The flow characteristic is divided into two regions depending on the critical pressure ratio$\left({\displaystyle \frac{p_{atm}}{p_{sep}}}\right)$. In the case of subcritical flow, when the pressure drop across the nozzle is lower than the critical pressure ratio, the mass flow rate can be calculated by Eq.(17)

$$\begin{gathered} W_{out,bleed}={\displaystyle \frac{C_D{A}_{T,bleed}{p}_{sep}}{\sqrt{R{T}_{sep}}}}{\left({\displaystyle \frac{p_{atm}}{p_{sep}}}\right)}^{{\displaystyle \frac{1}{\gamma }}}{\left\{{\displaystyle \frac{2\gamma }{\gamma -1}}\left[1-{\left({\displaystyle \frac{p_{atm}}{p_{sep}}}\right)}^{{\displaystyle \frac{\gamma -1}{\gamma }}}\right]\right\}}^{0.5} for \left({\displaystyle \frac{p_{atm}}{p_{sep}}}\right) \\ >{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma }{1-\gamma }}} \end{gathered}$$

(17)

where $\gamma$ is the ratio of specific heat capacities of gas, ${\displaystyle \frac{C_p}{C_v}}$, and the critical pressure ratio is 0.581. $C_D$ is the discharge coefficient of the nozzle, $A_{T,bleed}$ is the opening area of nozzle ($m^2$), and $R$ is the universal gas constant. In the case of critical flow, the mass flow rate is given by Eq.(18)

$$W_{out,bleed}={\displaystyle \frac{C_D{A}_{T,bleed}{p}_{sep}}{\sqrt{R{T}_{sep}}}}{\gamma }^{{\displaystyle \frac{1}{2}}}{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma +1}{2(\gamma -1)}}} for \left({\displaystyle \frac{p_{atm}}{p_{sep}}}\right)\le {\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma }{1-\gamma }}}$$

(18)

where $A_{T,bleed }$is the bleed orifice cross section [30]. This orifice cross section is used as a manipulated variable to control the flow rate of bleed and it is determined by the controller corresponding to the defined purpose.

3.2. Modeling the Air Supply System

The cathode or air supply subsystem supplies air with the required mass flow rate and pressure as determined by the supervisory controller. To simplify the control-oriented model, the cathode subsystem is represented as a linear multi-input multi-output system. This model takes the mass flow rate setpoint and air pressure setpoint as inputs and provides the corresponding outputs of mass flow rate and pressure. The cathode subsystem is linearized around the operating point 50 ${\displaystyle \frac{g}{s}}$ of mass flow rate, and $1.2 bar$ of pressure [25]. The linearized model is as follows.

$$\left[\begin{array}{c}{\dot{m}}_{air}\\ p_{air}\end{array}\right]=\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right]\left[\begin{array}{c}{\dot{m}}_{setpoint}^{air}\\ p_{setpoint}^{air}\end{array}\right]$$

(19)

$$G_{11}={\displaystyle \frac{2.797}{s(s+2.8)}},{ G}_{12}={\displaystyle \frac{0.0006592 s^2-0.001363 s+0.0007042}{s^3-2.979 s^2+2.959 s-0.9796}}$$

(20)

$$G_{21}=0, G_{22}={\displaystyle \frac{180.9442}{s(s+181.8496)}}$$

(21)

Humidifier Model

In this study, the membrane humidifier is modeled as a first-order system with a valve. The input of this model is the humidity ratio of the stack outlet (wet side),$H_w$, and the output is the humidity ratio of the inlet air (dry side) $H_d$.

$${\displaystyle \frac{H_d}{H_w}}={\displaystyle \frac{K_H{K}_v}{{\tau }_{H}s+1}}$$

(22)

where $K_H$ and ${\tau }_H$ are the gain and time constant, respectively. The control of humidification is done using a valve with an open rate $K_v$

$$K_v(t+1)=\left\{\begin{array}{c}{K}_v\left(t\right)+T_s\times Open rate, H_d<H_{d,setpoint}\\ K_v\left(t\right)-T_s\times Open rate, H_d>H_{d,setpoint}\end{array}\right.$$

(23)

where $H_d$, $H_w$ and $H_{d,setpoint}$ are the stack inlet air humidity ratio, stack outlet humidity ratio and stack inlet humidity ratio setpoint. The humidity ratio setpoint is determined by the supervisory controller.

3.3. Thermal Subsystem

The thermal subsystem is modeled as a first-order dynamic system, with the temperature of the fuel cell as its output. The simplified model is obtained based on the Inn-Balance experimentally validated Simulink model [26]. The proposed model is fit around the operating point with $T_{stack}=68℃$. The simplified dynamics is given by Eq.(24)

$${\displaystyle \frac{d{∆T}_{stack}}{dt}}={\displaystyle \frac{1}{{\tau }_{Th}\left(I_{stack}\right)}}({∆T}_{stack,setpoint}-{∆T}_{stack})$$

(24)

where ${∆T}_{stack}$ is the temperature over the inlet coolant flow temperature, which is assumed to be constant at $68\mathrm{℃}$. ${∆T}_{stack,setpoint}$ is the setpoint that is given by the NMPC. The nonlinear response time, ${\tau }_{Th}\left(I_{stack}\right)$, depends on the current and shows a polynomial behavior as follows:

$$\begin{gathered} {\tau }_{Th}\left(I_{stack}\right)=1.9556\times 10^{-18}{I_s}^7-4.3646\times 10^{-15}{I_s}^6+3.9468\times 10^{-12}{I_s}^5 \\ -1.8265\times 10^{-9}{I_s}^4+4.4546\times 10^{-7}{I_s}^3-5.2709\times 10^{-5}{I_s}^2 \\ +3.0746\times 10^{-3}{I}_s+1.1140\times 10^{-3} \end{gathered}$$

(25)

The temperature of the stack is given by Eq.(26)

$$T_{stack}={∆T}_{stack}+68\mathrm{℃}$$

(26)

The response time is found after fitting a step response to a linear first order system with unitary static gain. This is done for different step values of stack current$I_{stack}$.

3.4. Fuel Cell Stack Model

The next sections will provide an explanation of the control-oriented model of the stack used in the MPC. This model includes the gas channels dynamics, effect of nitrogen crossover through the membrane and water content in the membrane. Furthermore, the effect of liquid water in the cathode side and its effect on the electrochemical effective surface area is considered.

3.4.1. Cathode and Anode Channel

The mass balance equation characterizes the flow dynamics of gas species within the gas channels, leading to the determination of the concentration of the different gas species as follows [31].

$${\displaystyle \frac{\partial c_i}{\partial t}}=-{\displaystyle \frac{\partial }{\partial t}}\left(v{c}_i\right)-{\displaystyle \frac{{\dot{n}}_i}{{\delta }^j}}$$

(27)

$$v=-K^j(p_{in}-p_{out})$$

(28)

$$p=R{T}_{stack}\sum _i{c}_i$$

(29)

where $c_i$ is the concentration of the $ith$ gas, $v$ and $p$ are the gas velocity and the total pressure along the channel. The $i$ subscript is used for the gas species${ H}_2$, $O_2$, $N_2$ and $H_{2}O$. $T_{stack}$ is the temperature of the fuel cell. The parameters ${\delta }_j$ and $K_j$ are the channel height and the constant that relate the pressure drop with the gas velocity in the channels. The superscript $j$ is used for the cathode, $C$, and the anode, $A$. The term ${\dot{n}}_i$ comes from the transport of species to the GDL, whose direction is perpendicular to the gas channels. The pressures inside and at the outlet of the channel are $p_{in}$ and $p_{out}$, respectively. To simplify the model, the GDL is not modeled and the reactants consumption rate at the catalyst layer(CL) is modeled as follows.

$$W_{O_2,react}=M_{O_2}{\displaystyle \frac{n{I}_{st}}{4F}}$$

(30)

$$W_{H_2,react}=M_{H_2}{\displaystyle \frac{n{I}_{st}}{2F}}$$

(31)

where $M_{O_2}$ and $M_{H_2}$ are molar mass of oxygen and hydrogen, respectively. $F$ is Faraday constant and $n$ is the number of cells.

3.4.2. Nitrogen Crossover

The nitrogen crossover flux, $J_{N2,cross}$is determined by the nitrogen partial pressure gradient between the two sides of the membrane as Eq.(32)

$$J_{N2,cross}=K_{N2}{\displaystyle \frac{p_{N2,c}-p_{N2,a}}{{\delta }_m}}$$

(32)

where, $p_{N2,c}$, $p_{N2,a}$ are the partial pressure of nitrogen in the cathode catalyst layer (CCL) and the anode catalyst layer (ACL), respectively. ${\delta }_m$ is membrane thickness. $K_{N2}$ is the permeability of the membrane, which is determined by the temperature and water content in the membrane [32].

$$K_{N2}\left(W_m,T_{fc}\right)=(0.0295+1.2F_v-1.93{F_v}^2)\left({10}^{-11}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{E_{N2}}{R}}\right({\displaystyle \frac{1}{T_{ref}}}-{\displaystyle \frac{1}{T_{fc}}}\left)\right)$$

(33)

where $E_{N2}=24kJ{mol}^{-1}$ is the nitrogen molar energy and $F_v$is the volume fraction of water in the membrane that is given by Eq.(34)

$$F_v={\displaystyle \frac{W_m{V}_w}{V_m+W_m{V}_w}}$$

(34)

where $W_m$, $V_w$ and $V_m$are water contents of membrane, molar volume of liquid water and dry membrane, respectively.

3.4.3. Membrane Water Transport and Water Content

The flow rate of water vapor through the membrane, from the anode to the cathode, caused by back diffusions and electro-osmotic drag force can be expressed as Eq.(35) to (40)

$$W_{am}={\alpha }_{net}{N}_{cell}{\displaystyle \frac{I_{stc}{M}_{H2O}}{F}}$$

(35)

where ${\alpha }_{net}$ is the water transfer coefficient and it is given by Eq.(36)

$${\alpha }_{net}={n_d-{\displaystyle \frac{F{A}_{fc}}{I_{stak}}}{D}_w{\displaystyle \frac{{\rho }_{m,dry}}{t_m{M}_{m,dry}}}({\lambda }_{ca}-{\lambda }_{an})\alpha }_{net}$$

(36)

where $A_{fc}$ is the active area, ${\rho }_{m,dry}$ and $M_{m,dry}$ are the density and the weight per mole of the dry membrane. The parameter $t_m$ is the membrane thickness, $n_d$ the electro-osmatic drag coefficient, and $D_w$ the diffusion coefficient. They are given by Eqs.(37) and (38)

$$n_d=0.0029{{\lambda }_{an}}^2+0.0029{\lambda }_{an}-3.4\times {10}^{-19}$$

(37)

$$D_w=D_y\mathrm{e}\mathrm{x}\mathrm{p}\left(2416\left({\displaystyle \frac{1}{T_{ref}}}-{\displaystyle \frac{1}{T_{stack}}}\right)\right)$$

(38)

where $D_{\lambda }$ is a diffusion coefficient related to anode water content defined by Eq.(39) [3]

$$D_y=\left\{\begin{array}{c}\begin{array}{cc}{10}^{-10}& {\lambda }_{an}<2\end{array}\\ \begin{array}{cc}{10}^{-10}\left(1+2\left({\lambda }_{an}-2\right)\right)& { 2\le \lambda }_{an}<3\end{array}\\ \begin{array}{cc}{10}^{-10}\left(3-1.67\left({\lambda }_{an}-3\right)\right)& {3\le \lambda }_{an}<4.5\end{array}\\ \begin{array}{cc}1.25\times {10}^{-10}& {4.5\le \lambda }_{an}\end{array}\end{array}\right.$$

(39)

where ${\lambda }_{an}$ and ${\lambda }_{ca}$ are the water content of the membrane in the anode and the cathode side. They are calculated based on the water activity as Eq.(40) [3]

$${\lambda }_i=\left\{\begin{array}{cc}0.043+17.8a_{H2O,i}-39.85{a^2}_{H2O,i}+36{a^{23}}_{H2O,i}& {a^{23}}_{H2O,i}<1\\ 14+1.4({a^2}_{H2O,i}-1)& {{1\le a}^{23}}_{H2O,i}<3\\ 16.8 & {a^{23}}_{H2O,i}\ge 1\end{array}\right.$$

(40)

where $i$ stand for anode and cathode. The water activity is defined based on partial pressure of water and the saturation pressure as $a_{H2O,i}={\displaystyle \frac{P_{H2O}^i}{P_{sat}}}$.

3.4.4. Voltage Equation

The output voltage of the fuel cell stack is calculated by considering the Nernst equation as Eq.(41)

$$V=E_{nernst}-V_{act}-V_{ohm}$$

(41)

where $V$ is the stack output voltage, $E_{nernst},$ is the reversible Nernst voltage, $V_{act}$, is the activation loss and it is obtained using the Tefal equation and $V_{ohm}$ is the Ohmic voltage drop. Therefore, the output voltage is calculated according to Eq.(42) [30]

$$V=E_{nernst}-{\displaystyle \frac{R{T}_{stck}}{\alpha 2F}}\left(log\left({\displaystyle \frac{I_m}{I_0}}\right)-log\left({\displaystyle \frac{P_{O2}}{P_{O2,ref}}}\right)\right)-R_{Ohm}{I}_m$$

(42)

$$E_{nernst}=1.229-0.85\times {10}^{-3}\left(T_{stack}-298.15\right)+4.3085\times {10}^{-5}\left(ln\left(p_{H2}\right)+ln\left(p_{O2}\right)\right)$$

(43)

where $p_{H2}$,$p_{O2}$ are the partial pressures of hydrogen and oxygen. The exchange current density at the cathode $I_0$ is a function of the fuel cell temperature $T_{stack}$, oxygen pressure at the CCL, and the electrochemically active surface area at the CCL ($A_{ECSA}$). It is calculated as follows [31].

$$I_0=i_{0,ref}{\displaystyle \frac{A_{ECSA}}{A_{geo}}}{\left({\displaystyle \frac{p_{O_2}}{p_{O_2,ref}}}\right)}^{0.5}{e}^{\left[-{\displaystyle \frac{{\mathsf{\Delta }G}^{*}}{{RT}_{stack}}}\left(1-{\displaystyle \frac{T_{stack}}{T_{ref}}}\right)\right]}$$

(44)

where $i_{0,ref}$ is the intrinsic catalytic Pt activity at normal conditions$\left(T_{ref},p_{\mathrm{O}2,\mathrm{r}\mathrm{e}\mathrm{f}}\right)$. $A_{geo}$is the total surface area of the electrode, and ${∆G}^{\mathrm{*}}$is the Gibbs activation energy for the oxygen reduction reaction at the CCL. The electrochemically active surface area at the CCL ($A_{ECSA}$) is calculated as follows:

$$A_{ECSA}=A_{ref}{e}^{(k_{act}S-1)}$$

(45)

where $A_{ref}$ is the reference area multiplier and $k_{act}$ is the active area reduction coefficient. The liquid water in CCL is modeled with $S$.

3.4.5. Liquid Water and ECSA

In [33] the liquid water on the CCL is modeled with the mesoscopic pore filling effects of the CCL structure. The $S$ dynamics depends on the rate at which the liquid water is evaporated, $I_{evap}$, and the generation of liquid water on the CCL, ${\dot{n}}_{H_{2}O}^C$.

$$k_s{\displaystyle \frac{dS}{dt}}\left(z\right)=I_{gen}\left(z\right)-I_{evap}\left(z\right)={\dot{n}}_{H_{2}O}^C\left(z\right)-I_{evap}\left(z\right)$$

(46)

$$I_{evap}\left(z\right)=K_{evap}{\displaystyle \frac{s\left(z\right)\left(p^{sat}\left(z\right)-p^{vap}\right)}{RT{A}_{pore}}},if\left(p^{vap}\right(z)<p^{sat})$$

(47)

where $K_{evap}$ is the evaporation time constant,$p^{sat}$ is the saturation pressure and $p^{vap}$ is the vapor partial pressure at the GDL and $A_{pore}$ is the pore surface area per unit volume of the GDL/CL [34].

3.5. Anode Recirculation Loop Linear Model

In this section, the anode recirculation loop model is linearized around an operating point. This is done due to the highly nonlinear nature of the system’s dynamics. Moreover, it is necessary to discretize the model to enable its use in MPC. However, using a small sampling time increases the computational burden, particularly with a large prediction horizon. In this study, a slow NMPC approach is employed to determine the optimal setpoint for ${∆T}_{stack}$and the air inlet humidity ratio. Given the slow dynamics of temperature and humidity ratio, a large prediction horizon is utilized, and the linearized model of the anode recirculation system is applied to reduce the computational burden. To linearize a model, one defines the inputs and outputs, along with the operation points for linearization. The inputs and outputs are used to linearize the anode recirculation loop in this study are:

$$Input=\left[I_{stack},BG,p_{an}\right],Output=[H_{2in},N_{2in},{H_{2}O}_{in},p_{rm},H_{2tank}]$$

(48)

where, $I_{stack},BGand{p}_{an}$ are stack current, bleed gain rate and anode inlet pressure, respectively. The output vector includes, inlet hydrogen, $H_{2in}$, inlet nitrogen, $N_{2in}$, inlet water ${H_{2}O}_{in}$, return manifold pressure, $p_{rm}$ and finally hydrogen flow from the tank $H_{2tank}$Due to the complexity of the fuel cell system, the generated model is a large model with many states. However, MPC needs a reduced order model to have a low computational burden. Therefore, the balanced realization method is employed as a model reduction technique. The reduction is based on the Hankel singular values of the balanced model [35]. The linearized model of the anode recirculation loop is represented by the following transfer function matrixes.

$$\left[\begin{array}{c}{H}_{2in}\\ {\mathrm{N}}_{2in}\\ H_2{O}_{in}\\ p_{rm}\\ H_{2tank}\end{array}\right]=\left[\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ \begin{array}{c}{G}_{31}\\ G_{41}\\ G_{51}\end{array}& \begin{array}{c}{G}_{32}\\ G_{42}\\ G_{52}\end{array}& \begin{array}{c}{G}_{33}\\ G_{43}\\ G_{53}\end{array}\end{array}\right]\left[\begin{array}{c}{I}_{st}\\ \begin{array}{c}BG\\ p_{an}\end{array}\end{array}\right]$$

(49)

The state space matrixes of this linearized model are as follows:

$$A=\left[\begin{array}{c}-0.0007005 0.01161 0.03112 -0.005691\\ 0.01161 -0.2238 -0.9871 0.1526\\ 0.02662 -0.1445 -1.517 0.4096\\ 0.003876 -0.1038 -0.1672 -0.705\end{array}\right]$$

(50)

$$B=\left[\begin{array}{ccc}-0.002842& 0.3421& -0.09016\\ \begin{array}{c}0.08594\\ \begin{array}{c}0.1226\\ 0.02933\end{array}\end{array}& \begin{array}{c}-4.145\\ \begin{array}{c}-10.65 \\ -0.6559\end{array}\end{array}& \begin{array}{c}1.024\\ \begin{array}{c}-1.008\\ 1.261\end{array}\end{array}\end{array}\right]$$

(51)

$$C=\left[\begin{array}{cc}0.001217& \begin{array}{ccc}7.491e-5& -0.0001228& -0.0002258\end{array}\\ \begin{array}{c}1.945e-5\\ \begin{array}{c}-0.0001005\\ \begin{array}{c}-0.3538\\ 2.344e-5\end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}-6.73e-7\\ \begin{array}{c}-2.156e-5\\ \begin{array}{c}4.271\\ -0.0001574\end{array}\end{array}\end{array}& \begin{array}{c}3.797e-06\\ \begin{array}{c}-1.324e-5\\ \begin{array}{c}10.7 \\ -0.001162\end{array}\end{array}\end{array}& \begin{array}{c}-1.842e-05\\ \begin{array}{c}-0.0001181\\ \begin{array}{c}-1.422\\ 0.0001741\end{array}\end{array}\end{array}\end{array}\end{array}\right]$$

(52)

$$D=\left[\begin{array}{ccc}2.273\mathrm{e}-05& 0.001502& 0.00171\\ 1.845\mathrm{e}-07& 3.669\mathrm{e}-05& 3.799\mathrm{e}-05\\ \begin{array}{c}5.161\mathrm{e}-06\\ \begin{array}{c}-0.0005933\\ 9.785\mathrm{e}-05\end{array}\end{array}& \begin{array}{c}0.000241\\ \begin{array}{c}0.2264\\ -0.0001213\end{array}\end{array}& \begin{array}{c}0.0002766\\ \begin{array}{c}-1.032\\ 0.003369\end{array}\end{array}\end{array}\right]$$

(53)

4. Offline Generated Setpoints

The stack Simulink model, which has been validated through experiments, is simulated until it reaches a stable condition by testing all feasible combinations of valid inputs within a defined grid. This process creates an offline map of the optimal operating point corresponding to the load current. The input is the stack current and the outputs are the optimal setpoints. The compressor map, which links pressure ratio and mass flow rate, limits the setpoint values within the compressor’s operating space. Additionally, the compressor consumption is considered as a power loss in the system. The efficiencies are assessed for different setpoint groups, and the setpoint values that maximize efficiency are chosen as the optimal setpoints for each stack current. These offline maps for five operation points are shown in Figure 3.

The offline maps indicate that as the load current increases, there should be corresponding increases in mass flow rate, pressure, and temperature. Additionally, the bleed rate, which represents the percentage of the bleed valve’s opening, should also increase. This is because, as the load increases and more inlet air is introduced, the crossover of nitrogen also increases, leading to its accumulation in the anode recirculation loop. As a result, the bleed rate needs to be increased to open the bleed valve further and remove the accumulated nitrogen from the anode loop. Ultimately, this leads to an increase in the partial pressure of hydrogen. The generated water increases as the load current increase. Therefore, the setpoint of humidity ratio shows a decrease in the setpoint to prevent the overhydrate of the stack and flooding.

5. Nonlinear Model Predictive Control

This study aims to analyze various operating variables that impact the overall efficiency of the automotive fuel cell system. Additionally, it develops a controller to optimize their values. According to the described model, the system efficiency is influenced by the following variables: inlet air pressure, inlet air mass flow rate, inlet air humidity ratio, stack temperature, fuel pressure and bleed rate in the anode recirculation loop. To achieve the efficiency maximization, a nonlinear model predictive control(NMPC) strategy is proposed. This approach aims to determine the optimal values for these variables, ultimately improving the system’s efficiency. The NMPC utilizes the control oriented model (COM) of the system described in section 3 to predict its behavior and generate optimal setpoints. In order to implement dynamic equation in NMPC, finite difference method is used to solve numerically the dynamic equations. Applying the backward discretization method, the dynamic equations are discretized. Based on the system dynamics Eq.(1) - Eq.(46), the optimization problem using the NMPC is defined as follows.

$$\underset{U\left(k\right)}{\mathrm{Max}}J\left(X\left(k\right),U\left(k\right)\right)$$

(54)

$$X\left(k+i+1\right)=f\left(X\left(i+k\right),U(i+k)\right), 0<i\le N_p$$

(55)

$$U_{min}\left(i+k\right)\le U\left(i+k\right)\le U_{max}\left(i+k\right), 0\le i\le N_c$$

(56)

$$U\left(i+k\right)=U\left(k+N_c-1\right) N_c\le i\le N_p$$

(57)

$$X\left(k\right)=X_m\left(k\right)$$

(58)

where$N_c$, $N_p$ and $J$are the control and prediction horizon, $N_c\le N_p$, and cost function, respectively. The state vector $X\left(k\right)\in R^n$ represents the system’s dynamic states variables, which include pressure, concentration, CCL liquid water content, temperature, and humidity ratio. During each iteration of the NMPC, the initial values for the state variables are derived from the measured values $X_m\left(k\right)$. The input control vector $U\left(k\right)\in R^m$, on the other hand, acts as the independent optimization variables and is obtained during each iteration of the optimization problem.

The proposed NMPC architecture consists of three separate nonlinear model predictive controllers (NMPCs), each one with its own control-oriented model (COM), manipulated variables (MVs), constraints, outputs, and configuration. The reason for splitting the controller into three separate NMPCs is to lower the computational workload and eliminate unneeded calculations. Indeed, the five setpoints operate on distinct time scales. Then, by dividing the controller according to their distinct time constants, the controller becomes both more efficient and effective in terms of computation burden.

The optimization variable vector $U\left(k\right)\in R^m$ is divided into three: $U_{Slow}\left(k\right)$, $U_{BG}\left(k\right)$ and $U_{Fast}\left(k\right)$ as Eq.(59)

$$U_{Fast}\left(k\right)=\left[\begin{array}{c}{p}_{ca}\\ {\dot{m}}_{ca}\end{array}\right], U_{BG}\left(k\right)=\left[BG\right], U_{Slow}\left(k\right)=\left[\begin{array}{c}{∆T}_{stack}\\ {HR}_{air}\end{array}\right]$$

(59)

where $p_{ca}$, ${\dot{m}}_{ca}$ are inlet air pressure, and inlet air mass flow, respectively. The slow dynamic variables are ${∆T}_{stack}$, and ${HR}_{air}$, which is the humidity ratio of inlet air. The bleed gain rate, referred to as $BG$, represents the proportion of the bleed valve’s opening. This optimization problem is solved using the fmincon toolbox of MATLAB.

5.1. Fast Dynamic NMPC:

The fast dynamics NMPC is based on the model, COM-1 which mainly includes the air supply system, the hydrogen supply system, and the stack model. The designed fast dynamic NMPC calculates the setpoints for the inlet air mass flow rate and inlet air pressure. During the prediction horizon in fast NMPC, the three other controlled variables are considered constant, each one maintaining its previous measurement value as follows.

$$BG\left(k+i\right)={BG}_m\left(k-1\right), 0<i\le N_p^F$$

(60)

$$\Delta T\left(k+i\right)={\Delta T}_m\left(k-1\right), 0<i\le N_p^F$$

(61)

$$RH\left(k+i\right)={RH}_m\left(k-1\right), 0<i\le N_p^F$$

(62)

where $N_p^F$ is the fast dynamic NMPC prediction horizon. Given the slow-changing nature of the relative humidity and temperature, it is reasonable to assume that they remain constant over the short time frame of the fast dynamic prediction horizon. The other setpoint is the bleed rate, which also is constant over the fast dynamic prediction horizon. It is assumed that the bleed valve has a perfect controller that sets the valve opening percentage based on the given setpoint. During the fast prediction horizon, the position of the valve remains constant.

5.2. Slow Dynamic NMPCs:

The stack temperature, inlet air humidity ratio and bleed rate have larger time constants relative to the mass flow rate and pressure. Therefore, they have been separated and considered as the slow dynamics. The NMPC is designed to calculate the setpoints for the stack temperature and air humidity in one hand and the BG on the other hand. Due to the slow dynamic of temperature and humidity ratio, the prediction horizon, $N_p^S$, must be large enough. Therefore, it is assumed that the fast dynamic variables of inlet air mass flow rate and inlet air pressure are in the steady state conditions as follows.

$$\dot{m}\left(k+i\right)={\dot{m}}_{setpoint}\left(k\right), 1<i\le N_p^S$$

(63)

$$p\left(k+i\right)=p_{setpoint}\left(k\right), 1<i\le N_p^S$$

(64)

where $N_p^s$ is the slow dynamic NMPC prediction horizon. During the prediction horizon $N_p^S$, the ${\dot{m}}_{setpoint}$ and $p_{setpoint}$are calculated using the offline setpoint maps was already explained in section 4. Furthermore, a large prediction horizon increases the computational burden. Therefore, the COM-2 includes the stack dynamics in the steady state conditions, air supply system dynamic, thermal subsystem dynamic, humidifier dynamic and liner model of hydrogen recirculation system.

The slow dynamics of nitrogen accumulation in the anode recirculation loop make frequent adjustments to the bleed rate unnecessary. As a result, the bleed rate setpoint is updated less frequently than the setpoint of temperature and humidity ratio. This different control period is the reason to have two slow NMPC. As indicated in Table 1, the first slow NMPC is executed every $T_{slow1}$, updating the temperature and humidity ratio setpoint. The second slow NMPC is executed every $T_{slow2}$ to update the optimal bleed rate setpoint. The COM-3 for the second slow NMPC includes, stack dynamics, air supply system dynamics, hydrogen recalculation nonlinear model, thermal subsystem dynamic and humidifier dynamic.

5.3. Cost Function

The control objectives are as follows:

• Maximizing the overall system efficiency
• Minimizing the system degradation considering the manufacture limitations

The efficiency is defined as the ratio of the net generated energy to the total energy of the fed hydrogen. The efficiency, as the cost function, is defined over prediction horizon in Eq.(65)

$$J\left(k\right)=\eta ={\displaystyle \frac{{\displaystyle {\sum }_{i=1}^{Np}}{(IV}_{fc}\left(i\right)-P_s\left(i\right))}{HHV(m_{H2}^{i=1}-m_{H2}^{i=Np})+{\displaystyle {\sum }_{i=1}^{Np}}HHV.{\dot{n}}_{H_2}\left(i\right)}}$$

(65)

where ${\dot{n}}_{H_2}$is the $H_2$flow rate from the $H_2$tank, and $HHV$ denotes the higher heating value of the $H_2$. ${\mathrm{I}\mathrm{V}}_{fc}$ is the gross electrical power generated by the fuel cell, while $P_s$ represents the electrical power consumed by the subsystems. This auxiliary consumption is assumed to be the power consumption of the compressor, which is dependent on the air mass flow rate and pressure. The system’s stored hydrogen is determined by subtracting the initial hydrogen value in the various volumes, including the anode channel, water separator, supply and return manifold of the system,$m_{H2}^{i=1}$, from the final value over the prediction horizon $m_{H2}^{i=Np}$.

5.4. Fast Dynamic Variables Constraints

The manufacturing limitations and safety requirements, such as the compressor operation points map and the hard constraints related to the parameters are modeled as constraints. The defined setpoints have upper and lower limits and these limits must be considered in the control design. The air mass flow rate setpoint and the air pressure setpoint must be kept inside the compressor map margins. To increase the convergence time and decrease the computation burden, a box centered around the obtained optimal setpoint in the offline setpoint generator, i.e., ${\dot{m}}_{ref}$ and $p_{ref}$ is considered. The upper bound is defined as Eq.(66)

$$UB\left(k\right)=\left[{\dot{m}}_{ref}\right(k)+10,p_{ref}(k)+0.5]$$

(66)

The lower bound of the mass flow rate setpoint is defined as the upper bound. However, it is also important for the mass flow rate setpoint to satisfy the minimum stoichiometry requirement. Therefore, the lower bound of these two setpoints is defined as follows.

$$LB\left(k\right)=[\max \left({\dot{m}}_{ref}\left(k\right)-10, {\displaystyle \frac{{n_{cell}\lambda }_{o2min}{M}_{o2}}{0.21\times 4F}}{I}_{stack}\right),p_{ref}\left(k\right)-0.5]$$

(67)

The box considered in Eq.(66) and (67) defines the upper and lower constraints. However, in certain load current scenarios, the vertices of the box may fall outside the range of the compressor map. In such cases, linear inequalities are employed to correct this deviation. The obtained linear inequality is as Eq.(68)

$$n_{geq}\left[\begin{array}{c}{\dot{m}}_{ref}^3\left(k\right)\\ {\dot{m}}_{ref}^2\left(k\right)\\ {\dot{m}}_{ref}^1\left(k\right)\end{array}\right]-p_{ref}\left(k\right)\le -l_{geq}$$

(68)

The process is repeated for the upper bound. If the upper vertices are outside the permissible range, a tangent line to the surge line on the left side of the compressor map is determined as Eq.(69)

$$n_{leq}\left[\begin{array}{c}{\dot{m}}_{ref}^2\left(k\right)\\ {\dot{m}}_{ref}^1\left(k\right)\end{array}\right]-p_{ref}\left(k\right)\ge l_{leq}$$

(69)

where $n_{geq}\mathrm{a}\mathrm{n}\mathrm{d}$ $n_{leq}$ are the coefficient of the defined line in the lower and upper margins, respectively.

$$n_{geq}=\left[-4.3984\times {10}^{-7}2.1305\times {10}^{-4}-1.2157\times {10}^{-3}\right],l_{geq}=1.0171$$

(70)

$$n_{leq}=\left[4.4201\times {10}^{-5}-1.8396\times {10}^{-3}\right], l_{leq}=1.0288$$

(71)

5.5. Slow Dynamic Variable Constraints

The slow dynamic setpoints are treated as independent variables, with their upper and lower limits defined. These setpoints encompass the inlet air humidity ratio, thermal difference temperature, and bleed rate. The upper and lower limits of the humidity ratio are established based on the operating range of the inlet air pressure, temperature, and relative humidity. To ensure stable operation and prevent issues like dryness or flooding, the relative humidity is maintained within a range of 30% to 55%. By considering the relationship between the humidity ratio and relative humidity, as expressed in Eq.(72), and taking into account the air pressure range, the upper limit of the humidity ratio can be determined as (0.01<HR<0.12) [30].

$$HR=0.62{\displaystyle \frac{RH\times P_{sat}}{\left(p-RH\times P_{sat}\right)}}$$

(72)

where $P_{sat}$ is the vapor saturation pressure and depends on the temperature. Due to the limitations on the maximum and minimum pumping mass flow of the coolant circuit, it is not possible to achieve all possible temperature references. Figure 4 illustrates the temperature range achievable for each load current.

To determine the boundaries, certain adjustments were made. The inlet coolant temperature was set to $68°C$ and the reference for the outlet coolant temperature difference was modified. The lower boundary was obtained by setting $∆T= 0 °C$. Similarly, the upper boundary was obtained by setting $∆T=12 °C$. Observing the Figure 4, it shows that at lower currents, the fuel cell heat is unable to reach the desired $+12 °C$ setpoint. Conversely, at higher currents, the coolant circuit cannot achieve the targeted $+0 °C$ setpoint. Consequently, considering the stack current, both the upper and lower limits of the setpoint need to be adjusted between the limits.

The bleed rate, $BG$, is the percentage of the opening and it is a number between 0% to 100% opening.

5.6. Integration of NMPC with the Fuel Cell System

The block diagram of the proposed structure is shown in Figure 5. The block diagram shows the three NMPCs, their COMs, inputs and outputs, prediction and control horizon, COM sampling time, and the relevant updating time of the setpoints.

The stack current is a disturbance for the system and needs to be predicted within the prediction time horizon. A simple method is used to predict the current within the prediction horizon. This method relies on calculating the rate of change between two consecutive values of current during a $T_{predict}$ second, and it is employed to create the future current profile, as outlined in Eq.(73):

$$I_{stack\_predict}\left(k+i+1\right)=I_{stack}\left(k+i\right)+slop{\mathrm{T}}_s$$

(73)

$$slop={\displaystyle \frac{I_{stack}\left(k\right)-I_{stack}\left(k-T_{predict}\right)}{T_s}}$$

(74)

Due to difference in COM of three NMPCs, three different sampling time, ${\mathrm{T}}_{s1},{\mathrm{T}}_{s2},{\mathrm{T}}_{s3}$ are used. Therefore, corresponding to the type of NMPC that is used, the relevant sampling time,${\mathrm{T}}_{si}(i=\mathrm{1,2},3)$, is used in the Eq.(73).

6. Simulation Results

In this section, the purpose is to investigate how the online optimal setpoint generator, which in this study is a NMPC, updates the operating parameters setpoint in order to maximize the overall system efficiency. The load profile is obtained from the Artemis highway cycle. Furthermore, the performance of the balance of plant subsystem local controllers in terms of setpoint tracking is considered and evaluated. In this research, a PI controller, ( $Kp=2\times {10}^{-4}$, $K_I=1$), is used to regulate the anode inlet pressure, maintaining a constant value above the cathode inlet pressure. This approach is adopted for safety reasons and to prevent physical damage to the membrane. Specifically, a pressure difference of $∆p=0.3$ bar is maintained between the anode inlet and cathode inlet pressures.

Characteristics of the proposed model predictive controllers are summarized in Table 1.

6.1. Bleed Rate Effect on Fuel Cell Efficiency

In order to minimize hydrogen losses and maintain efficiency, precise control of the bleed rate is necessary. On one hand, increasing the bleed rate leads to higher hydrogen loss and decreased efficiency. On the other hand, it helps reduce nitrogen accumulation, increase hydrogen partial pressure, and ultimately enhance efficiency. In Figure 6, the effect of bleed rate for a constant load current and constant air mass flow rate, pressure and humidity ratio is shown. The simulation was done in long enough time that the states reaches their steady state value. This simulation shows the importance of bleed rate control. As indicated by the vertical red dashed lines at t = 289 and t = 358, the system efficiency reaches its maximum value when the bleed gain is at its optimal value. Increasing the bleed gain beyond this optimal value leads to a decrease in system efficiency due to increased hydrogen loss and recirculation flow. Despite the rise in hydrogen partial pressure, the system efficiency diminishes. Conversely, if the bleed gain value is decreased from the optimal value, nitrogen (N2) begins to accumulate, resulting in a decrease in hydrogen partial pressure and ultimately leading to a decrease in efficiency.

6.2. Generated Optimal Setpoints for a Dynamic Load Using NMPC

In this section, the designed structure is tested by applying a dynamic load current. The objective is to assess how the model predictive controller generates the five setpoints for the load current in real-time. Figure 7 shows the generated setpoint for the fast dynamic setpoints, specifically, the inlet air mass flow rate and pressure. The other three setpoints, bleed gain rate, humidity ratio and temperature difference are shown in Figure 8. The setpoint behavior indicates that as the load current increases, the setpoints generally rise as anticipated, with the exception of the humidity ratio, which should decrease to prevent flooding in the stack.

6.3. Comparison of Off-Line and NMPC Strategies in Terms of the System Efficiency Improvement

In order to compare the proposed method and analyze its effectiveness, the average efficiency improvement is assessed by comparing the NMPC with offline setpoint maps obtained through steady state analysis in section 4. In the comparison, it is evident that the average efficiency, calculated based on Eq.(75), is higher when using the model predictive control. Nevertheless, the computational time is longer with the model predictive control. Despite the increased computational burden, the real time measurement of the variables has the effect of disturbance and system uncertainty. Therefore, the model predictive control method is more robust in real-time than the offline setpoint generator, which relies on a constant steady state model. Figure 9 shows the average efficiency comparison between the offline setpoint generator and MPC.

$${\eta }_{ave}={\displaystyle \frac{{\int }_{t=t_0}^{t_{ave}}{(IV}_{fc}\left(\tau \right)-P_s\left(\tau \right))d\tau }{HHV\left(m_{H2}^{t=t_0}-m_{H2}^{t=t_{ave}}\right)+{\int }_{t=t_0}^{t_{ave}}{(IV}_{fc}\left(\tau \right)-P_s\left(\tau \right))d\tau }}$$

(75)

The averaging is calculated during $t_{ave}=1s$. Because MPC relies on model predictions within a prediction horizon, having accurate load current predictions for the future is crucial for achieving the best results. However, in automotive applications, load current is not known in advance and depends on driver actions. Nevertheless, researchers are working to predict future load currents based on historical load data [36]. These predictions are more suitable in vehicles doing regular stablished trip, as buses. Figure 10 illustrates the average efficiency when assuming that load current is known in advance and compares it with two other scenarios. The simulation results demonstrate that the average efficiency is better when load current for the future is known. Figure 10 shows the average efficiency in the steady state condition after the transient behavior of the simulation.

7. Conclusions

In this study, a nonlinear model predictive controller (NMPC) was developed to optimize the control variables of a fuel cell system with the objective of enhancing its average efficiency. The fuel cell system includes components such as the anode recirculation loop with the ejector and water separator. Recognizing the varying time scales of the selected control variables, our proposed NMPC was structured into three distinct parts, each tailored to its own prediction horizon and sampling time. We categorized air pressure and mass flow rate as fast dynamic input variables, while temperature and humidity ratio were designated as slow dynamic input variables. To minimize hydrogen loss and maintain efficiency, a third NMPC was designed to control the setpoint of the bleed valve, which regulates the percentage of valve opening.

A comparative analysis of the system’s average efficiency between our proposed structure and an offline setpoint generator was conducted. The simulation results indicated that the NMPC consistently demonstrated superior average efficiency when compared to the offline setpoint generator. Additionally, the average efficiency is better when assuming that load current is known in advance in comparison with the scenario that load current is considered constant during the prediction. This results highlights the advantage of load current predictions within the prediction horizon over the constant load approach for optimizing efficiency.