Figure 1 shows that within a PEMFC, a cell comprises essential components, encompassing bipolar plates, a proton exchange membrane (PEM), gas diffusion layers (GDLs), and catalyst layers. At the anode, the catalytic decomposition of hydrogen gas generates protons and electrons. Protons traverse the PEM to the cathode, where they combine with oxygen to form water and generate heat. Simultaneously, electrons liberated at the anode traverse an external circuit to the cathode, generating an electric current that powers the load.

The format of the formulas involved in the article has been completely modified.

The expression of voltage for an individual cell within a PEMFC is given as follows: V cell = E Nernet − V ohm − V con − V act , (1) where E Nernet is the reversible voltage drop, V ohm represents the Ohmic voltage drop, V con indicates the concentration voltage drop, and V act corresponds to the activation voltage drop. The term E Nernst can be expressed by the following formula (Wang et al., 2017): E Nernst = 1.229 − 0.85 × 10 − 3 T st − 299 + 4.3085 × 10 − 5 ⋅ T st ln P H 2 + 1 2 ln P O 2 , (2) where T st represents the temperature of the stack, P H 2 denotes the partial pressure of hydrogen, and P O 2 indicates the partial pressure of oxygen. Protons and electrons produce V ohm when encountering electrical resistance, which can be described as (Wang et al., 2014) V ohm = I st R M + R C , (3) where I st represents the stack current, R M stands for the equivalent resistance associated with proton transport, and R C denotes the equivalent resistance related to electron transport. The calculation procedure for determining R M is outlined as follows: R M = 181.6 1 + 0.03 I st A + 0.062 T st 303 2 I st A 2.5 l λ m − 0.634 − 3 I st A exp 4.18 T st − 303 T st A , (4) where A is the area of the membrane, l stands for the water content of the membrane, and λ m denotes the thickness of the membrane. The decrease in concentration is given as follows: V con = − R T st 2 F ln 1 − I st I max , (5) where I max represents the maximum current, R denotes the gas constant, and F stands for the Faraday constant. The activation voltage drop is as follows: V act = − α 1 + α 2 + α 3 T st ⁡ ln C O 2 + α 4 T st ⁡ ln I st , (6) where α₁, α₂, α₃, and α₄ represent the parameters of the fuel cell, C O 2 denotes the concentration of oxygen and indicates that the concentration of oxygen can be determined by the temperature and partial pressure of oxygen, as is known from C O 2 = P O 2 101325 5.08 × 10 6 ⁡ exp − 498 T st . (7)

To facilitate the construction and analysis of the model, the following assumptions are made:

(1) Assuming that the chemical energy of the gas involved in the fuel cell reaction is exclusively converted into electrical and thermodynamic energy without dissipation into other forms of energy.

As shown in the Figure 2, the fuel cell temperature management model is established, mainly including the fuel cell, circulating water pump, cooling fan, water tank, and temperature sensor and controller.

Based on the thermal balance equation Q = CM Δ T , the heat equilibrium relationship of the fuel cell is as follows: C st M st dT st dt = Q gen − Q dis , (8) where C st represents the specific heat capacity of the reactor [kJ/(kg K)], M st denotes the mass of the reactor (kg), T st signifies the operating temperature of the reactor (K), Q gen represents the thermal power of the reactor (kW), and Q dis is the stack dispersion of thermal power (kW). The formula for Q gen can be expressed as follows: Q gen = Q tot − P st , (9) where Q tot is all chemical energy (kW) of the reactant itself per unit time, and P st denotes the output power of the system. Q tot = ∆ H N H 2 reacted , (10) N H 2 reacted = N cells I st 2 F , (11) where N H 2 reacted is the molar rate (mol/s) of hydrogen, I st is the actual output current of the stack, N cells is the number of pieces of the cell, and F is the Avogadro constant. P st = V stack I st , (12) Q dis = Q gas + Q cool + Q atm , (13) where Q gas is heat dissipation of gas in the system, Q cool is heat dissipation of the cooling water, and Q atm signifies radiation heat dissipation power of the reactor. Q cool = W cl C H 2 O I T st − T st in , (14) Q cool = AK T rad − T atm , (15) where A is the radiator heat exchange area ( m 2 ), K is the radiator heat transfer coefficient (W⁄ m 2 K ), T rad is the radiator heat exchange temperature (K), and T atm is the ambient temperature (K). For the radiator, the heat transfer temperature of the radiator is the arithmetic mean of the inlet and outlet temperatures of the radiator. T rad = 1 / 2 T rt + T rt in , (16) where T rt represents the radiator outlet temperature and T rt in indicates the radiator inlet temperature.

