Efﬁcient Two-Step Parametrization of a Control-Oriented Zero-Dimensional Polymer Electrolyte Membrane Fuel Cell Model Based on Measured Stack Data

: This paper proposes a new efﬁcient two-step method for parametrizing control-oriented zero-dimensional physical polymer electrolyte membrane fuel cell (PEMFC) models with measured stack data. Parametrizations of these models are computationally intensive due to the numerous unknown parameters and the typically nonlinear, stiff model properties. This work reduces an existing model to decrease its stiffness for accelerated numerical simulations. Subdividing the parametrization into two consecutive subproblems (thermodynamic and electrochemical ones) reduces the solution space signiﬁcantly. A parameter sensitivity analysis further reduces each sub-solution space by excluding non-signiﬁcant parameters. The method results in an efﬁcient parametrization process. The two-step approach minimizes each sub-solution space’s dimension by two-thirds, respectively three-fourths, compared to the global one. An achieved R 2 value between simulation and measurement of 91% on average provides the required accuracy for control-oriented models.


Introduction
PEMFCs are promising candidates for replacing internal combustion engines in mobile applications. However, fuel cell (FC) driven vehicles are still far from having a significant market share because of various challenges. One of the challenges is to control the FC during transient operations, especially regarding avoiding harmful operating conditions. Another one is to improve the FC system's efficiency via optimal control. Moreover, experimental expenditures during development are high but are reducible by using simulations instead. One of the first steps towards resolving the named and other challenges is obtaining a proper FC model. PEMFCs are promising candidates because they offer the rapid startup and low operating temperature required for automotive applications. However, it requires appropriate water management to address liquid water formation. Solid oxide FCs use non-noble metal catalysts, which results in low raw material costs, but the high operating temperature leads to thermal stress and precludes this FC for transient applications. Molten carbonate FCs can reform a wide variety of fuel sources, but the long startup time eliminates this FC for anything but continuous-power applications. Phosphoric acid FCs provide easy water management, but the low power density is not appropriate for portable applications. Alkaline FCs have the highest demonstrated operating efficiency of any FC system, but the intolerance to carbon dioxide forces the use of carbon dioxide removal equipment and pure oxygen and hydrogen [1].
In general, modeling approaches are distinguishable into three groups. First, black-box modeling means fitting an artificial model to replicate the measured input and output data

Fuel Cell Model
This section briefly portrays the used PEMFC model to demonstrate the proposed two-step parametrization method. Additionally, this section describes the model reduction for accelerated numerical simulations. The proposed method is, of course, not limited to the described FC model. Any model separable into a submodel with ODEs and a submodel without ODEs is utilizable. The model does not have to be analytically differentiable because the parameter sensitivity analysis is numerically approximable. However, numerical approximations are less efficient and accurate evaluable than analytical solutions.

Model Description
Ritzberger et al. [10] developed the model, and it is an adapted version of the Pukrushpan et al. [8] model. The difference is that the adapted model has the property of analytical differentiability, which is the reason for its selection. This property is beneficial for various control methodologies [20,21], and parameter sensitivity analyses (Section 3.2). The model is a zero-dimensional physical PEMFC model, and Figure 1 gives a schematic overview. It utilizes, amongst others, mass balances, linear nozzle equations, diffusion equations, electrochemical equations, and the ideal gas law [22][23][24]. However, it does not consider energy balances and thus cannot model the FC temperature over time. The model treats the measured temperature as an input. This temperature is assumed to be the uniform temperature for the whole FC. The model consists of a cathode, anode, membrane, and electrochemical submodel, which are interconnected. The model equations and their derivations are not the focus of this work. More information in this regard is available in [10].
The cathode submodel has four mass states (oxygen mass m ca,O 2 , nitrogen mass m ca,N 2 , vapor mass m ca,vap , and liquid water mass m ca,liq ) and two pressure states (supply manifold pressure p ca,sm , and exit manifold pressure p ca,em ). The anode submodel also has four mass states (hydrogen mass m an,H 2 , nitrogen mass m an,N 2 , vapor mass m an,vap , and liquid water mass m an,liq ) but only one pressure state (exit manifold pressure p an,em ). The membrane submodel only has the membrane water activity state a m . The model inputs are the air mass flowṁ ca,in , the hydrogen mass flowṁ an,in , the anode supply manifold pressure p an,sm , the relative humidity of the air mass flow ϕ ca,in , the purging signal α purge , the FC temperature T, the atmospheric pressure p atm , and the current I. The outputs are the cathode supply manifold pressure p ca,sm , the cathode exit manifold pressure p ca,em , the anode exit manifold pressure p an,em , and the voltage U. Hence, the following equations describe the nonlinear FC state-space model: Here, x nr = x nr (t) denotes the non-reduced state vector, u = u(t) the input vector, y = y(t) the output vector, θ ∈ R 25 the parameter vector, f nr the non-reduced system function, g nr the non-reduced output function, and t the time [25]. The respective vectors are structured as follows: V denotes the volumes, k the nozzle or mass flow coefficients, τ m the water activity time constant, 2 the membrane conductivity parameter, R c the ohmic contact resistance, E Processes 2021, 9, 713 5 of 18 the energy, K the intrinsic exchange current parameter, and CD the combined diffusion coefficient. Compared to the model described by Ritzberger et al. [10], this one has an additional output, the cathode exit manifold pressure p ca,em . It is measured, therefore considering it as a supplementary output has the advantage that it yields additional insight into the system. Furthermore, the pressure difference between the exit manifold and the environment is too big to be sufficiently described by a linear nozzle equation. Hence nonlinear nozzle equations (derived from the Bernoulli equation) replace the linear ones at each exit manifold.

