Effective Parameterization of PEM Fuel Cell Models—Part II: Robust Parameter Subset Selection, Robust Optimal Experimental Design, and Multi-Step Parameter Identification Algorithm

In this work, we use a physics-based model that has been developed recently for automotive applications.^37–39 Particularly, the model is pseudo-2D and captures non-isothermal and two-phase phenomena and has been extensively validated with performance and resistance data obtained experimentally from state-of-the-art fuel cell stacks.³⁹ Further details about the model and its validation can be found in the original publications.^37–39

The sensitivities of several of the model predictions to a variety of parameters were examined in the first part of this work using an extended local analysis.⁴⁰ Particularly, cell voltage, high frequency resistance (HFR), and normalized membrane water crossover were studied for their sensitivity to 50 of the model parameters. The local analysis was conducted about several points in the parameter space to obtain a more global picture of the model sensitivities. Additionally, the sensitivities were generated using a large library of 729 different operating conditions with variations in temperature, pressure, relative humidity, and stoichiometric ratio. The parameters considered included geometrical and structural features of the cell, along with a number of fitting coefficients. The results were analyzed for sensitivity magnitudes and their relative independence, as they have implications for parameter identifiability. It was also demonstrated that many of the model parameters are not identifiable from cell-level measurements. This unidentifiability is often due to insensitivity of a parameter or multiple parameters having the same impact on the predicted outputs, which points to opportunities for model refinement, e.g., by re-parameterizing the model.^41–43 It was further shown how additional measurements can enhance identifiability of model parameters by increasing their corresponding impact on the predicted outputs and reducing the correlation between the effect of different parameters (see the first part of the study for further discussion⁴⁰).

Having analyzed model sensitivities, our exclusive focus in this part of the study is on parameter identifiability. Accordingly, we note that while we had considered 50 parameters for sensitivity analysis, many of these parameters, such as those describing the channel geometry, are known for a given stack. Although these parameters were included in our sensitivity analysis to ensure that our results are not constrained to a particular cell design, they are not suitable candidates for identification. Therefore, we narrow down the list of considered parameters in this part to focus on those that may be considered for model parameterization. Specifically, we limit our analysis here to a set of 31 parameters that are given in the Appendix Table A·I, which also provides the bounds for each parameter. Note that the specified bounds, which are determined based on the uncertainty of the corresponding parameters reported in the literature, ensure the parameters maintain their physical significance. Regarding the outputs of interest, here we assume that voltage and HFR are measured in the experiments as they often are in real-world testings. Therefore, our discussion here is based on the combined sensitivity of the voltage and HFR predictions. A summary of the sensitivity results is shown in Fig. 1 for the 31 parameters that are subject to this part of the study.

Zoom In Zoom Out Reset image size

Attempting to identify all of the model parameters using only voltage and HFR data results in an ill-conditioned inverse problem.⁴⁴ This necessitates the need for parameter subset selection so that the identification algorithm can focus exclusively on the parameters that are practically identifiable with the available measurements. We recall that a parameter is practically identifiable when it has an observable impact on the output predictions that is uniquely distinguishable from the impact of other parameters.¹⁷ Therefore, a large sensitivity magnitude alone does not guarantee identifiability of a parameter; the corresponding sensitivity vector must also be independent from that of other parameters. Accordingly, this section provides an algorithm that selects an identifiable subset of parameters. This subset selection procedure may also be viewed as a regularization technique to alleviate the ill-conditioning of the inverse problem.^45,46 Briefly, the proposed method is as follows:

• 1.  
  Inspect the parameter set for highly collinear pairs and remove the parameter with smaller sensitivity magnitude from such pairs.
• 2.  
  Choose the number of parameters, n_I, to be included in the subset. This is done using singular value decomposition.
• 3.  
  Find the locally optimal subset of parameters at each point in the parameter space, for which sensitivity information has been calculated.
• 4.  
  Find the globally optimal subset of parameters for the entire parameter space. The local solutions from the previous step are used to guide the optimization.

Below, each step is described in detail.

We start by inspecting the parameter sensitivities for highly collinear pairs. Given two sensitivity vectors ${{\boldsymbol{s}}}_{i}^{k}$ and ${{\boldsymbol{s}}}_{j}^{k}$ for parameters i and j obtained at point k in the parameter space, their collinearity index is defined as:

