Data-driven parameterization of polymer electrolyte membrane fuel cell models via simultaneous local linear structured state space identification

Data-driven parameterization of a 0D dynamic PEMFC stack model. Computational efficient parameterization methodology. Simultaneous identification of analytic local linear models. Parametric output sensitivitybased identifiability analysis. Experimental validation on a commercial 30 kW PEMFC stack. a r t i c l e i n f o


Introduction
Polymer Electrolyte Membrane Fuel Cells (PEMFCs) are regarded as a promising candidate towards replacing the internal combustion engine and enabling an environmentally friendly, sustainable mobility. However, to achieve a wide-spread market penetration, significant improvements, especially regarding the specific costs and durability [1], have to be made in order to be competitive with the internal combustion engine. To extend the lifetime, degrading operating conditions (e.g. fuel starvation, excessive pressure differences across the membrane, flooding or dry-out) have to be prevented [2,3], which requires precise control of the balance-of-plants components, especially during transient operation. Therefore, advanced control methods are required. A prerequisite for this, however, is the availability of a real-time computable, yet sufficiently accurate, fuel cell model. To strike a balance between complexity and computational efficiency, control oriented fuel cell models are typically of low spatial order (0D/ 1D), with the governing dynamic equations derived from the transient mass balances of species, and an algebraic electrochemical model describing the relation between the thermodynamic states of the fuel cell, current and voltage [4e11]. Control oriented models have to describe the specific fuel cell system with sufficient accuracy, which is rarely achieved by adhering to a white-box modeling approach (i.e. all the model parameters can be derived from fundamental constants and geometric considerations). For this reason, the model parameters are usually determined with the help of experimental data by minimizing the discrepancy between simulated and measured output signals [12,13].
The experimental parameterization of a PEM fuel cell model poses a considerable challenge: 1. The resulting model is typically highly non-linear and numerically stiff. 2. A large number of parameters (with varying levels of significance) are to be determined. 3. The associated parameter optimization problem is typically non-linear, non-convex with multiple local minima, and (possibly) ill-posed.
Owing to these difficulties, heuristic optimization techniques (such as differential evolution, genetic algorithms, particle swarm optimization) are frequently used [14e16]. Although satisfactory results in terms of the fitted model output can thereby be obtained, the heuristic optimization of non-linear, dynamic models leads to a substantial computational effort. Additionally, the optimization problem being illposed is often a symptom of an over-parameterized model [17,18] resulting in a non-unique solution. When applying heuristic optimization techniques, these symptoms of an over-parameterized model are neglected, giving a false sense of security, that a high accuracy in terms of the model output must be the result of having obtained the correct model parameters. The present work aims to accelerate the task of model parameterization and, additionally, to analyze the significance (and therefore identifiability) of model parameters with respect to the available experimental data. For that, the following estimation procedure is proposed: The nonlinear fuel cell model is analytically linearized such that for any given steady state input a local linear model is obtained. As linear dynamic systems can be solved very efficiently, the computation time is drastically reduced as compared to the original non-linear system. Since the linearization is done analytically, the resulting linear dynamic model is still depending on the same physical parameters as the non-linear model. Written in state space form, the estimation of such parametric local linear models termed structured state-space identification [19e21].
Since the PEM fuel cell model is highly non-linear, it is unreasonable to assume that all parameters describing the non-linear model can be consistently estimated based on a single local linear model. It is therefore proposed to identify multiple linear structured state space models simultaneously. Note, that this method of identifying a non-linear system via multiple local linear approximations is closely related to the multimodel approach [22] or local linear model networks [23]. However, within the classical multi model framework, one is typically only concerned with describing the input/output behavior of a non-linear system with sufficient accuracy, e.g. black-box modeling. Therefore the physical meaning of states and estimated parameters in the resulting model are normally not accessible. Additionally, each local model depends on its own parameter vector, and the local model structure needs to be determined from the data as-well, whereas for the analytic linearization approach, the local model structure is determined by the non-linear model, and each model shares the same physical parameter vector and state space.
To facilitate the estimation of local linear operating point models, the experiments need to be designed accordingly: At each operating point, the fuel cell system is dynamically excited using the available input channels with small enough i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x amplitudes, as to justify the assumption of linearity around an operating point. So instead of having a single experiment covering the complete operating space (suitable for the direct estimation of a non-linear model), several local experiments with reduced excitation amplitudes around multiple operating points are carried out instead. A schematic representation thereof can be seen in Fig. 1. To facilitate an a-priori analysis of the significance and identifiability of the model parameters with respect to the available data (and local linear models), a semi-analytic algorithm for the computation of parametric output sensitivity, and subsequently the Fisher Information Matrix (FIM), is developed. With the knowledge of the FIM, the identifiability of parameters can be concisely analyzed [17,24].
In summary, the proposed identification method reduces the needed computational effort to identify a PEMFC model by multiple orders of magnitude by identifying multiple local linear structured state-space models instead of a single nonlinear one. The analytically linearizable model allows an additional a priori semi-analytic parameter identifiability analysis.
The remaining paper is structured as follows: In section Fuel cell stack model, the control-oriented fuel cell stack model is presented. In section Simultaneous local linear structured state space estimation, the parameterization methodology utilizing the simultaneous estimation of local linear structured state space models is discussed in detail. Section parametric sensitivity and fisher information regards the identifiability analysis of parameters. Section Experi mental setup gives a brief overview of the experimental setup that is used to conduct the experiments on the 30 kW PEM fuel cell stack. The obtained results are discussed in section Results.