Model Reduction
The model reduction goal in this work is to reduce the model's stiffness for accelerated numerical simulations. Lambert [26] stated that stiffness occurs when (numerical) stability requirements, rather than those of accuracy, constrain the (simulation time) step length. The pressure dynamics of the model are faster by multiple orders of magnitude compared to the mass dynamics. Thus the pressure dynamics are assumed to be steady-state at all times. The ODE for the cathode supply manifold pressure p ca,sm results from a combination of the ideal gas law and the mass balancė where R ca,sm denotes the cathode supply manifold's mass-specific gas constant,ṁ ca,sm,cm = k ca,sm,out (p ca,sm − p ca,cm ) the mass flow between the cathode supply and center manifold, k ca,sm,out the cathode supply manifold's outflow nozzle coefficient, and p ca,cm the cathode center manifold pressure. Assuming steady-state, the ODE (4) transforms into p ca,sm =ṁ ca,in k ca,sm,out + p ca,cm .
The exit manifold pressures of the cathode p ca,em and anode p an,em are addressed in a similar fashion leading to the reduced nonlinear FC state-space model given bẏ Here, f denotes the reduced system function, g the reduced output function, and x = x(t) the reduced state vector. The reduced state vector is structured as x = [m ca,O 2 , m ca,N 2 , m ca,vap , m ca,liq , m an,H 2 , m an,N 2 , m an,vap , m an,liq , a m ] T , and it does not contain the three pressure states (p ca,sm , p ca,em , and p an,em ) anymore.The input vector u, the output vector y, and the parameter vector θ remain unchanged. The longest integration step length of the reduced model, which still yields a stable solution, is about 50% longer than for the non-reduced model, leading to a roughly 36% shorter simulation time [27].

Key Idea
The two-step parametrization method's key idea is to subdivide the given nested state-space model (6) and (7) into two submodels. The nested model is separable in the following way: The model's parameter vector is divisible as well: On the one hand, the subdivision yields the thermodynamic submodel and on the other hand, it results in the electrochemical submodel The thermodynamic submodel (12) and (13) is only dependent on the thermodynamic parameter vector θ th , which also holds for the system function f according to Equation (9). Therefore the electrochemical parameter vector θ el does not affect the states x and the thermodynamic outputs y th , but only the electrochemical output y el . Thus the electrochemical submodel (14) utilizes the full parameter vector θ = [θ th , θ el ] T . Note that the thermodynamic submodel consists of ODEs, and the electrochemical one does not. The two-step parametrization method exploits the described properties.