To validate the precision of the electrochemical model, we selected two different working temperatures, 333 K and 335 K, with specific parameters as shown in Table 1. By contrasting the simulation results with the experimental data, as depicted in Figure 3, it is evident that the simulation results align with the experimental data trend, confirming the reliability and effectiveness of the simulation model. Additionally, Figures 4A,B display the validation curve for the PEMFC current and temperature model and illustrate that with a step change in current, the temperature promptly responds to the variation in current, showing a high degree of agreement between the simulation curve and experimental data. This observation underscores the reliability and effectiveness of the temperature model.

Effective temperature control is imperative for optimizing fuel cell performance and ensuring long-term reliability. This study introduces a collaborative MPC strategy for managing fuel cell pumps and radiators and compares it with the traditional PID control and water pump MPC. Table 2 outlines the three control strategies, and Figure 5 illustrates their application to the fuel cell thermal management system. Figure 6A depicts PID control, a classic technique combining proportional, integral, and derivative components for temperature regulation. By contrast, Figure 6B showcases MPC employing predictive models to optimize inputs while considering constraints and dynamic control (Bressel et al., 2020; Rui et al., 2020). Implementing MPC for the fuel cell’s water pump and radiator allows a performance comparison against PID control. This analysis aims to identify the optimal temperature control strategy, enhancing fuel cell system efficiency and stability.

In fuel cell system research, we linearized the model to understand its dynamic response at various operating points. Focusing on crucial input parameters, such as cooling fan flow rate and coolant flow rate, and outputs like reactor input temperature and temperature difference, we derived two transfer functions through linear approximation. These functions effectively depict how coolant and fan flow rates impact reactor temperature. Using the temperature and flow rate as operating points allows the precise capture of the system’s dynamic behavior, which is a valuable tool for optimizing fuel cell performance. Linearizing the transfer function simplifies the complexity of non-linear systems, making system behavior more comprehensible. This aids in controller design, allowing the application of classical methods tailored for linear systems, ultimately optimizing system performance.

MPC excels over PID by offering predictive capability, handling constraints, and optimizing multivariable control. Its adaptability to system changes, optimal control approach, and effectiveness in handling non-linearities make it a superior choice for complex and dynamic systems. The pulse response model was utilized to forecast the output at future time points. y m k + j = ∑ i = 1 n h i u k + j − i , (17) where p is the prediction time domain (1 < j < p), h i i s  the finite impulse response , i represents the time step of the predicted time domain, and k is a certain point in the control process. The state of u(k + j) remains constant when i > (m − 1). u k + 1 = u k + m − 1 , (18) where m should be less than p, and m is the control time domain. The p-step forecasted values for the upcoming output values are y m k + j = ∑ i = 1 n h i u k + j − i , (19) where the value of j is taken from 1 to m − 1. y m k + j = ∑ i = 1 j − m + 1 h i u k + m − 1 + ∑ i = j − m + 2 n h i u k + j − i . (20) where the value of j in the above formula is the integer from m to p. The control effect can be categorized as either known or unknown.

The controlled effects when known are U 1 k = u k − n + 1   u k − n + 2 ⋯ u k − 1 1 × n − 1 T . (21)