$$\begin{eqnarray}&&{\psi }_{i,j}^{k}=\cos \left({\phi }_{i,j}^{k}\right)=\displaystyle \frac{\left|{{{\boldsymbol{s}}}_{i}^{k}}^{{\mathsf{T}}}{{\boldsymbol{s}}}_{j}^{k}\right|}{{\parallel {{\boldsymbol{s}}}_{i}^{k}\parallel }_{2}{\parallel {{\boldsymbol{s}}}_{j}^{k}\parallel }_{2}},\end{eqnarray} \tag{ 1 }$$

where ${\phi }_{i,j}^{k}$ is the angle between the sensitivity vectors and the collinearity index ${\psi }_{i,j}^{k}$ varies between 0 and 1, indicating orthogonal and collinear sensitivity vectors, respectively. This collinearity index is calculated for all parameter pairs at all points in the parameter space, for which we had generated sensitivity data in the first part of the study.⁴⁰ The median of the calculated indices is then used as the representative collinearity index for each parameter pair over the entire parameter space:

$$\begin{eqnarray}&&{\overline{\psi }}_{i,j}=\mathop{\mathrm{med}}\limits_{k}\,{\psi }_{i,j}^{k},\end{eqnarray} \tag{ 2 }$$

where k is the index for the sample point in the parameter space. Only one parameter from a pair whose median collinearity index is high can be practically identified.²⁰ Therefore, for pairs with ${\overline{\psi }}_{i,j}$ exceeding a threshold value of ψ_th, the parameter with lower sensitivity magnitude is discarded. Here we use ψ_th = 0.9, which results in removal of i_0,ca and E_rev,PtO from the list of parameters due to their significant collinearity with α_ca and ω, respectively. We also note that a large threshold is chosen to ensure a certain degree of flexibility is maintained in the resulting set of parameters. More specifically, smaller thresholds can lead to removal of parameters that can have a noticeably different impact than their collinear counterparts only at a few selective operating conditions.

In the second step we determine the number of parameters, n_I, that are needed to capture the variations in the data due to perturbations to the remaining parameters. Singular value decomposition⁴⁷ is a powerful tool for this purpose. In particular, at each sample point in the parameter space we calculate the singular values of the sensitivity matrix. Let these singular values be ordered as follows:

$$\begin{eqnarray}&&{\sigma }_{1}^{k}\gt {\sigma }_{2}^{k}\ \gt ...\gt \ {\sigma }_{{n}_{\theta }}^{k},\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{i}^{k}$ denotes the i-th largest singular value of the sensitivity matrix obtained at the k-th sample point in the parameter space and n_θ is the number of parameters under consideration after the lower sensitivity parameters in collinear pairs have been removed (n_θ = 29). We then seek to find the smallest r such that:

$$\begin{eqnarray}&&{\overline{\sigma }}_{r}^{k}=\displaystyle \frac{{\sum }_{j=1}^{r}{\sigma }_{j}^{k}}{{\sum }_{i=1}^{{n}_{\theta }}{\sigma }_{i}^{k}}\geqslant {\overline{\sigma }}_{\mathrm{th}},\end{eqnarray} \tag{ 4 }$$

where ${\overline{\sigma }}_{r}^{k}$ monotonically increases with r and provides a qualitative measure as to how much of the total variations in the sensitivity directions are captured by the first r singular vectors. Moreover, ${\overline{\sigma }}_{\mathrm{th}}$ is a threshold value that determines the number of selected singular values. Choosing a threshold of ${\overline{\sigma }}_{\mathrm{th}}=0.9$, we have found that, on average, 12 singular values are required to satisfy the above inequality at the sample points in the parameter space. Accordingly, we focus on selecting a subset of size n_I = 12 for identification.

The third step of the algorithm is finding the locally optimal subset of parameters for identification. To this end, we must determine which of the model parameters have the largest unique impact on the predicted outcomes. This is done using the sensitivity matrix given by:

$$\begin{eqnarray}&&{{\bf{S}}}^{k}={\left.\displaystyle \frac{\partial \tilde{{\boldsymbol{y}}}}{\partial {\boldsymbol{\theta }}}\right|}_{{\boldsymbol{\theta }}={{\boldsymbol{\theta }}}^{k}},\end{eqnarray} \tag{ 5 }$$