Parameter Sensitivity Analysis
The FC model used in this work has numerous unknown parameters, which raises the question of whether each parameter is unambiguously identifiable. The parameter identifiability is, in general, strongly dependent on the model structure and the available measurement data. A parameter sensitivity analysis can aid in answering these questions, and the Fisher information matrix (FIM) F is well-established for conducting such analysis [28]. Under the assumption of Gaussian prediction errors with zero mean values and time-independent covariances, the Cramér-Rao inequality holds [29]: The inequality says that the inverse of F is the lower bound of the parameter covariances. The first step for obtaining the FIM is computing the state parameter sensitivities . . , n θ } denotes a parameter, and n θ the number of parameters. They are obtainable from solving the following first-order ODE: The second step is to calculate the output parameter sensitivities ψ i = dy/dθ i with The model described by Ritzberger et al. [10] has analytic derivatives available, and this work utilizes MATLAB R2020b's Symbolic Math Toolbox [30] to compute them. All output parameter sensitivities ψ i merged into one matrix yields the output parameter sensitivity matrix ψ(t) = [ψ 1 (t), ψ 2 (t), . . . , ψ n θ (t)]. The FIM is determinable with ψ(t) taken at the sample times t k for k ∈ {0, 1, . . . , n k }, where n k + 1 is the number of sample instants. Finally, using the sampled output parameter sensitivity matrix ψ(t k ), the FIM is computable with where Σ e denotes the prediction error covariance matrix, and F ∈ R n θ ×n θ holds. Under the assumption of a perfect model, the prediction error covariance Σ e is identical to the measurement noise covariance and thus assumed to be known. According to the Crámer-Rao inequality (15), the inverse of the FIM is the lower bound of the absolute parameter variances. Thus directly analyzing F would mean that the parameter's physical units and their magnitudes bias the analysis. Using the normalized dimensionless FIM F norm for a physical unit independent analysis resolves this issue [31][32][33]: The identifiability of the parameters is theoretically derivable from the inverse of F norm . Unfortunately, this is not always feasible because the FIM is often ill-conditioned or even singular. Using a singular value decomposition (SVD) of F norm resolves this issue and enables further analysis [32][33][34]: U denotes the matrix containing the left singular vectors, Σ = diag(σ 1 , σ 2 , . . . , σ n θ ) the singular value matrix, and V the matrix containing the right singular vectors. The left and the right singular vector matrix are identical because F norm is a symmetric matrix. The analysis utilizes the right singular vector matrix V = [v 1 , v 2 , . . . , v n θ ], where v l for l ∈ {1, 2, . . . , n θ } denotes a singular vector corresponding to the singular value σ l . The singular values are interpretable as the amount of information, and the corresponding singular vectors determine the direction in the parameter space. The Euclidean norm of the vectors is 1. Thus the relative direction share of the vector component v l,i is v 2 l,i [35]. The relative direction v 2 l,i multiplied with the singular value σ l is the amount of information showing into the parameter θ i 's direction. Summing up all the information of one parameter from each singular value yields the parameter's total information The parameter's total information σ θ i indicates if a parameter is identifiable with the given model and measurement data. The most significant parameters θ ms are determinable by sorting the parameters according to their total information in descending order. The most significant parameters are determined by summing up the first n θ ms parameters until the following inequality holds: Note that γ ∈ [0, 1] denotes the threshold and is adaptable for different purposes. This work uses γ = 0.99999. In this case, the most significant parameters θ ms describe ≥99.999% of F norm , and the least significant ones θ ls describe ≤0.001%, which makes the latter negligible for parametrization.

Procedure
The proposed method parametrizes the model in two consecutive main steps:

1.
Thermodynamic submodel (a) Parameter sensitivity analysis with respect to the thermodynamic parameters θ th yields a subset with the most significant parameters θ th,ms , where θ th,ms ⊆ θ th holds.
Parametrization with respect to the most significant parameters θ th,ms yields the optimized parameters θ th,opt . The least significant parameters θ th,ls are kept constant at their initial values.