The controlled effects when unknown are U 2 k = u k   u k + 1 ⋯ u k + m − 1 1 × m T , (22) Y M k = H 1 U 1 k + H 2 U 2 k . (23)

Because there are various disturbances in the actual control process, the system is controlled at all times, and the actual output value y (k) will have a certain error with the output prediction value y m (k) of the system prediction model. The error between the actual output value and predicted output value is expressed by the following formula: e k = y k − y m k . (24)

Correction against the predicted value of the model yields the following formula: y P k + j = y m k + j + β j y k − y m k , (25) y m k = ∑ i = 1 N h i u k − i . (26)

Then, the P-step prediction value of the system can be expressed as follows: Y P k = Y m k + β e k , (27) where β is the output error correction gain. β = β 1   β 2   ⋯ β P T , (28) Y P k = y P k + 1   y P k + 2   y P k + p 1 × P T . (29)

The predictive control does not require rapid tracking of the setpoint, rather it prompts a gradual convergence of the output toward the setpoint along a specific trajectory. The reference track is calculated using both the setpoint and current measurements of the process output and is expressed as follows: y r k + j = α j y k + 1 − α j y d , (30) α = e − T S T , (31) where T S is the sampling time and T represents time constants governing the reference track. Y r k = y r k + 1   y r k + 2   y r k + p 1 × p T . (32)

The target functions for the optimization control are min   J = ∥ Y P k − Y r k ∥ Q 2 + ∥ U 2 k ∥ R 2 , (33) = Y P k − Y r k T Q Y P k − Y r k + U 2 T k R U 2 k . (34)

In the Y P k , J = H 1 U 1 k + H 2 U 2 k + β e k − Y r k T Q H 1 U 1 k + H 2 U 2 k + β e k − Y r k ] + U 2 T k R U 2 k . (35)

Solving for the optimal control rate ∂ J ∂ U 2 k = 0 gets the optimal control rate: U 2 k = H 2 T Q H 2 + R − 1 H 2 T Q Y r k − H 1 U 1 k − β e k , (36) Q = diag q 1   q 2 ⋯ q P , (37) R = diag r 1   r 2 ⋯ r M , (38) where Q is the output tracking weight, and R is the input move weight.

The optimal control of time K is U 2 k = D T Y r k − H 1 U 1 k − β e k , (39) D T = 1   0 ⋯ 0 1 × M H 2 T Q H 2 + R − 1 H 2 T Q . (40)

The MPC method is used for water pump control and radiator control. When applying the MPC control to the water pump: T r , d k + j = α j T std k + 1 − α j T std , set , (41) where T std represents the temperature difference between the inlet and outlet of the PEMFC during the operation, while T std , set is the designated value for this temperature difference.

The predicted temperature T P , d expression is as follows: T P , d k + j = T m , d k + j + β j T std k − T m , d k . (42)

The optimal control objective function for the water pump is given as follows: min   J = ∥ T P , d k − T r , d k ∥ Q 2 + ∥ W cl ∥ R 2 , (43) where W cl means that the best control variable is the coolant flow.

When applied to radiator control: T r , in k + j = α j T in k + 1 − α j T in , set , (44) where T in k represents the inlet temperature of the PEMFC and T std , set is the designated set value for the inlet temperature. The predicted expression for the stack inlet temperature is as follows: T P , in k + j = T m , in k + j + β j T in k − T m , in k . (45)

The optimal control objective function for the radiator is given as follows: min   J = ∥ T P , in k − T r , in k ∥ Q 2 + ∥ W air ∥ R 2 (46) where W air means that the best control variable is the radiator air volume.