Fuel cell stack model
In this section, the fuel cell stack model used for the experimental parameterization is described. The described PEMFC model is based on the control-oriented zero-dimensional model presented in Ref. [4]. The switching between condensation and evaporation of water is adapted due to the necessity of analytical differentiability. Additionally, the electrochemical model is exchanged, and the one proposed in Ref. [25] is used. The reason is that the latter gives better insights into the internal states. Further details in this regard are discussed in the following subsections. The overall structure of the model considered in this work can be seen in Fig. 2. Both the cathode and the anode are divided into three consecutive zero dimensional volumes, interconnected with linear nozzle equations. The first volume of the anode and cathode respectively is denoted as the supply manifold (sm) and the last volume as exit manifold (em). The center manifolds describe the respective volumes of cathode and anode (ca/an), whose species concentrations are used to determine the electrochemical reaction. The model dynamics are described by the transient mass balance of the interconnected zero dimensional volumes. Water vapor and nitrogen are exchanged between the center manifold of cathode and anode. The membrane humidity is described by the average humidity of cathode and anode, filtered by a first order system with an associated time constant as to approximate the transient drying and wetting. The internal dynamic states and their derived quantities are then used in the electrochemical model to describe the relation between current and voltage. The temperature of the fuel cell stack is assumed to be uniform and externally controlled. Gas recirculation in the anode and purging is considered. It is to be noted that the main focus in this work lies in the experimental parameterization. The presented model is to be understood as a candidate model in explaining the experimental data and not as a proposition and validation of modeling assumptions.

Cathode
The total mass flow of humid air flowing into the cathode-side supply manifold is characterized by its total mass flow rate W ca in and relative humidity 4 ca in . Assuming a constant mass fraction  of nitrogen in standard dry air and using basic thermodynamic calculations (see for example [4]), the mass fractions of species y ca N2; sm , y ca O2; sm and y ca vap; sm and the specific gas constant of the gas mixture inside the supply manifold R ca sm can be calculated. The transient change of pressure inside the supply manifold is then given by the ideal gas law and mass balance as Thereby, T denotes the uniform stack temperature, V ca sm the volume of the supply manifold, and W ca sm the total mass flow between the supply manifold and the cathode, which is described using the linearized nozzle flow equation The oxygen consumption and the water vapor generation due to the electrochemical reaction is given by and W vap; gen ¼ n cells I M H2O 2F ; respectively, where n cells denotes the number of cells in the stack, I the current, and F the Faraday constant. The nitrogen crossover through the membrane is considered via The net water vapor flow through the membrane W vap, mem is defined later when discussing the membrane sub-model. For total mass flow between cathode volume and exit manifold again, a linear nozzle equation is utilized The overall pressure inside the cathode is given by the sum of partial pressures: with where R denotes the universal gas constant, M k the molar mass of species k and V ca the cathode volume. Different approaches have been proposed to deal with the phase change and generation of liquid water. In Ref. [4], it is assumed that any water mass exceeding the saturated vapor mass is instantaneously condensed into liquid water (and vice-versa). The partial pressure of water vapor is then fixed to the saturation pressure. It is evident, that this approach introduces numerical difficulties and discontinuities to the model due to the infinitely fast evaporation and condensation. Additionally, when transitioning from unsaturated to the saturated region, the state-space of the model is effectively reduced due to the algebraic constraint of constant water vapor pressure. In contrast, transient phase change relations have been considered in Refs.
[26e28], which are preferable as thereby liquid water and water vapor can be treated as individual states in the model.
The phase change rate is then expressed by where k cond and k evap are the condensation and evaporation rate constants, and J the transition function between condensation and evaporation. Typically, a discontinuous hard switching is implemented as a transition function However, due to the intended analytic linearization of the model, continuous differentiability is required. Therefore the hard switching transition function is approximated using the following sigmoid function: Thereby, the constant k sig can be used to define the width of the transition area and therefore the trade-off between approximation error and numerical stiffness.
The exit manifold pressure is given similar to the supply manifold: The concentration of species inside the exit manifold is thereby assumed to be the same as in the cathode volume.
where p atm denotes the atmospheric pressure. In the case of an actuated back-pressure valve, fa bp 2Rr0 a 1g is the time dependent valve opening position, which can be used as an additional degree of freedom in controlling the cathode pressure.