where $\tilde{{\boldsymbol{y}}}$ is the vector of output predictions and ${\boldsymbol{\theta }}$ is the vector of parameters and the sensitivity matrix is evaluated at the k-th sample point in the parameter space. While our sensitivity calculations are slightly different,⁴⁰ the structure of the resulting matrix is the same. To select a parameter subset, we want to choose columns of the sensitivity matrix that ensure the corresponding parameters can be identified with the required precision. The parameter estimation precision is determined by the corresponding covariance (${{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}$), which is constrained by the inverse of the FIM through the Cramer-Rao bound⁴⁸:

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{\boldsymbol{\theta }}}\geqslant {({{\bf{F}}}^{k})}^{-1},\end{eqnarray} \tag{ 6 }$$

where F^k is the FIM at the k-th sample point in the parameter space, which is given by:

$$\begin{eqnarray}&&{{\bf{F}}}^{k}={{{\bf{S}}}^{k}}^{{\mathsf{T}}}{{\bf{S}}}^{k}.\end{eqnarray} \tag{ 7 }$$

Note that the definition of the FIM also includes the covariance of the measurement noise, which, for the subset selection procedure, can be assumed to be identity without loss of generality.²⁰ With the Cramer-Rao bound and the FIM defined, we recall that the objective of this step is to select a subset of parameters, whose corresponding covariance is minimal. This can be achieved by minimizing the lower bound given by the FIM. To this end, a scalar metric of the FIM is optimized. A number of scalar metrics have been suggested based on different criteria.⁴⁹ Among them, the D-optimality criterion is most commonly used,²⁰ which maximizes the determinant of the FIM obtained from selected columns of the sensitivity matrix. Specifically, such an FIM is given by:

$$\begin{eqnarray}&&{{\bf{F}}}_{I}^{k}={{{\bf{S}}}_{I}^{k}}^{{\mathsf{T}}}{{\bf{S}}}_{I}^{k}={\left({{\bf{S}}}^{k}{\bf{L}}\right)}^{{\mathsf{T}}}\left({{\bf{S}}}^{k}{\bf{L}}\right),\end{eqnarray} \tag{ 8 }$$

where L is the selection matrix defined as:

$$\begin{eqnarray}&&{\bf{L}}=\left[{{\boldsymbol{e}}}_{{i}_{1}},{{\boldsymbol{e}}}_{{i}_{2}},\ldots ,{{\boldsymbol{e}}}_{{i}_{{n}_{I}}}\right].\end{eqnarray} \tag{ 9 }$$

In the above equation, ${{\boldsymbol{e}}}_{j}$ denotes the jth standard unit vector and {i₁, i₂, ..., ${i}_{{n}_{I}}$} is the set of selected parameter indices. We can now define the D-optimal selection criterion explicitly:

$$\begin{eqnarray}&&{{\varphi }_{D}^{k}}^{* }=\max {\varphi }_{D}({{\bf{F}}}_{I}^{k})=\max \,\mathrm{log}\,\det ({{\bf{F}}}_{I}^{k}),\end{eqnarray} \tag{ 10 }$$

where ${\varphi }_{D}$ is the D-optimality criterion. The D-optimal solution, ${{\varphi }_{D}^{k}}^{* }$, maximizes the determinant of the FIM, which has the effect of minimizing the volume of the confidence ellipsoid for the parameter estimates. To find ${{\varphi }_{D}^{k}}^{* }$, an integer program must be solved²⁰:

$$\begin{eqnarray}\begin{array}{cl}\mathop{\mathrm{maximize}}\limits_{{i}_{1},{i}_{2},\ldots ,{i}_{{n}_{I}}} & \,{\varphi }_{D}({{\bf{F}}}_{I}^{k})\\ \mathrm{subject}\ \mathrm{to} & \,{{\bf{F}}}_{I}^{k}={\left({{\bf{S}}}^{k}{\bf{L}}\right)}^{{\mathsf{T}}}\left({{\bf{S}}}^{k}{\bf{L}}\right)\\ & \,{\bf{L}}=\left[{{\boldsymbol{e}}}_{{i}_{1}},{{\boldsymbol{e}}}_{{i}_{2}},...,{{\boldsymbol{e}}}_{{i}_{{n}_{I}}}\right]\\ & \,{i}_{j}\in \{1,\ldots ,{n}_{\theta }\},\,j=1,\ldots ,{n}_{I}.\end{array}\end{eqnarray} \tag{ 11 }$$

This is a combinatorial problem, whose solution may be obtained by exhaustive search for small problems.⁵⁰ However, the search space grows combinatorially and this approach becomes infeasible even for moderately-sized problems. Therefore, clustering methods have been utilized to reduce the size of the problem.^20,51 Optimization routines such as genetic algorithm (GA) may also be used, albeit with no guarantee of convergence to the true optimal solution.¹⁹ Alternatively, the popular orthogonal method^14,52,53 can be used to obtain a sub-optimal solution. Particularly, Chu et al.¹⁸ showed that the orthogonal method is a forward selection algorithm that yields a sub-optimal solution for ${\varphi }_{D}^{* }$. Forward selection in this context refers to the fact that at each iteration, the algorithm only picks the best parameter to be added to the already selected list of parameters. In that sense, the orthogonal method leads to a greedy algorithm.⁵⁴ Variations of this algorithm have been suggested to consider multiple parameters at each stage, achieving more optimal results.¹⁸ The advantage of the original orthogonal method lies in its simple implementation. For example, the QR algorithm in MATLAB can be readily used for this purpose.¹⁴ Accordingly, we use the orthogonal method to find an initial sub-optimal solution that is then refined through GA.

Since only local sensitivity data are used in the above process, the resulting subset of parameters is only locally optimal and not applicable throughout the entire parameter space. A more robust approach is therefore needed. Two such approaches have been commonly used in the literature based on the R-optimality and ED-optimality criteria.²⁶ In particular, the R-optimal selection criterion is based on a worst-case scenario and yields a selection that maximizes the minimum of the D-optimality criterion over the parameter space^25,26:

$$\begin{eqnarray}&&{\varphi }_{R}^{* }=\max \,\mathop{\min }\limits_{k}\,{\varphi }_{D}({{\bf{F}}}_{I}^{k})=\max \,\mathop{\min }\limits_{k}\,(\mathrm{log}\,\det ({{\bf{F}}}_{I}^{k})),\end{eqnarray} \tag{ 12 }$$

On the other hand, the ED-optimality criterion maximizes the expected value of the D-optimality criterion over the parameter space^24,26:

$$\begin{eqnarray}&&{\varphi }_{{ED}}^{* }=\max \,{\mathbb{E}}\left[{\varphi }_{D}({{\bf{F}}}_{I})\right]=\max \,{\mathbb{E}}\left[\mathrm{log}\,\det ({{\bf{F}}}_{I})\right],\end{eqnarray} \tag{ 13 }$$

Either criterion may be used to select a robust subset of parameters that is less dependent on the nominal parameter values. We have found the R-optimality criterion to be sensitive to outlier points in the parameter space. Therefore, we use the ED-optimality criterion in this work. The corresponding optimization problem is an integer program similar to that presented in Eq. 11, which is solved using GA. In this case, the local solutions obtained during the previous step are used in the initial GA population to guide the optimization. The GA algorithm is run multiple times and the best solution is selected.

Using the procedure outlined in this section, we have selected a subset of 12 parameters for identification, which are highlighted in Fig. 1. The selected parameters cover a range of physical characteristics that impact the voltage and HFR predictions under a variety of operating conditions. For example, the selected subset includes parameters that impact the reaction kinetics, ohmic resistance, and mass transport losses. These are all parameters that may be inferred from voltage and HFR measurements alone. There are other parameters that could affect the internal conditions, but have no discernible impact on the measured signals, and therefore are not selected by the subset selection algorithm. Moreover, the selected subset may be affected by a range of modeler's decisions including the original parameter set considered (out of which a subset is to be selected), the bounds used for each parameter, and the operating conditions included in the predefined library. We also emphasize that this selection is optimal in an average sense, meaning that a better solution might exist in some regions of the parameter space. However, using the ED-optimality criterion helps ensure some degree of robustness to initial assumptions about the parameter values.

Having selected an identifiable subset of parameters, we now consider the experimental data used for identification. Importantly, we note that the data used for this purpose must be informative to ensure successful identification. This requirement motivates model-based optimal experimental design (OED), wherein the model is used to find the operating conditions that render the model predictions most sensitive to the parameters of interest. This idea is schematically shown in Fig. 2, where a numerical example is also illustrated.

Zoom In Zoom Out Reset image size

The OED procedure requires solving an optimal control problem to obtain the most informative inputs. However, this approach is not computationally tractable for a large-scale model. Accordingly, methods such as input vector parameterization have been suggested,²⁶ which searches over a discrete number of inputs to find the optimal solution. Alternatively, selecting an optimal set of conditions from a predefined library has been proposed as a computationally efficient method to solve the OED problem.³² Since we have generated sensitivity data for a library of conditions in part I,⁴⁰ we use the latter approach to design maximally informative experiments for parameter identification.

Note that when working with a library of operating conditions, the OED problem is very similar to the parameter subset selection discussed earlier and the same procedures apply. Specifically, the sensitivity matrix obtained at the k-th sample point in the parameter space using the entire library of operating conditions consists of blocks that correspond to the individual experiments:

$$\begin{eqnarray}&&{{\bf{S}}}_{I}^{k}=\left[\begin{array}{c}{({{\bf{S}}}_{I}^{k})}_{1}\\ {({{\bf{S}}}_{I}^{k})}_{2}\\ \vdots \\ {({{\bf{S}}}_{I}^{k})}_{{n}_{l}}\end{array}\right]\end{eqnarray} \tag{ 14 }$$

Here ${({{\bf{S}}}_{I}^{k})}_{i}$ denotes the sensitivity of the output predictions to the selected subset of parameters (set I) under the i-th set of operating conditions from the library of n_l conditions. We recall that in this work we have n_l = 729 conditions.

With the above block structure in mind, we note that while parameter subset selection is concerned with optimally choosing the columns of the sensitivity matrix to contain the most information, the OED focuses on choosing the best ${({{\bf{S}}}_{I}^{k})}_{i}$ blocks of the matrix for that same purpose. Therefore, D-optimality criterion is also commonly used for experimental design.^22,32,35,36 The ED-optimality criterion may likewise be used for robust design.^24,26 We note, however, that the orthogonal method is no longer applicable, since we are operating on the rows of the sensitivity matrix rather than its columns. Instead, the local D-optimal solution can be obtained through convex relaxation of the underlying integer program³²:

$$\begin{eqnarray}\begin{array}{cl}\mathop{\mathrm{maximize}}\limits_{{\eta }_{1},{\eta }_{2},\ldots ,{\eta }_{{n}_{l}}} & \,\mathrm{log}\ \det ({\tilde{{\bf{F}}}}_{I}^{k})\\ \mathrm{subject}\ \mathrm{to} & \,{\tilde{{\bf{F}}}}_{I}^{k}=\sum _{j=1}^{{n}_{l}}{\eta }_{j}{\left({{\bf{S}}}_{I}^{k}\right)}_{j}^{{\mathsf{T}}}{\left({{\bf{S}}}_{I}^{k}\right)}_{j}\\ & \,0\leqslant {\eta }_{j}\leqslant \displaystyle \frac{1}{{n}_{u}},\,j=1,\ldots ,{n}_{l}\\ & \,\sum _{j=1}^{{n}_{l}}{\eta }_{j}=1,\end{array}\end{eqnarray} \tag{ 15 }$$

where n_u is the experimental budget, i.e., the number of experiments to be selected from the library. Determining the experimental budget requires consideration of the actual experimental cost, as well as the computational cost required to evaluate the model for the selected experiments. Additionally, the number of selected experiments must provide enough information for parameter identification. Taking these into consideration, here we set the number of experiments to n_u = 6.

The solution to the above convex program is a set of n_u experiments that are D-optimal at sample point k in the parameter space. In other words, this solution is only locally optimal. Therefore, to find an experimental design that is less dependent on the initial parameter values, the ED-optimality criterion of Eq. 13 is used.

Briefly, the robust OED procedure is as follows:

• 1.  
  Choose the number of experiments, n_u, to be selected from the library.
• 2.  
  Solve the convex program 15 to obtain the D-optimal design at the sampled points in the parameter space.
• 3.  
  Use GA to obtain the ED-optimal design considering all of the sampled points. The D-optimal designs from the previous step are used in the initial population to guide the GA.

The experimental design resulting from this procedure is shown in Table I. The identification results based on this robust OED procedure are compared with another experimental design using Latin Hypercube Sampling (LHS). Importantly, we use the same number of experiments (n_u = 6) for the LHS design. The resulting LHS design is given in Table II and the identified parameters using both experimental design procedures are compared below.

So far, we have selected a subset of parameters for identification and designed informative experiments for that purpose. The last critical piece of the model parameterization framework is the identification algorithm. A good identification algorithm should be sufficiently constrained. An algorithm with too many degrees of freedom can lead to over-fitting and poor parameter identification.⁴⁴ Regularization methods such as ridge regression^45,46 and LASSO⁵⁵ are commonly used to constrain the identification process and formulate a well-conditioned problem. However, in the context of parameter identification, application of these methods typically requires prior information about the parameter values. Since we assume no prior information about the parameters, here we focus on using the structure of the data to constrain the identification procedure. We also recall that the parameter subset selection approach, described above, indeed regularizes the problem to a certain extent and the identification algorithm of this section is developed for further regularization.

The underlying principle of the algorithm proposed here is that the cell response shows varying sensitivity to parameter perturbations depending on the load as seen in Fig. 3. For instance, a parameter that affects reactant transport resistance might only become identifiable at higher loads. Therefore, we use this particular sensitivity structure to identify the selected parameters in multiple steps. Particularly, we divide the polarization curve and the corresponding HFR data into 3 sections based on the load (see Fig. 3): (1) low current density or kinetic region, (2) medium current density or ohmic region, and (3) high current density or mass transport region. Alternatively, the regions may be defined based on the sensitivity data, for example, by classifying the regions based on average voltage sensitivity (low, medium, and high) under different conditions. However, here we utilize the existing knowledge in the fuel cell literature, where the three aforementioned regions are often well-defined. Encoding this knowledge into the framework allows for simpler formulation.

Zoom In Zoom Out Reset image size

Following the above division, the selected parameters are grouped based on the corresponding sensitivities of the predicted outputs in each region. Specifically, the first parameter group consists of those parameters that are identifiable in the low current density region of the polarization curve, while the second and third groups include the parameters that become identifiable at medium and high current densities, respectively. Identifiability in this context is measured by the average output sensitivity to the parameter of interest; a parameter is deemed identifiable in a particular region if the corresponding sensitivity is greater than a threshold value S_th, which is chosen to be 0.1 in this work. This grouping procedure based on output sensitivities is shown in Fig. 4.

Zoom In Zoom Out Reset image size

Following the above grouping of parameters, the proposed multi-step identification algorithm is as follows:

• 1.  
  Identify the parameters in the first group using data obtained at low current densities. The parameters in the other two groups are fixed at their nominal values.
• 2.  
  Identify the parameters in the second group using data obtained at medium current densities. Parameters in the third group are kept at their nominal values, while those that were identified in the first step are further refined by allowing them to vary within a contracted search space.
• 3.  
  Identify the parameters in the third group using data obtained at high current densities. Previously identified parameters are allowed to vary, but their respective search spaces are contracted. Particularly, the search spaces for the parameters in the first group are shrunken further compared to their ranges in the second step of the algorithm.
• 4.  
  Refine all identified parameters using the entire data set. The search space for all parameters are further contracted in this step to constrain the problem. This step is only used to refine the identified values.

The general idea is that the output sensitivities to parameter perturbations typically increase as the load is increased. Therefore, it is helpful to start the identification process by only using the data at low current densities and as our confidence in the parameter estimates increases, move to higher loads and identify the remaining parameters. The increased sensitivity at higher loads also means that it is critical to allow previously identified parameters to be refined as the load is increased. Such cumulative fitting enables more accurate parameter estimation.^31,32

Alternative implementations

As described above, the proposed model parameterization framework has three main components, namely, parameter subset selection, optimal experimental design, and the multi-step identification procedure. In this section, we demonstrate the utility of each of these components. To this end, we use the model to generate synthetic experimental data with a set of nominal parameter values (${{\boldsymbol{\theta }}}^{* }$). The synthetic data are then used to identify the model parameters starting from perturbed values. Particularly, the optimization cost function is the weighted norm of the error:

$$\begin{eqnarray}&&J({\boldsymbol{\theta }})={w}_{V}\displaystyle \frac{\parallel {{\boldsymbol{e}}}_{V}({\boldsymbol{\theta }}){\parallel }_{2}}{\sqrt{{n}_{V}}}+{w}_{R}\displaystyle \frac{\parallel {{\boldsymbol{e}}}_{R}({\boldsymbol{\theta }}){\parallel }_{2}}{\sqrt{{n}_{R}}},\end{eqnarray} \tag{ 16 }$$

where ${{\boldsymbol{e}}}_{V}$ and ${{\boldsymbol{e}}}_{R}$ denote the prediction error vectors for voltage (measured in mV) and HFR (measured in mΩ · cm²), respectively, ${n}_{V}$ and ${n}_{R}$ are the number of elements in ${{\boldsymbol{e}}}_{V}$ and ${{\boldsymbol{e}}}_{R}$, and ${w}_{V}$ and ${w}_{R}$ are the corresponding weights. Since the HFR varies in a smaller range compared to voltage, here we use ${w}_{V}$ = 1 and ${w}_{R}$ = 10. Given that we apply large perturbations to the nominal values and the above cost has many local minima, a gradient-based optimizer would fail to find the optimal solution. Use of a global optimizer is thus warranted. We have experimented with both GA and particle swarm (PS) algorithms and found that with the same level of effort to tune the optimizer hyper-parameters, the latter can find a more optimal solution. Therefore, we use PS to identify the selected parameters.

To highlight the importance of each part of the framework, we consider 4 different cases for identification:

• 1.  
  Only the selected parameters (see Fig. 1) are perturbed from their nominal values. The OED conditions (Table I) are used along with the multi-step identification procedure to identify the parameters.
• 2.  
  Only the selected parameters (see Fig. 1) are perturbed from their nominal values. The LHS conditions (Table II) are used along with the multi-step identification procedure to identify the parameters.
• 3.  
  Only the selected parameters (see Fig. 1) are perturbed from their nominal values. The OED conditions (Table I) are used along with a single-step identification procedure, where all of the selected parameters are identified using the entire data set.
• 4.  
  All parameters, except those that were removed during the collinearity analysis, are perturbed from their nominal values. The OED conditions (Table I) are used along with the systematic identification procedure to identify the parameters.

In all of the above cases, the applied parameter perturbations are random. In particular, the perturbations to the selected parameters are drawn from ${ \mathcal U }(-0.3,0.3)$ (i.e., uniform distribution between −0.3 and 0.3), while the perturbations to the parameters that are not included in the identification (case 4) are drawn from ${ \mathcal U }(-0.15,0.15)$. We also note that these perturbations are applied to the scaled parameter values and would amount to significant deviations from the nominal values on the original parameter axes. Additionally, it should be pointed out that the nominal values do not coincide with any of the sample points in the parameter space that are used to generate the sensitivity data, select the parameter subsets, and design the optimal experiments. These measures are taken to properly examine the robustness of the proposed framework.

The identified parameter values for the first two cases are provided in Table III along with the nominal and perturbed values. Moreover, the identification results for all of the above scenarios are shown in Fig. 5 in terms of the relative error of the identified parameters. Particularly, Fig. 5a compares the results from the first two scenarios and highlights the utility of OED. As seen in the figure, the maximum relative error in the parameter estimates is only 8.4% with the OED (for n_e), whereas the LHS design results in a maximum relative error of about 170% (for ξ_diff,mb). Moreover, the mean relative errors are 3.9% and 28.7% for the optimal and LHS experimental designs, respectively. In addition to these relative errors, the model voltage and HFR predictions are shown in Fig. 6. Particularly, the predictions with the nominal and identified parameter values are shown for both OED and LHS cases. Two main observations follow: first we note that the OED conditions lead to more pronounced variations in the predictions; therefore, they are more informative. Second, we observe that for both sets of conditions, the identified parameter values reproduce the synthetic experimental data very well. This result emphasizes that a good fit to experimental data does not readily imply an accurate parameter identification. These observations, combined with the relative parameter estimation errors mentioned above, demonstrate the importance of having maximally informative experimental data for parameter identification, which can be obtained through model-based OED.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To show the effectiveness of the proposed multi-step identification algorithm, Fig. 5b compares the relative errors obtained using the proposed algorithm with those obtained from a single-step algorithm. The single-step algorithm uses the entire data set to identify all of the selected parameters. Note that in both cases we use the data based on OED conditions of Table I. The figure shows that the maximum and mean relative errors for the single-step case are 33.8% and 10.5%, respectively. This showcases the improved identification capability that is achieved by constraining the problem using the data structure.

Lastly, the subset selection approach is verified by comparing the identification results of case 1 with those of case 4, wherein the parameters that are not included in the selected subset are also perturbed from their nominal values. These results are shown in Fig. 5c. We note that if the subset selection algorithm is successful, then the parameters that are not included in the subset are not expected to have a significant impact on the output predictions. Therefore, perturbations to these parameters should not lead to considerable bias in the identified parameters. We further note that the parameters that were left out of the selected subset due to their collinearity with another parameter are not perturbed in case 4, as their perturbations can result in bias in the estimates of their collinear counterparts. As can be seen in the figure, perturbing the parameters outside of the selected subset does not have a significant impact on the relative errors in the parameter estimates. In particular, the maximum and mean relative errors are 11.6% and 3.8%, respectively, which are on par with those of case 1. Therefore, the subset selection algorithm has been successful.

These results highlight the importance of three key elements in successful parameter identification: selecting a suitable subset of parameters for identification, optimally designing experiments to generate maximally informative data for identification, and constructing a sufficiently regularized identification algorithm. Moreover, by extending the analysis beyond the vicinity of a point in the parameter space, these procedures can be made robust to initial assumptions about the nominal parameter values. We also emphasize that while a particular fuel cell model is used as an example to demonstrate the various steps involved in the parameter identification framework, the proposed procedures are generally applicable to different types of models, including other fuel cell models.

Since effective model parameterization is crucial to successful predictive modeling, these results can have profound implications for model-based design and control strategies. In addition, parameter identification plays a significant role in inferring physical characteristics of a fuel cell from experimental data. For example, parameter identification with cell-level measurements such as voltage has been used to diagnose fuel cell degradation issues,¹⁰ determine catalyst activity,^8,56 and evaluate mass transport losses.^9,30 Numerous electrochemical impedance spectroscopy (EIS) studies provide further examples of such efforts.⁵⁷ Against this backdrop, the methods presented in this work also gain critical importance in further helping improve the current understanding of various aspects of fuel cell performance and durability by enabling more accurate parameter identification.

A framework for systematic parameter identification for PEM fuel cell models is developed. The framework consists of three main components: (1) selecting a robust subset of parameters for identification, (2) using model-based design of experiments to obtain a robust set of maximally informative data for identification of the selected parameters, and (3) using the structure of the data to formulate a regularized identification problem. This framework is based on the extended local sensitivity analysis that is explained in the first part of the study. Particularly, sensitivities of model predictions to parameter perturbations are used to determine the parameters that have the largest unique impact on the predicted outputs. This is formalized through the D-optimality criterion. Importantly, rather than maximizing the local D-optimality criterion, we maximize its expected value over a number of sample points in the parameter space to ensure the resulting subset is less dependent on the nominal parameter values. The same procedure is applied to select the most informative operating conditions from a predefined library of conditions. Lastly, a multi-step identification procedure is proposed to regularize the optimization problem. The effectiveness of each of the components in the framework is verified through a case study, wherein synthetic experimental data generated with nominal parameter values are used to identify the selected parameters. The results demonstrate the efficacy of the proposed subset selection, optimal experimental design, and the multi-step identification procedures in improving parameter identification results. While the methods of this work are applied to a full cell model, they can also be utilized for sub-models that focus on a specific phenomenon. Therefore, we believe that adoption of such methods would help ensure that the models are properly parameterized and their full predictive capabilities are realized.

Financial and technical support by Ford Motor Company is gratefully acknowledged.

The complete list of model parameters considered in this work is provided below. Particularly, the parameter descriptions and their units and corresponding ranges are given. Importantly, note that all parameters considered here are bounded to specific ranges, which are determined based on the corresponding uncertainty in the parameter values reported in the literature.