Figure 7A illustrates a step current varying from 100 A to 220 A, showcasing step changes in current over time. Figures 7B,E depict temperature variations in the stack and fluctuations in temperature difference between the stack's inlet and outlet under different control strategies. It is evident that for both the stack inlet and outlet temperatures, the MPC collaborative control strategy exhibits superior performance. Comparative analysis with the traditional PID control and water pump MPC reveals significantly improved overshoot and stability time, as detailed in Table 3. The MPC strategy notably enhances control of the temperature difference between the reactor inlet and outlet, reducing fluctuation by over 50% and achieving a shorter settling time. This heightened control performance underscores the effectiveness of the MPC strategy in maintaining stability and precision.

In reactor thermal management, the controller directly regulates the coolant flow rate and cooling fan air volume as crucial control variables, ensuring rapid response for desired reactor inlet and outlet temperatures. Figures 7C,D reveal that under MPC, the coolant flow and radiator air volume adeptly track current changes, optimizing heat dissipation and maintaining temperature balance. Compared to PID-controlled gradual flow changes, pump MPC and cooperative MPC show quicker responses, with the cooperative MPC demonstrating superior super-harmonic response and stabilization time.

To comprehensively compare control strategies, an increasing load was chosen, and current variations are shown in Figure 8A. Temperature control outcomes are presented in Figures 8B,E. Under a constant speed load change, cooperative MPC exhibits minimal temperature fluctuation, effectively tracking current changes and maintaining temperature near the 0.15 K target. By contrast, PID control shows a 0.78 K temperature fluctuation, while water pump MPC surpasses PID by approximately 0.23 K, highlighting MPC’s superior stability. The temperature difference fluctuation under MPC is significantly better than in PID, emphasizing MPC’s enhanced control performance.

Data in Figures 8C,D highlight MPC’s ability to accurately track current changes for efficient heat dissipation and system temperature balance. The PID controlled flow following response is obviously slow, so the reactor outlet temperature under PID control fluctuates longer. This emphasizes MPC’s advantages in achieving rapid and precise thermal management, contributing to improved system performance and stability.

This subsection evaluates the MPC controller’s ability to track targets under changed conditions, using the load current as shown in Figure 9A. At 1,400 s, the outlet control target increases from 343 K to 348 K, and the inlet control target increases from 338 K to 343 K. Figure 9E highlights MPC’s superior control effectiveness in the temperature difference analysis. Temperature control outcomes are presented in Figure 9B.

In Figures 9C,D, the coolant flow rate and heat dissipation air volume respond to the step change in the control target temperature. The initial rapid decline in the coolant flow rate reduces heat removal, leading to an increase in outlet coolant temperature. This results in decreased cooling air volume and an increase in inlet coolant temperature. As the temperature reaches the set value, the controller smoothly transitions into the operational mode. The pump flow surges, and the cooling fan flow increases for efficient heat dissipation. Under MPC, the flow stabilizes rapidly, while under PID control, it gradually decreases before reaching a steady state. After approximately 50 s, the inlet coolant temperature approaches the new set point, and the outlet coolant temperature stabilizes after a brief fluctuation. In summary, MPC effectively controls variable targets.

Considering the substantial impact of ambient temperature on PEMFC, we aimed to compare MPC’s temperature regulation performance across varied ambient temperatures. Figure 10A depicts stack temperature variations under different ambient temperatures. The curve indicates that higher ambient temperatures result in faster reactors reaching the target temperature. Significant differences were observed before 800 s, as highlighted in Figure 10B. The lower temperature curves exhibit smaller oscillations. Our study found a notable similarity between stack inlet temperature and stack temperature performance, which is shown in Figures 10C,D.

Figure 10E compares cooling air flow rate curves as ambient temperature increases from 288.15 K to 308.15 K. Higher ambient temperatures show increased cooling air flow, compensating for reduced heat exchange with the environment. The MPC controller consistently maintains stack inlet temperature within the target range under diverse ambient temperatures. Figure 10F presents a comparison of coolant flow rate curves, where ambient temperature dictates the required coolant flow for the reactor.