Anode
The modeling of the anode side is analogous to that of the cathode when considering open end operation. However, when the fuel cell system to be identified is subjected to a dead-end operation, or includes an anode recirculation loop, as depicted in Fig. 2, the following factors have to be taken into account: When disregarding the effects of leakages in the piping, the anode pressure dynamics exhibit an integrating behavior. If now the simulation is driven by the measured hydrogen mass-flow, any measurement noise, bias, or un-modeled disturbance lead to a random walk or continuous drift in pressure due to the integrating behavior. As such, a pressure driven anode-side simulation is to be preferred. The measured hydrogen mass-flow is then only used as to determine the mass fractions of species inside the anode supply manifold: and y an vap; sm ¼ The total mass flow between anode supply manifold and center manifold is given as Mass balance of the individual species inside the anode center manifold leads to _ m an H2 ¼ y an H2; sm W an sm À W H2; react À y an H2 W an ; _ m an N2 ¼ y an N2; sm W an sm þ W N2; cross À y an N2 W an ; _ m an vap ¼ y an vap;sm W an sm À W vap;mem À y an vap W an ; _ m an liq ¼ W an phase À W an liq : The overall pressure and partial pressures of species in the anode are calculated analogously as on the cathode side ( (12) and (13)).
Mass flow between anode center manifold and exit manifold is given by The anode exit manifold pressure is described by The mass flow due to purging W purge is described similar to the back pressure valve (18), however with a respective nozzle coefficient k purge and purge valve actuation signal a purge . The mass flow due to internal leakage in the anode (28) is, on the one hand, physically motivated due to the fact that the anode piping and membrane are not perfectly impermeable and, on the other hand, due to the aforementioned numerical reasons of the anode exhibiting integrating behavior if modeled as being a hundred percent gas tight.

Membrane and water exchange
Water vapor is exchanged between cathode and anode Thereby, A cell denotes the cell area and N v the total molar flux of water vapor through the membrane. With d m as the membrane thickness, the total molar flux can be described by the superposition of electro-osmotic drag and back diffusion where the electro-osmotic drag coefficient n d (l m ) and the diffusion coefficient of water through the membrane D w (l m ) are typically understood to be functions of the membrane water content l m . Numerical expressions and approximations have been obtained by fitting experimental data, for example by Dutta et al. [29]. Thereby the water diffusivity is described by a set of piecewise linear functions. However, this expression would violate the requirement of continuous differentiability for the subsequent parameterization approach. The hard switching could again be numerically smoothed as to achieve higher order differentiability. However, as pointed out and summarized in Ref. [30], numerous alternative experimental results and fitting functions for the electro-osmotic drag and diffusion coefficient have been obtained and proposed in literature, with drastic differences in structure and numerical values. As a result, adhering to a parsimonious modeling principle, simple linear relations for the electro-osmotic drag and diffusivity coefficient are therefore used [31,32] in this work. The membrane water content l m with respect to the water activity a m is described by the following commonly found relation [33] l ¼ 0:043 þ 17:81a m À 39:85a 2 m þ 36a 3 m : The membrane water activity a m is typically obtained by averaging the water activity in cathode, and anode [4,5,27]. However, to account for the dynamic wetting and drying of the membrane, a simple first order dynamic with a characteristic time constant t m is introduced with The water vapor concentrations in (30) are described by where r mem, dry denotes the membrane's dry density and M mem, dry the membrane dry equivalent weight.
The influence of the membrane water content l and temperature on the ionic conductivity is considered by the following empirical equation [4]: The ohmic resistance is then given by where d m denotes the membrane thickness, A cell the effective cell area, and R 0 the contact ohmic resistance.

Electrochemical model
The electrochemical model links the internal thermodynamical states of the 0D-FC model with the current and voltage. A highly accurate electrochemical model is sought, especially, as the measured voltage and current can give invaluable insights into the (otherwise unmeasured) internal states. For control oriented fuel cell models, a variety of electrochemical models have been proposed in literature. In Ref. [27], a simple electrochemical model, equating the Nernst potential, activation, concentration, and ohmic losses, is present, however, the model exhibits little dependency on the thermodynamic states of the fuel cell as lumped parameters of the limiting and exchange current density are fitted. In order to account for this fact, simple electrochemical models are often extended using a black-box/data-driven approach [4,34], where the parameters are extracted for several operating conditions, and the dependency of the electrochemical parameters to the thermodynamic states of the cell are afterwards fitted and interpolated. Highly accurate results have thereby been obtained, however as an obvious drawback, a vast amount of experimental data is thereby required. In this work, the electrochemical model developed in Ref. [25] is applied, as it exhibits suitable extrapolation capabilities during different operating conditions on a reasonable amount of fitting parameters. However, the effect of liquid water accumulation in the gas diffusion layer, reducing the effective diffusivity [35], is not considered. The model is derived based on the following major assumptions [25]: 1. The problem is isothermal, however the electrodes could have different temperatures.

Charge transfer coefficient in electrochemical reactions
is assumed to take the value of 0.5. 3. Concentration of electrons and protons is approximately uniform in the electrodes. 4. Catalyst layers are infinitely thin, and a uniform potential is assumed for them. 5. The gasses are treated as ideal. 6. The diffusion constants are pressure insensitive. 7. The diffusion system is always bicomponential and therefore described by scalar diffusion constant. 8. There is no liquid water in the GDLs and channels. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x 9. There is no convective transport in the GDLs. 10. The membrane assumes no gas crossover. 11. The membrane is ideally and constantly hydrated at all times.
Following, the voltage per cell V cell is calculated with In (37) k B denotes the Boltzmann constant, e 0 the elementary charge, Z the number of electrons transferred in the electrochemical reaction, Dg 0 the specific Gibbs free energy, T 0 the reference temperature, Ds the specific entropy, and E 0 the activation energy. ca and an indicates the corresponding electrode. C describes the unitless concentration, which can be obtained by dividing the corresponding partial pressure by the atmospheric pressure. I L,ca denotes the cathode limit current with the relation where F is the Faraday constant, and CD ca is the combined diffusion coefficient for the cathode. Furthermore, I ca 0 describes the cathode intrinsic exchange current with the relation where K ca is the cathode intrinsic exchange current parameter. The proceeding for the anode side is analogous. Finally, the stack voltage V can be evaluated with

Overview of model properties
Equations (1)e(34), together with the electrochemical model, can be written as a non-linear state space system of the form y ¼ gðx; u; qÞ: The model exhibits a direct feed-through, i.e. from current to voltage. The input vector u, state vector x and output vector y are depending on the fuel cell stack system setup and the available measurements. For the subsequent parameter estimation based on the experimental data, they are given as Since anode recirculation is considered, the pressure in the supply manifold p an sm is used as an input. The mass flow from the hydrogen tank into the anode W an in is only used to determine the concentrations inside the supply manifold, where the pure, dry hydrogen from the tank (e.g. y vap, tank ¼ 0) is mixed with the recirculated mass flow.
The state derivative function (41) and output function (42) are continuously differentiable as required by the parameterization method.

Simultaneous local linear structured state space estimation
In this section, the estimation procedure for the fuel cell model is described. At first, the commonly applied minimization of the simulation error of the non-linear model is discussed. Afterwards the approach of analytic local linearization is introduced, leading to the simultaneous estimation of local linear structured state space models.

Parameter estimation via output error minimization
In order to estimate the unknown parameter vector q of the model described by equations (41) and (42), a non-linear optimization problem is formulated based on the error minimization between simulated and measured outputs. A suitable scalar objective function thereof constitutes Thereby, Tr (,) denotes the trace of a matrix, i.e. the sum of diagonal elements. The matrix e Y2R NÂny contains the measured output signals, with N sampled time instances of the n y channels. The matrix Y(q) denotes the corresponding simulated model outputs (subject to the same input signals as in the experiment). The matrix Q y is used to weigh the importance of different outputs with respect to each other. The  (44) is used to penalize deviations from the initial parameter vector b q 0 . On the one hand, the choice of Q q can be motivated by a practical point of view, in the sense, that it reflects the knowledge and reliability of a-priori assumed parameters (expert's knowledge). On the other hand, by introducing the second term in (44), a regularization of the subsequent optimization problem is achieved, giving additional degrees of freedom as to influence the bias-variance trade-off [36,37].
The optimal parameter estimate is formally stated as However, an analytic solution is, in general, not obtainable. Therefore, iterative optimization algorithms are applied Thereby, 4 opt (,) denotes an arbitrary optimization algorithm, and I the number of iterations after which the algorithm is stopped. At each iteration, in order to evaluate Jð c q i Þ, the non-linear dynamic model (42)e(43) needs to be simulated which can formally be stated as YðqÞ ¼ 4 sim ðx 0 ðqÞ;Ũ; qÞ: Note that the initial state vector x 0 (q) is parameter dependent and requires a solution of the implicit equation fðx 0 ; e u 0 ; qÞ ¼ 0: Ideally, the final parameter vector b q I serves as a suitable approximation of the optimal solution b q * . However, it is obvious that the obtained result is heavily influenced by the choice of the optimization algorithm 4 opt , and applied solver for the dynamic model 4 sim . In practice, the following pitfalls can prohibit the convergence to a suitable solution: 1. The initial parameter vector q 0 is improper. 2. The current parameter vector during optimization q i leads to numerically stiff model and the solver 4 sim (,) does not terminate or fails completely. Likewise, the solver to obtain the steady state solution x 0 (q i ) from (48) can fail. 3. Convergence to local minimum of the objective function (45) or an unsuitable optimization algorithm 4 opt (,) with insufficient settings is used.
It is evident that a certain iterativeness due to trial-anderror cannot be avoided. Therefore it is all the more important to accelerate every single parameterization run as much as possible. Especially with regard to the dynamic simulation, the switch to local linear models enables the use of very efficient linear system solvers and, therefore, a drastic increase in computational speed. Additional, improper candidate parameter vectors during optimization, leading to undesired model properties such as instability, can be easily detected and penalized.

Structured state space model
A local linear operating point model is given by where the system matrices are defined as The inputs and states of each local model are considered as deviation around its operating point Dx ¼ x À x 0 ðqÞ; whereas the local model output (50) is given in its absolute values, for which the steady state output needs to be computed via By calculating the partial derivatives in (51)e(54) analytically, as opposed to a numerical approximation, the functional dependency of the linear system matrices with respect to the parameter vector q is retained. To carry out the partial differentiation of the non-linear state space model programmatically, several software packages are readily available. In this work, the MATLAB symbolic toolbox [38] was used. These analytic calculations are done once in a preprocessing step of the non-linear model, and the resulting matrix functions are stored. For the numeric simulation during optimization, a specific local linear model is efficiently generated by inserting the numerical values of the operating point and parameter vector into these preprocessed state space matrix functions.

Simultaneous estimation of multiple structured state space models
As the fuel cell model (41)e(42) is highly non-linear, it is unlikely that a single local linear model, together with a single local experiment, is sufficient for estimating the complete parameter vector q. It is therefore proposed to estimate multiple local linear structured state-space models simultaneously. An appropriate objective function can thereby be given as l m J m ðqÞ þ ðq À q 0 Þ T Q q ðq À q 0 Þ; Y m the corresponding measured outputs of the local experiment. Note that by increasing the number of local linear models, the number of parameters to be identified does not increase. To carry out the minimization of (58), it is advisable to have multiple optimization algorithms at hand. In the present work, a hybrid optimization approach is used as they have been found to outperform purely global optimization techniques for the parameterization of large transient models [39].

Parametric sensitivity and Fisher information
An important measure for determining the quality and significance of estimated parameters is its associated parameter covariance Thereby, q U denotes the "true" parameters of the system and b q its estimates. Assuming a consistent estimator, unlimited experimental data, and a correct model structure, it follows that b q /q U and covð b q Þ/0. However, in the practical case of limited experimental data, the knowledge of the parameter covariance allows for a statistical assessment of confidence of the obtained parameters. A direct numerical approximation of (60) is in general not feasible, however, there exists a lower limit defined by the Cram er-Rao inequality [40]: where F denotes the Fisher information matrix, which, for a dynamical system, can be given as [41]: with S e as the error covariance matrix and The matrix 4(k, q) is commonly referred to as the (local) parametric output sensitivity. Note that (62) and (63) holds for any model structure, e.g. for both, the non-linear and the local linear structured state space model. An obvious choice for computing (63) would be a finite difference approximation of the derivative, but for the present fuel cell model, this has been found to be numerically sensitive due to the large spread of parameter values, the stiffness of the model, and the trade-off between the numerical approximation and round-off error when computing the difference approximation. However, for the identification via multiple local models, the parametric state space structure can be effectively exploited as to robustly and efficiently compute the parametric output sensitivity.
With (62) the confidence interval of the estimated param- where z a/2 is the z critical value at significance level a/2 and ðF À1 Þ ii is the i-th diagonal element of the inverse Fisher information matrix.
Analytic computation of the parametric sensitivity for structured state space models As a first step, in order to account for the discrete nature of the experimental data, the structured state space model (49)e(50) is discretized yðkÞ ¼ CðqÞDxðkÞ þ DðqÞDuðkÞ þ y 0 ðqÞ; (66) with the discrete system matrix and input matrix Under the assumption of constant inputs during the sampling interval Dt, (65) is the exact solution of (49). The parameter sensitivity (63) for the discrete structured state space model can then be written as Thereby, x i (k) ¼ dDx(k)/dq i has been introduced as the parametric state sensitivity, which can be recursively computed by In order to evaluate (69) and (70), the matrix derivatives therin need to be computed. While dC/dq i and dD/dq i can be readily carried out analytically in a preprocessing step of the model, the matrix derivatives vA k /vq i and vB k /vq i in (70) must be considered separately, due to the intermediate step of discretization. However, an elegant solution therefore has been given in Ref. [42], where it has been shown that the discrete system matrices (67) and (68), as well as its derivatives with respect to the parameters, can be efficiently calculated by i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x 6 6 6 4 Thereby, the evaluation of the discrete matrix derivatives is shifted to the computation of the continuous matrix derivatives, which can readily be done in a pre-processing step of the model.

Fisher information based identifiability analysis
In order to determine the lower bound of the parameter covariance (61), the Fisher information matrix needs to be inverted. To analyze the invertibility and well-posedness of the inversion, it is advised to conduct a singular value decomposition of the Fisher information matrix [43] where U and V are unitary matrices and with s 1 ! / ! s np as the singular values. In the case of noninvertibility, the Fisher information matrix is rank deficient, i.e. some singular values are identical to zero. The singular value decomposition can therefore be factored into This effectively means that the Fisher information can be, without any loss, described by a parameter space of reduced dimensionality. Consequently, V 1 represents the identifiable subspace of the parameters and V 2 the unidentifiable parameter subspace. To investigate which parameters are causing the problems, the column vectors V 2 (e.g. base vectors of the unidentifiable parameter subspace) need to be analyzed. In practice, a singular Fisher information matrix for a first principle model can arise when the effect of a parameter algebraically cancels out, or when the output is only effected by a combination of parameters as in the trivial case: Regardless of how the input signal u(t) is chosen, it is not possible to uniquely determine the influence of q 1 and q 2 on the output. The estimation problem is therefore non unique, and the Fisher information is singular.
In the first case of an irrelevant parameter, the associated singular vector in the unidentifiable parameter subspace V 2 will directly point to the parameter in question: e.g. with a singular vector (76) of it can be concluded that the second element of the parameter vector can be disregarded. In the second case, of effects of parameters appearing in combinations, the associated singular vector will have non zero entries for the parameters in question. In order to obtain a unique solution, it is advised to either fix all parameters in question but one to a constant value, or, if possible, to lump the parameters together. If, with the described measures, the unidentifiable subspace is reduced to zero, it is still possible for the Fisher information to be ill-conditioned, leading to large parameter covariances. The steps above can therefore be repeated, for singular values below a defined threshold s min .

Experimental setup
For the generation of the measurement data, a PEM-fuel cell test bench was set up, which is shown in Fig. 3. The test bench consists of a 30 kW PEM fuel cell stack (96 cells in series, 409 cm 2 cell area), an anode module, the cooling circuits, the air supply system and the control cabinet. The fuel cell stack is connected to a battery simulator (dynamic DC/AC converter) for the load set point control. During operation, the electric fuel cell stack power is fed into the electric grid via the battery simulator. The balance of plant components (BOP e air compressor, pumps, control cabinet, etc.) are supplied with energy via the electric grid.

Media supply
The reactants of the stack are supplied by the anode module (hydrogen) and air supply system. On the hydrogen side, the gas bottle pressure is reduced and supplied to the anode module. Another pressure controller guarantees a further reduced constant hydrogen pressure inside the anode. A hydrogen recirculation pump provides a constant volume flow through the anode. No external humidification is used for either the hydrogen or air side (internal humidification inside the stack). In addition to the hydrogen purge valve in the anode module, a water separator is used to periodically remove excess water from the hydrogen gas. On the air side, an air filter prevents dust and damaging gases from entering the air supply system. A turbo compressor pumps the air through an intercooler to the fuel cell stack (no external humidification of the reactant air). The intake air has a constant temperature of 20 C and a constant relative humidity of 45 %. After the stack, the electronically controlled air throttle valve can be utilized to vary the backpressure of the stack.

Cooling circuits
The stack, the air intercooler and the air compressor are liquid cooled. All the other components are passively cooled. For the stack cooling, a cooling circuit with deionized cooling fluid is used (constant speed and volume flow of the coolant pump). Tap water is used to cool the heat exchanger of the stack cooling circuit, the intercooler, and the air compressor. For the temperature control of each component, a flow control valve on the tap water side is utilized. The target stack coolant inlet temperature and air inlet temperature is 55 C, the compressor uses a constant coolant volume flow at all times, which suffices for the necessary cooling demand.

Test bench control
The fuel cell system test bench control software was specifically developed in LabView for full control and flexibility of each operating parameter. During standard operation, all operating parameters are set constant (target temperatures, hydrogen pressure, hydrogen recirculation pump speed, coolant pump speed) except the air mass flow, for which a constant stoichiometric factor of 1.5 is set. Numerous sensors are integrated into the system at all relevant positions (temperature, pressure, mass flow, voltage, current, and humidity sensors) to allow for a comprehensive measurement analysis. A more detailed description of the test bench set up is given in Refs. [44,45], where steady state experiments were performed.

Experimental tests and operating conditions
The experimental tests were performed by varying the operating parameters within a range where a stable operation is still feasible. In this manner, measurement data under different operating conditions was obtained, and the effects of individual parameter changes were recorded. The standard operating conditions are shown in Table 1. The varied parameters are electric current, air mass flow (which also changes the cathode supply manifold pressure), and anode pressure. The step size of each variation is randomly set and low enough not to interfere with the proper operation of system. A more detailed description of the experimental tests is given together with simulated and measured data in Section Results. Fig. 4 shows the measurement data of three local experiments and the corresponding fit of the local linear structured state space models. The results were obtained by applying the method described in Section Simultaneous local linear structured state space estimation. First, suitable measurement data is analyzed to determine local linear operating points. Second, the analytical, linear state-space system is derived. Third, a plausible initial parameter vector b q 0 is chosen based on literature values and expert's knowledge. Finally, the proposed parametrization problem is optimized with MAT-LAB's gradient-based optimizer fmincon. This procedure is, of course, iterative. The outputs to which the model has been fitted is the stack voltage, the cathode supply manifold pressure, and the anode exit manifold pressure. The fuel cell stack was excited by applying sequential steps to the current, the air mass flow, and the reference value of the controlled anode pressure. The current and the air mass flow signals are additionally given in Figs. 4 and 5 for reference. In general, a good  i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x agreement between the local models and the experimental data is evident. However, in the third experiment, the goodness of fit between measured and fitted voltage is slightly less. One possible explanation could be that the effects of liquid water, whose presence has been experimentally observed at the exhaust port of the fuel cell stack, are currently not present in the model, i.e. the reduction of diffusivity in the GDL due to flooding. The cathode pressure behavior, excited by changing the mass flow of inlet air, is sufficiently reproduced for all experiments. The periodic pressure spikes in the anode exit manifold result from purging. It can be seen that changes in the cathode pressure have a significant effect on the stack voltage, whereas the change in anode pressure of 200 mbar (second experiment at 50 s) has no evident effect on the stack voltage. An apparent advantage of the estimation approach via multiple local linear models lies in the significant increase of computational efficiency. Using a standard personal laptop (Intel i7 CPU @ 2.60 GHz), the required computation time to simulate the non-linear model for the 120 s of measurement data in the first experiment of Fig. 4, amounts to 30.1 s. Thereby, MATLABs ODE15s solver for stiff systems is used with the largest possible tolerance as to produce a stable solution. However, the simulation of the derived linear structured state space model for the same dataset requires 0.0281 s. As the simulation is embedded into the optimization, this decrease in computation time of three orders of magnitude directly leads to a proportional decrease of the overall time required for the optimization.

Results
After the parameter estimation via local linear structured state space estimation, the resulting parameter vector is inserted into the original non-linear model. To cross-validate the resulting non-linear model, a global experiment has been conducted. Thereby, the stack current is increased from 0 A to 300 A in a sequence of rate-limited current steps. The air mass flow throughout the experiment is constant and set such that the oxygen stoichiometry at the maximum current is l O2 ¼ 1:1. During experiment, the anode set point pressure is changed from 400 mbarG to 625 mbarG. A comparison between measurement and parameterized model is shown in Fig. 5 (a). The corresponding internal states are visible in Fig. 5 (b). The oxygen mass inside the cathode decreases with increasing current. At maximum current of 300 A (t ¼ 135 s) the oxygen mass inside the cathode is almost reduced to zero (as it is expected due to the low stoichiometry) almost leading to oxygen starvation. The displacement of hydrogen in the anode due to water vapor is visible. The accumulated nitrogen in the anode is discharged when a purging event takes place (at t ¼ 100 s and t ¼ 142 s). The membrane humidity increases with increasing current due to the generation of water vapor. Liquid water starts to form in the cathode, and the membrane water activity is saturated. When reducing the current, the liquid water evaporates, the membrane humidity decreases, and oxygen mass is increased.

Identifiability analysis
In order to assess the identifiability of parameters, the analysis presented in section Analytic computation of the parametric sensitivity for structured state space models and Fisher information based identifiability analysis has been conducted for the full set of model parameters, on the datasets seen in Fig. 4. The singular values of the Fisher information matrices, calculated for the increasing number of experiments, is shown in Fig. 6. The horizontal red line denotes the Fig. 4 e Fit of the local linear structured state space models to their respective measurement data. The columns represent the local experiments and the rows the measured output signals, to which the local models are fitted. The current and the air mass flow are given for reference.
tolerance, based on the machine precision, for which all values below are treated as effectively zero. It can be seen that the lowest singular value stays significantly below the tolerance for all experiments. Inspecting the corresponding rightsingular vector revealed that it unambiguously points in the direction of the anode supply manifold volume V an sm . The issue that the parameter V an sm cannot be identified is obvious from the fact that the anode side of the model is pressure driven, as discussed in section Anode. By switching from mass-flowdriven to pressure-driven, the differential equation coupling the inflow to pressure in the anode supply manifold is omitted, and as such V an sm has no influence on the model. Although obvious in this example, problems like this arise naturally in complex, high-dimensional models with many parameters, and it highlights the benefits of an identifiability analysis. The second lowest singular value, which was bumped slightly over the tolerance threshold by including additional experiments, can be associated with the parameter ε 2 in (35), describing the temperature dependence of the ionic conductivity. As all the experiments were conducted with an actively controlled stack temperature, such that the deviation of coolant outlet temperature where within 5 C, it can be concluded that the change in temperature during experiment is not sufficient as to determine ε 2 with reasonable significance.
The most significant parameters can be determined in an analogous manner. By evaluating the dominant direction in the parameter space of the singular vectors associated to the largest singular values, the following 10 most significant parameters (in descending order) have been determined: V ca , k ca , R 0 , V an , q 1,elchem , k an sm , q 2,elchem , k ca em , k ca sm , V ca em .
In general, a parameter could be nonsignificant if it is not relevant in the model (e. g. V an sm ), the model input is not exciting enough (e. g. ε 2 ), or both. High parameter significances can be interpreted in a vice versa manner. With this information, remodeling, or better design of the experiment can be considered. The parameter significances can also be interpreted from the PEMFC point of view. Focusing on optimizing the most significant parameters' real-world counterpart (e. g. V ca , k ca , and R 0 ) will improve the system output more significantly than optimizing the least significant ones. Optimizing the latter will lead to barely any change in the outputs. Therefore with the given parameter significance, engineering resources can be applied in a more targeted way.
To verify the order of significance of parameters, the identification via local linear structured state space models has been redone iteratively, with an increasing number of parameters to be optimized n q ¼ 1, …, 10. At each iteration, the additional parameter added corresponds to the order stated above. Each optimization is terminated after 50 iterations. Table 2 shows the relative objective function J rel ¼ J 50 /J 0 .
With only one parameter (V ca ) no effective improvement of the objective function can be made. Including the second parameter (k ca ) for optimization enables a sufficient degree of freedom as to match the overall pressure levels on the cathode side of the model vs. measurement data and thus results in a considerable reduction of the objective function. With inclusion of the third parameter R 0 , the DC voltage levels can be matched to the measurement data, resulting in another significant reduction of objective function. The parameters following thereafter give additional degrees of freedom as to fit the transient dynamics in the measured response data. But i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x as the errors due to transient mismatch are relative small as compared to the errors arising due to steady state errors, the improvement of the objective function is not as pronounced as with the inclusion of the first parameters. It is to be noted that the preference of reduction of the steady-state error, due to its disproportionately large magnitude as compared to the mismatch of transient responses, can lead to problems regarding convergence to local minima. This can be remedied by a suitable choice of a pre-filter for the experimental data [41], or by decoupling and separate weighting of steady-state and dynamic errors.
The parameter significances are strongly dependent on the model inputs, and therefore the conducted significance analysis is limited to the given experiments and is not generalizable. For example, with a more significant variation of the temperature during the experiment, the significance of temperature-dependent parameters will increase.
In the context of this work, the singular values of the Fisher information are used to analyze the significance of parameters with respect to given experimental data. As an outlook, it is to be noted that the inverse approach is of equal interest: How to design an experiment such that the identifiability of parameters is improved? Model based optimal design of experiment techniques are typically based on the maximization of a scalar criterion derived from the singular values of the Fisher Information Matrix [46], for which the algorithm presented in section Analytic computation of the parametric sensitivity for structured state space models can be directly utilized.

Conclusion
In this work, a parameter estimation methodology for the identification of highly non-linear, dynamic fuel cell stack models was presented. In order to increase the computational efficiency, a parameterization scheme based on analytic linearization has been conducted. The local linear parameterization was carried out with respect to three local experiments on a 30 kW stack. As an apparent advantage of the parameterization via multiple local linear models, a decrease in computation time of three orders of magnitude was observed. The proposed identifiability and parameter