2.
Electrochemical submodel (a) Solve thermodynamic submodel using the optimized thermodynamic parameters θ th,opt and store the resulting model states x for further usage.
Parameter sensitivity analysis with respect to the electrochemical parameters θ el yields a subset with the most significant parameters θ el,ms , where θ el,ms ⊆ θ el holds. (c) Parametrization with respect to the most significant parameters θ el,ms yields the optimized parameters θ el,opt . The least significant parameters θ el,ls are kept constant at their initial values.
Combining both solutions lead to the full optimized parameter vector Assuming that the result is near the optima, the entire model with its optimized most significant parameters θ ms,opt = θ th,ms,opt , θ el,ms,opt T can be further refined with iterative methods to approach the optima while the least significant parameters are kept constant θ ls = [θ th,ls , θ el,ls ] T at their initial values. The two-step parametrization method requires measurement data (according to the model's inputs and outputs), and a model separated into a submodel with ODEs and a submodel without it. The optimization goal is to minimize an objective function J, which contains the weighted squared errors between the simulated and measured output signals and a regularization term [19]: Here, y(t k , θ) denotes the model output at sampling instant t k for k ∈ {1, 2, . . . , n k }, y * (t k ) the measured output, Q y the output weighting matrix, θ 0 a plausible initial guess of the parameter vector, and Q θ the regularization matrix. The weighting matrix Q y takes the individual weighting of each output and the different output magnitudes into account. The last term in the objective function J (23), also called regularization, penalizes the deviation of the parameter vector θ from its initial guess θ 0 and takes the different parameter magnitudes into account. The regularization is also beneficial if the parameters' uncertainty differs, for example, if the approximate values for some parameters are derivable from literature. The goal is to minimize the objective function J. In this case, the optimization problem is stated as follows: Solving the optimization problem yields the optimized parameter vector θ opt . For optimal results, the parameter space has to be constrained. The parameter bounds, θ i,min and θ i,max , are derivable from physical considerations and expert knowledge.
First, parametrizing the thermodynamic submodel with its parameter vector θ th reduces the solution space's dimension by 8 25 = 32% compared to the global one. With the determined optimized thermodynamic parameter vector θ th,opt and a given input u, the state trajectories x do not change anymore. Therefore, for the electrochemical submodel's parametrization, the thermodynamic submodel only needs to be solved once, and the resulting states x are stored for further usage. Second, parametrizing the electrochemical submodel with its parameter vector θ el reduces the solution space's dimension by 17 25 = 68% compared to the global one. Here only the electrochemical parameters are optimized, and the thermodynamic ones are kept constant. A parameter sensitivity analysis further reduces each sub-solution space's dimension by excluding non-significant parameters, see Section 3.4. This proceeding significantly simplifies the optimization problem. The thermodynamic submodel's parametrization solves ODEs in every iteration, and the electrochemical submodel's parametrization does not, which makes the latter considerably faster.

Validation of Method
This section shows the validation of the two-step parametrization method. Simulating the model with a plausible initial guess of the parameter vector θ 0 and using the measured inputs u (described in Section 4 and depicted in Section 5 under parametrization data) yields the simulated model outputs y. The initial values x(t 0 ) are obtainable by assuming steady-state at the first sample instant. The data sequence is appropriately selected so that at sampling instant t 0 , the steady-state assumption holds. The simulated model outputs with additional Gaussian noise replace the measured outputs y * as reference data to validate the method. In this case, the proposed method has to deliver the used parameter θ 0 on average if it is an unbiased estimator.
Conducting the two-step parametrization method described in Section 3 yields the parameters' total information, shown in Figure 2a,b. Due to the floating-point arithmetic computation, the result's accuracy depends on the used hardware [36]. Therefore, the parameters with a total information value less than the hardware accuracy (indicated by the yellow line) are not identifiable. The most significant thermodynamic parameters in descending order are k ca,em,out , k ca,sm,out , k ca,em,in , k an,sm,out , V ca,cm , k cond , k an,em,in , k perm , k an,em,out , and V an,cm . The least significant ones in ascending order are V ca,sm , V an,em , V an,sm , V ca,em , k evap , k an,leak , and τ m . Only parametrizing the most significant thermodynamic parameters reduces the dimension of the sub-solution space by 15 25 = 60% compared to the global one. The most significant electrochemical parameters in descending order are E ca,act , E an,act , 2 , K ca , K an , and R c . The least significant ones in ascending order are CD an , and CD ca . Only parametrizing the most significant electrochemical parameters reduces the dimension of the sub-solution space by 19 25 = 76% compared to the global one. The least significant parameters are kept constant at their initial values at all times. fore, the parameters with a total information value less than the hardware accuracy (in-dicated by the yellow line) are not identifiable. The most significant thermodynamic parameters in descending order are k ca,em,out , k ca,sm,out , k ca,em,in , k an,sm,out , V ca,cm , k cond , k an,em,in , k perm , k an,em,out , and V an,cm . The least significant ones in ascending order are V ca,sm , V an,em , V an,sm , V ca,em , k evap , k an,leak , and τ m . Only parametrizing the most significant thermodynamic parameters reduces the dimension of the sub-solution space by 15 25 = 60% compared to the global one. The most significant electrochemical parameters in descending order are E ca,act , E an,act , 2 , K ca , K an , and R c . The least significant ones in ascending order are CD an , and CD ca . Only parametrizing the most significant electrochemical parameters reduces the dimension of the sub-solution space by 19 25 = 76% compared to the global one. The least significant parameters are kept constant at their initial values at all times. parameters in descending order are k ca,em,out , k ca,sm,out , k ca,em,in , k an,sm,out , V ca,cm , k cond , k an,em,in , k perm , k an,em,out , and V an,cm . The least significant ones in ascending order are V ca,sm , V an,em , V an,sm , V ca,em , k evap , k an,leak , and τ m . Only parametrizing the most significant thermodynamic parameters reduces the dimension of the sub-solution space by 15 25 = 60% compared to the global one. The most significant electrochemical parameters in descending order are E ca,act , E an,act , ǫ 2 , K ca , K an , and R c . The least significant ones in ascending order are CD an , and CD ca . Only parametrizing the most significant electrochemical parameters reduces the dimension of the sub-solution space by 19 25 = 76% compared to the global one. The least significant parameters are kept constant at their initial values at all times. The two submodels with their most significant parameters are subject to parametrization. MATLAB R2020b's genetic algorithm optimizer parametrizes each submodel with respect to the most significant parameters θ ms 100 times [37]. Heuristic algorithms increase the probability of finding the global minimum in a high dimensional solution space compared to iterative methods. Every iteration randomly initializes the genetic algorithm optimizer and the Gaussian noise for the model outputs. The optimizer's population size is 200, the generation limit is 20, and the remaining options remain unchanged. Figure 2c,d illustrate the relative distribution of the estimated parameters θ ms,opt . The conclusion is that the estimated parameters are equal to the used parameters θ 0 on average, validating the proposed method. The respective spreads of the estimates correlate with the total information. The more information exists, the smaller the spread. The correlation is not Total information of thermodynamic σ θ th,i (obtained from F th,norm ) and electrochemical parameters σ θ el,i (obtained from F el,norm ), respectively. Plots (c,d): Relative distribution (100 independent estimations) of optimized most significant parameters of thermodynamic θ th,ms,opt,i and electrochemical submodel θ el,ms,opt,i , respectively. Plot (e): Relative distribution (100 independent estimations) of optimized parameters of entire model θ opt,i with usual "single-step" approach. θ 0,i denotes true parameter value, and red + symbol indicates outliers in the box plots.
The two submodels with their most significant parameters are subject to parametrization. MATLAB R2020b's genetic algorithm optimizer parametrizes each submodel with respect to the most significant parameters θ ms 100 times [37]. Heuristic algorithms increase the probability of finding the global minimum in a high dimensional solution space compared to iterative methods. Every iteration randomly initializes the genetic algorithm optimizer and the Gaussian noise for the model outputs. The optimizer's population size is 200, the generation limit is 20, and the remaining options remain unchanged. Figure 2c,d illustrate the relative distribution of the estimated parameters θ ms,opt . The conclusion is that the estimated parameters are equal to the used parameters θ 0 on average, validating the proposed method. The respective spreads of the estimates correlate with the total information. The more information exists, the smaller the spread. The correlation is not Total information of thermodynamic σ θ th,i (obtained from F th,norm ) and electrochemical parameters σ θ el,i (obtained from F el,norm ), respectively. Plots (c,d): Relative distribution (100 independent estimations) of optimized most significant parameters of thermodynamic θ th,ms,opt,i and electrochemical submodel θ el,ms,opt,i , respectively. Plot (e): Relative distribution (100 independent estimations) of optimized parameters of entire model θ opt,i with usual "single-step" approach. θ 0,i denotes true parameter value, and red + symbol indicates outliers in the box plots.
The two submodels with their most significant parameters are subject to parametrization. MATLAB R2020b's genetic algorithm optimizer parametrizes each submodel with respect to the most significant parameters θ ms 100 times [37]. Heuristic algorithms increase the probability of finding the global minimum in a high dimensional solution space compared to iterative methods. Every iteration randomly initializes the genetic algorithm optimizer and the Gaussian noise for the model outputs. The optimizer's population size is 200, the generation limit is 20, and the remaining options remain unchanged. Figure 2c,d illustrate the relative distribution of the estimated parameters θ ms,opt . The conclusion is that the estimated parameters are equal to the used parameters θ 0 on average, validating the proposed method. The respective spreads of the estimates correlate with the total information. The more information exists, the smaller the spread. The correlation is not perfect but will converge to the Cramer-Ráo bound (15) by increasing the population size, the generation limit, and the number of estimations.
For comparison, Figure 2e depicts the relative distribution of the estimated parameters θ opt,i with a usual "single-step" approach for the entire model. The genetic algorithm optimizer minimizes the objective function J (23) in a single step for the entire model without any further considerations. The resulting spreads of the parameters are much higher than those from the two-step method, which means a more complicated search for the optima.

Experimental Setup
A PEMFC system test bench was set up for the acquisition of the measurement data, which are required for the parametrization. In Figure 3, the schematic overview of the test bench is shown. The main components of the test bench are a 30 kW PEMFC stack, a hydrogen and an air supply system, the cooling circuits as well as the control system. The experimental tests were conducted by setting the load point of the FC stack via a dynamic DC/AC inverter (battery simulator). Either the current or the voltage level can be controlled by the battery simulator. Further, the battery simulator supplies the electric power, which is generated by the FC stack during operation, into the electric grid. Additionally, a power supply was designed to feed the balance of plant components (e.g., air compressor) with energy from the electric grid.

Media Supply
The reactants of the stack are fed by the hydrogen and air supply system. The hydrogen of high purity (99.999%) is supplied by a high pressure (300 bar) gas cylinder bundle. In the first stage, the highly compressed hydrogen is reduced to a medium pressure level of 6 bar by means of a pressure reducer. Subsequently, before the hydrogen enters the anode of the stack, the pressure controller further decompresses the hydrogen and ensures a constant pressure inside the anode. The anode of the stack operates in the flow-through mode. Therefore, a hydrogen recirculation pump is used to guarantee a constant volume flow through the anode. Besides, the stack operates without any device for humidification, either of hydrogen or air. Only the internal humidification of the reactants inside the stack is utilized. Complementary to the hydrogen purge valve, a water separator is installed at the anode side to remove the excess water periodically from the hydrogen gas. The air, which is required to supply the necessary oxygen for the electrochemical reaction, is taken in from the conditioned test bench room. For all experimental tests in this paper, the room is operated at a constant temperature of 23°C and constant relative air humidity of 50%. The turbo compressor sucks in the ambient air from the room via an air filter, which is designed to mechanically and chemically clean the intake air from impurities. After the air compressor, the air flows through an intercooler to the cathode side of the FC stack. In order to vary the backpressure of the stack, an electronically controlled air throttle valve is implemented at the process air exhaust of the FC stack.

Cooling Circuits
Two different cooling circuits of the test bench environment are used to keep the stack, the air compressor, and the air intercooler at an appropriate temperature. All the other components of the system are designed so that passive cooling is sufficient for them. The cooling circuit of the FC is thermally connected via a heat exchanger to the 6°C cooling circuit of the test bench environment. Moreover, in the cooling circuit of the FC stack, a deionized coolant is used. During the operation of the FC stack, the values of the volume flow and the speed of the coolant pump are maintained constant. The 40°C cooling circuit of the test bench environment is utilized to cool the air compressor and the intercooler. In order to guarantee a sufficient cooling of the components, flow control valves in the cooling circuits of the test bench environment are integrated. The target stack coolant inlet temperature is 55°C. The air inlet temperature is kept at a value of 40°C. The air compressor is supplied at all times with a constant volume flow of coolant, which is sufficient for the necessary cooling demand.

Test Bench Control System
Dedicated software was developed in LabVIEW and implemented on a NI Com-pactRio to monitor and control the FC system. With this specific software, the possibility has been created to vary each operating parameter of the FC system in an admissible range. Furthermore, with the control system the acquisition of measurement data is realized with a sampling rate of 10 Hz. In order to provide a detailed experimental investigation, a variety of sensors (temperature, pressure, mass flow, voltage, current, and humidity) are installed at all relevant positions. Further details with respect to the test bench set up are given in [38,39].

Experimental Tests and Operating Conditions
During the experimental tests, the load point of the stack and the air mass flow, which is controlled according to the electric current and to a constant air excess ratio of 1.5, were varied. In order to obtain measurement data under different operating states as well as operating conditions of the FC system, the load point was adjusted arbitrarily either by setting the stack voltage or the stack current level in a range, in which a stable operation of the FC is still guaranteed. All other operating parameters of the FC system (target temperatures, anode pressure setpoint, hydrogen recirculation pump speed, coolant pump speed) were kept constant under all load and operating conditions. The standard values of the constant operating parameters can be seen in Table 1. For more details regarding the experimental tests, please refer to Section 5, in which simulated and measured data are presented and discussed.

Results
This work proposes a validated two-step parametrization method, and it demonstrates the method by parametrizing the presented model using the measurement data obtained from the FC test bench. The proceeding is almost identical to the one described in Section 3.4. The differences are that the measured outputs y * now serve as reference data, the optimizer's population size is 1000, and the model gets parametrized only once. Therefore all conclusions stated there still hold. Confidentiality agreements do not allow the publication of the numerical values of the optimized parameters. Figure 4a-d show the parametrization results. Simulating the model with the optimized parameter vector θ opt and using the measured inputs u described in Section 4 yields the simulated model outputs y. The figures additionally depict the measured outputs and four model inputs as reference: the mass flow of airṁ ca,in , the FC temperature T, the anode supply manifold pressure p an,sm , and the current I. The parameterized model achieves an R 2 value of 93% on average. For validation data, the model reaches an R 2 value of 91% on average (see Figure 4e-h), providing the required accuracy for control-oriented applications.

Discussion
The simulated outputs fit very well. However, the simulated pressures deviate more from the measurements than the voltage. The simulated anode exit manifold pressure p an,em (Figure 4c,g) has deviations, which seem to be load-dependent. During low currents, the measured pressure is higher and vice versa. A reason could be the recirculation flow. The recirculation pump runs at a constant speed and should provide a constant flow during steady-state. The flow is not measured and is assumed constant in the model. However, steady-state conditions are not always given, especially during purging. Measuring the recirculation flow and adding it as an input to the model could minimize the deviations. The simulated cathode pressures deviate more from the validation data (Figure 4e,f) than from the parametrization data (Figure 4a,b). The maximum deviation is 5%, and a reason may be the inflow air temperature, which differs up to 10°C between the two data sets. Incorporating the inflow air temperature into the model may improve the agreement between the simulated and measured cathode pressures. The model's simulated voltage response (Figure 4h) replicates the undershooting only in a moderate way for the validation data. One reason could be the underrepresented voltage undershooting in the parametrization data (Figure 4d). Using a data set with more often recurring voltage undershooting may improve the model's behavior in this regard. Figure 2b depicts the total information of the electrochemical parameters σ θ el,i . The combined diffusion coefficients, CD ca and CD an , have relatively low information. The reason is that the experiments only cover the safe operating region. The named coefficients determine the limit current and are only well identifiable in the unsafe limit current region. According to Figure 2a, the supply and exit manifold volumes (V ca,sm , V ca,em , V an,sm , and V an,em ) are not identifiable. The reason is that only the pressure states needed them, which do not exist anymore. Thus the named volumes are not used in the reduced model at all. The named aspects only concern the model and the experiments, but not the parametrization method itself.
A drawback of the two-step method is that by conducting separate parameter sensitivity analyses for the two submodels, slightly different conclusions may follow compared to the complete model analysis. The reason is that the electrochemical model's output contributes additional information on the thermodynamic parameters. However, the advantage of the reduced solution space dimension compensates for this drawback. The stated conclusions derived from the FIM only hold in a region near the initial guess of the parameters θ 0 . Hence the FIM is only a local parameter sensitivity analysis. With another parameter vector, the conclusions from the newly derived FIM may be different. In general, parameters could be nonsignificant if the excitation is ill-suited (e.g., CD ca ), it is not relevant in the model (e.g., V ca,sm ), or both. High significancies are interpretable in a vice versa manner. As discussed, the parameter identifiability and the derived conclusions are strongly dependent on the available measurement data and the model structure. With this knowledge, remodeling or a more targeted design of the experiment is possible to improve identifiability. From the engineering point of view, this information is utilizable to distribute engineering resources efficiently. Focusing on improving the significant parameter's real-world counterparts affects the system outputs substantially more than optimizing the less significant ones. Optimizing the latter will lead to hardly any changes in the output. Furthermore, more measured signals (e.g., anode and cathode outflow relative humidity, species concentrations, recirculation flow) and more "exciting" system inputs would significantly affect the obtainable results by increasing the Fisher information of the parameters, leading to better parametrization results.

Discussion
The simulated outputs fit very well. However, the simulated pressures deviate more from the measurements than the voltage. The simulated anode exit manifold pressure p an,em (Figure 4c,g) has deviations, which seem to be load-dependent. During low currents, the measured pressure is higher and vice versa. A reason could be the recirculation flow. The recirculation pump runs at a constant speed and should provide a constant flow during steady-state. The flow is not measured and is assumed constant in the model. However, steady-state conditions are not always given, especially during purging. Measuring the recirculation flow and adding it as an input to the model could minimize the deviations. The simulated cathode pressures deviate more from the validation data (Fig-Figure 4. Plots (a,e): Cathode supply manifold pressure p ca,sm (blue, yellow and red), and air mass flowṁ ca,in (purple). Plots (b,f): Cathode exit manifold pressure p ca,em (blue, yellow and red), and FC temperature T (purple). Plots (c,g): Anode exit manifold pressure p an,em (blue, yellow and red), and anode supply manifold pressure p an,sm (purple). Plots (d,h): Stack voltage U (blue, yellow and red), and stack current I (purple). The plots depict parameterization and validation data, respectively, and y denotes output and u input.

Discussion
The simulated outputs fit very well. However, the simulated pressures deviate more from the measurements than the voltage. The simulated anode exit manifold pressure p an,em (Figure 4c,g) has deviations, which seem to be load-dependent. During low currents, the measured pressure is higher and vice versa. A reason could be the recirculation flow. The recirculation pump runs at a constant speed and should provide a constant flow during steady-state. The flow is not measured and is assumed constant in the model. However, steady-state conditions are not always given, especially during purging. Measuring the recirculation flow and adding it as an input to the model could minimize the deviations. . Plots (a,e): Cathode supply manifold pressure p ca,sm (blue, yellow and red), and air mass flowṁ ca,in (purple). Plots (b,f): Cathode exit manifold pressure p ca,em (blue, yellow and red), and FC temperature T (purple). Plots (c,g): Anode exit manifold pressure p an,em (blue, yellow and red), and anode supply manifold pressure p an,sm (purple). Plots (d,h): Stack voltage U (blue, yellow and red), and stack current I (purple). The plots depict parameterization and validation data, respectively, and y denotes output and u input.
Compared to the usual "single-step" approach, the two-step method simplifies the parametrization process significantly by reducing the solution space's dimension. As depicted in Figure 2e, the parameters resulting from a single-step parametrization of the entire model (without any further considerations) have a much larger spread, which implicitly means a worse fit of the simulation data with the measurements. Figure 4 additionally shows the simulation results with parameters obtained from the single-step approach. The simulations achieve an average R 2 value of 33% for the parametrization data and −62% for the validation data, which is much worse than the proposed two-step method (93% and 91%, respectively). However, increasing the optimizer's population size and the number of generations the single-step approach would yield similar results as the two-step method. By doing so, the single-step approach searches more thoroughly through the bigger solution space of the entire model for the optima, but at the expense of exponentially growing computation time. Therefore, considering limited resources and the same optimizer options, the two-step method delivers more accurate results than the single-step approach, making the former superior in this case. A comparison of the proposed two-step method with methods consisting of more than two steps is, in this case, not meaningful because this work utilizes pre-existing data. Multiple-step methods [14,18] require appropriately designed experiments, which is usually not the case for pre-existing data. If the data allows the utilization of a method with more than two steps, then the multiple-step method will likely deliver better results than the two-step method considering similar resources. The fact that the two-step method does not need appropriately designed experiments compensates for this drawback.

Conclusions
This paper presents an efficient two-step parametrization method for FCs. A separation of the model into two submodels decreases the solution space drastically. The method parametrizes the submodels in two consecutive steps. A parameter sensitivity analysis further reduces each sub-solution space, which simplifies the search for the optima. This work demonstrates the efficient method by parametrizing an FC model with measurements and illustrates the parameter sensitivity analysis's advantages. The parameterized model's outputs replicate the validation data excellently.
The parameter sensitivity analysis is a powerful tool to determine the parameters' identifiability with a given model and measurement data. The inverse approach is also promising for future research: developing a "design of experiment controller" with a given model to maximize the parameter's identifiability in real-time.
Author Contributions: Conceptualization, and methodology, Z.P.D. and S.J.; investigation, data curation, Z.P.D. and C.S.; software and visualization, Z.P.D.; writing-original draft preparation, Z.P.D. and C.S.; writing-review and editing, supervision, and project administration, S.J. All authors have read and agreed to the published version of the manuscript.

Abbreviations
The following abbreviations are used in this manuscript:

FC
Fuel cell FIM Fisher information matrix ODE Ordinary differential equation PEMFC Polymer electrolyte membrane fuel cell SVD Singular value decomposition

Nomenclature
The following symbols are used in this manuscript: