Online estimation of the electrochemical impedance of polymer electrolyte membrane fuel cells using broad-band current excitation

Journal of Power Sources 405 (2018) 150–161

Contents lists available at ScienceDirect

Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour

Online estimation of the electrochemical impedance of polymer electrolyte membrane fuel cells using broad-band current excitation

T

D. Ritzbergera,∗, M. Striednigb, C. Simonb, C. Hametnera, S. Jakubeka a b

TU Wien, Getreidemarkt 9/E325, 1060, Vienna, Austria Chair of Technical Electrochemistry, Department of Chemistry and Catalysis Research Center, Technical University of Munich, D-85748, Garching, Germany

HIGHLIGHTS

impedance spectrum is described by a discrete time domain model. • Electrochemical online estimation algorithm is derived based on generalized total least squares. • An circuits are synthesized from the discrete time domain model. • Equivalent is done using experimental measurements from a 5 cm single cell. • Validation • Time resolved impedance information during transient operations are obtained. 2

ARTICLE INFO

ABSTRACT

Keywords: Electrochemical impedance spectroscopy Polymer electrolyte fuel cell Fault detection Generalized total least squares Online estimation

In this work, an online algorithm to estimate the electrochemical impedance spectrum of polymer electrolyte membrane fuel cells (PEMFCs) is presented. The electrochemical impedance, constituting the linear dynamic relation between current and voltage, is described by a discrete time filter model. The unknown model parameters are estimated from streaming time domain data using a recursive algorithm based on the generalized total least squares (GTLS) technique. This allows for flexible noise assumptions, leading to unbiased parameter estimates in the presence of noisy current and voltage measurements. Therefore, the robustness of the proposed method is significantly increased. Broad-band current excitation signals are utilized as to sufficiently stimulate the fuel cell for the online parameter estimation. The model order is selected using Akaike's information criterion (AIC), leading to a parsimonious model. Equivalent circuits are automatically generated using Foster's method of analog circuit synthesis. The applicability and benefits of the proposed algorithm regarding fuel cell diagnostics are demonstrated using measurement data of a small scale single cell (5 cm2 active area) subjected to a broad-band excitation.

1. Introduction To achieve a sustainable propulsion concept in automotive applications, polymer electrolyte fuel cells (PEMFCs) are widely recognized as a promising candidate. Significant progress has been made as to reach a wide-spread market deployment, yet some challenges still remain. One of them is the durability, which can not compete with internal combustion engines [1,2]. Whereas one focal area to increase durability constitutes the field of material sciences [3], the deployment of suitable online diagnostics, prognostics and health management systems has received a substantial amount of interest as to prolong the lifetime of PEMFC systems [4,5]. Electrochemical impedance spectroscopy (EIS) has proven itself to be an important tool in fuel cell diagnostics [6,7], as it enables insights into the internal processes of fuel cells and a variety of faults have been

∗

successfully detected and isolated. Nevertheless, several draw-backs limit its application for online diagnostics. 1.1. Drawbacks of EIS The electrochemical impedance, which describes a linear dynamic relation between current and voltage around a certain operating point of the fuel cell, is typically obtained by single sinusoidal excitations over a wide range of excitation frequencies. Especially due to excitations in the low frequency range, this procedure is lengthy and time consuming, during which the fuel cell has to be operated under stationary conditions. After the experimental data have been recorded, the obtained impedance spectrum is then interpreted using a parametric model, most commonly by the means of equivalent circuits [8–11].

Corresponding author. E-mail address: [email protected] (D. Ritzberger).

https://doi.org/10.1016/j.jpowsour.2018.08.082 Received 28 May 2018; Received in revised form 23 August 2018; Accepted 27 August 2018 Available online 29 October 2018 0378-7753/ © 2018 Elsevier B.V. All rights reserved.

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

Model parameters are found offline by using suitable optimization algorithms, minimizing the difference between the model in the frequency domain and the experimentally determined impedance. By applying this procedure periodically, fault information can be inferred from changes of equivalent circuit parameters over time [8,12]. Fast tracking of faults is thereby difficult due to the aforementioned significant experimental time and off-line parameter estimation.

emerging faults. However, since the error to be minimized considers only differences appearing one-step-ahead, it is apparent that prediction error methods emphasize high frequency model components in the estimation as an error build-up over time cannot be captured [19]. On the other hand, minimizing the simulated output error captures the error build-up over time, increasing the low frequency model fit. But even for linear models the resulting optimization problem is highly non-linear, requiring computationally expensive optimization algorithms. Initial parameter estimates are required, and due to the existence of local minima, convergence to the optimum can not be guaranteed. Therefore, output error minimization schemes are typically found in applications where high fidelity simulation models are in demand and the execution time of the parameter estimation plays a subordinate role. An in depth discussion about the addressed topics of system identification are found in Ljung's fundamental work [16].

1.2. Fuel cell diagnostics in the time domain As an alternative to EIS, several methods have been proposed that extract diagnostic features from time domain signals with the aim to significantly shorten the required experimental time. One common procedure constitutes the current interrupt technique [13]. Thereby, the electric load is disconnected from the fuel cell and the resulting voltage response is analyzed. Since the ohmic losses disappear immediately, the ohmic resistance can be extracted from the voltage jump following the disconnection of the electric load. Additionally, characteristic time constants are evident from the relaxation of the voltage back to the open circuit voltage [14]. From a system identification and data-driven modeling point of view, the current interrupt technique resembles the commonly applied step response experiment from which important system properties can be inferred. A more elaborate use of the current interrupt test can be found in Ref. [15]. Thereby, the obtained step response data are used to fit the parameters of an equivalent circuit model in the time domain. Apparent disadvantages of the current interrupt test and subsequent model identification are, that the disconnection of the load imposes a severe interference on the fuel cell system, and that due to the magnitude of the current interrupt the voltage response is most likely non-linear, e.g. the Kramers-Kronig relation is violated. Therefore, fitting equivalent circuits, whose impedance spectra are only valid in proximity around an operating point, is not entirely justified. Instead of disconnecting the load to obtain a step response, it is advantageous to superimpose a sufficiently small excitation signal around an operating point. A commonly applied excitation signal for system identification and parameter estimation constitutes the pseudo random binary sequence (PRBS), which consists of a series of steps between two fixed amplitude levels. PRBS excitations are favorable input signals for the system identification, as they are persistently exciting [16], e.g. provide a sufficient stimulation over a broad frequency range. In Ref. [17] a nonparametric identification method by the means of wavelet transformation is presented to extract the electrochemical impedance from PRBS data, having the distinct advantage that no assumptions about an underlying time domain model are necessary. However, the trade-off typically lies in the high uncertainties of nonparametric models. An extensive overview about non-model based diagnostic methods is given in Ref. [18]. Nonparametric methods play an important role in the preidentification of dynamic systems, as to obtain insights into possible model orders and structural properties. However, by a-priori making the (correct) structural assumptions, the accuracy of the obtained model may be significantly improved.

1.4. Impedance estimation via time domain system identification The aforementioned general considerations also hold in the context of PEMFC model identification. In Ref. [21], an output error based parameter identification of an equivalent circuit model, consisting of one resistance and 2 R C elements in series, is proposed. To overcome local minima and find suitable initial parameter estimates, a computational expensive hybrid global/local optimization is proposed. A good overall fit to experimental data (in the time and frequency domain) for the resulting model is obtained. In Ref. [22], a similar identification procedure is proposed. The chosen equivalent circuit model includes transport phenomena which are modeled by a second order approximation of the Warburg impedance, leading to a wide spread of characteristic times. Due to the chosen output error approach, satisfactory results over the whole frequency range are obtained. In Ref. [23], a prediction error based, nonrecursive least squares method is applied for the estimation of a low order model. Due to the high frequency emphasis of the estimation approach, mass transport phenomena are not captured. However, the parameter estimates obtained correspond well to changing operating conditions with significantly less computational effort, highlighting the applicability for fault diagnostics. Further discussion about advantages and disadvantages of linear vs. non-linear optimization (e.g. prediction error vs. output error) for model-based PEM fault detection methods is found in Ref. [24]. For the online impedance estimation of lithium-ion batteries, prediction error based recursive least squares estimation schemes have been proposed in Refs. [25–27]. Thereby, fast response time and tracking of system parameters during transient operations is confirmed. A main limitation regarding the practical application of least squares based online impedance estimation schemes constitutes the fact, that in the presence of significant noise in voltage and/or current measurements, the obtained parameter estimates become biased [28]. The reason for this is, that the noise assumption that leads to the convenient least squares estimation, is inherently incorrect for dynamic, autoregressive systems. Thereby it is assumed that only the independent variable (output vector) is affected by noise where as the dependent variables (regressor) are noise-free. But since for discrete time domain models past (noisy) outputs are also inputs to the system, the least squares assumption is always violated and biased parameter estimates are the result, which is a well-known phenomenon [20,29]. To obtain un-biased parameter estimates in the presence of noisy model inputs, the total least squares (TLS) algorithm has been proposed and is thoroughly discussed in Ref. [30]. As an extension, the generalized total least squares (GTLS) addresses the problem when some inputs are affected by noise, whereas some are noise-free [31–33]. This is a useful assumption in many practical cases. Moreover, the LS and TLS solution are obtained as special cases of the GTLS estimation procedure.

1.3. Prediction error vs. output error estimation techniques The fitting of the time domain model to the measurement data can be achieved by either considering the one-step-ahead prediction error, or the simulation output error [19]. Depending on the intended use of the model the following considerations can give guidelines for an appropriate choice: For linear model structures (in combination with certain noise assumptions), the minimization of the one-step-ahead prediction error leads to the least squares (LS) estimation problem. Inherent advantages are its computational efficiency and that no initial estimate is required. Owing to the simplicity and the availability of recursive algorithms, estimations based on the one-step-ahead prediction error are commonly applied for online parameter estimation schemes for fault detection [20]. Parameter updates are obtained at each time step leading to a fast reaction to

1.5. Proposed methodology As a main contribution of this work, the GTLS estimation is reformulated as a recursive algorithm suitable for online equivalent 151

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

circuit parameter estimation. Special considerations regarding the computational efficiency are thereby made. Due to the more flexible noise assumptions, the estimation bias can be avoided as opposed to the commonly found recursive least squares algorithm. The superimposed excitation signal used is an amplitude-modulated pseudo random binary sequency (APRBS). The advantage of amplitude modulation lies in the fact that non-linearities in the output signal, which are to be avoided, are more easily detectable. Another key aspect in this work is the data-driven synthesis of equivalent circuits. Typically, the common workflow comprises the apriori selection of an equivalent circuit model [21–24,28]. Afterwards, the model is discretized to enable the fitting to the discrete time measurement data. The relations between estimated discrete model parameters and desired physical equivalent circuit parameters need to be derived analytically. It is obvious, that for increasing complexity of the equivalent circuit models, this analytic derivation becomes tedious and hinders the flexibility. In this work, an alternative approach is proposed, where the starting point is an arbitrary discrete time difference equation with a pre-defined order. For the order selection, statistical tools commonly found in the area of data-driven modeling and system identification, such as Akaike's information criterion (AIC) [34,35] are used, leading to a parsimonious model. Ideally, such a model only consists of model terms which are significant for the present time-domain data. After obtaining the parameter estimate, Foster's method [36] is applied, which synthesizes an equivalent circuit whose impedance matches that of the estimated model. The remainder of the paper is structured as follows: At first the methodology of prediction error based impedance estimation in the time domain is discussed. The occurance of parameter estimation bias due to measurement noise for the standard LS approach is shown, motivating the TLS and GTLS estimation techniques. Afterwards, the GTLS algorithm is adapted to a recursive and computationally efficient online algorithm. The topic of parsimonious model order selection via AIC and subsequent equivalent circuit synthesis is addressed. In the results section, statistical analysis of the proposed estimation algorithm in the context of impedance estimation is done using simulated equivalent circuit data. Following this, experimental fuel cell data of a 5 cm2 single cell under transient operating conditions are analyzed.

synthesizing an analogue circuit from (1) proposed in this work is discussed in section 2.4. 2.2. Parameter estimation in the time domain Rewriting Eqn. (1) in Matrix form for N data samples leads to with

Y=

xT (k ) = [

… an y (k

n) + b0 u (k ) + b1 u (k

n),

(4)

y (k

1) …

y (k

n ) u (k ) u (k

1) … u (k

n)].

(6)

˜ = Y + Y, Y the least squares parameter estimate is obtained as

ˆLS = arg min Y ˜

X

2 2

˜. = (XT X) 1XT Y

(7)

In the remainder of this work, a tilde ˜ is used to denote vectors and matrices subjected to Gaussian noise, and a caret ˆ for estimated/reconstructed quantities. From (5) it can be directly seen, that by assuming noisy output measurements y˜ (k ) , the regressor is inherently corrupted by noise as well,

˜ = X + X, X

(8)

as each row contains lagged output measurements. This directly violates the assumption of a noise-free regressor and biased estimates are obtained as a result:

= =

˜ TX ˜ ) 1X ˜ TY ˜} {(X T ˜ X ˜ ) 1X ˜ T (X ˜ + Y)} {(X T ˜ 1˜ T ˜ + {(X X) X Y}.

(9)

Thereby, { } denotes the expected value operator. Note, that if the least squares assumption is valid, e.g. X = 0 , the second term of the right hand side (9) vanishes, as { Y} = 0 and unbiased estimates are obtained. For X 0 , this is in general not the case. The occurance of biased equivalent circuit parameter estimates in the presence of measurement noise when applying RLS techniques for the online impedance estimation of lithium-ion batteries has been emphasized in Ref. [28]. To address this issue, the proposed online estimation procedure in this work relies on the more flexible generalized total least squares estimation technique.

1) (1)

with the unknown parameter vector

= [ a1 … an b0 b1 … bn ]T .

xT (N )

(5)

It is assumed that the electrochemical impedance to be estimated is derived from a purely linear system and that the Kramers-Kronig relation [37] is valid. This justifies the use of linear dynamic models, which are given in the discrete time domain as the following difference equation:

1)

.

Assuming that the dependent vector is subject to Gaussian noise with zero mean and equal variance,

{ ˆLS} =

+ …+ bn u (k

,X=

xT (n + 1) xT (n + 2)

Thereby, the regressor xT (k ) is given as

2.1. Discrete time domain model of PEMFC impedance

a1 y (k

y (n + 1) y (n + 2) y (N )

2. Methodology

y (k ) =

(3)

Y=X ,

(2)

2.3. Generalized total least squares

Thereby, y (k ) denotes the discrete model output, e.g. voltage, and u (k ) the model input, e.g. current. The index k is used as the discrete time index according to t = k t with t beeing the sampling time. The model order n defines maximum number of lagged inputs and outputs considered. Eqn. (1) is a valid discrete time domain representation of any equivalent circuit model, consisting of a network of resistors, capacitors and inductances, with arbitrary structure. However, the effects of constant phase elements can thereby not be modeled, as they would lead to fractional order derivates in the time domain [38,39]. Commonly applied equivalent circuits and their corresponding difference equations can be found in Refs. [21–24,28]. The alternative approach of

In the case, that both the output vector and regressor matrix are subject to noise, the TLS algorithm [30] is the proper extension to the LS estimation problem (7). The more general case, that only some of the columns of the regressor matrix are subjected to noise, is considered by the GTLS problem. A practical application for this assumption is, for example, that the measurement noise variance of the input signal is several orders of magnitude smaller than the output measurement noise variance. Therefore a noise decorrelation can easily become ill-conditioned. 152

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

2.3.1. Parameter estimation from batch data with TLS and GTLS N ×1 With the deterministic output vector Y and deterministic N ×M regressor matrix X the linear model (3) can be written as

[X Y]

1

W

ˆ Tn (k ) = wTn (k ) w

Wn

||Wn

ˆ ||2F W

(12)

ˆ ). Image(X

ˆ n||2F , W

ˆ ). Image(X

1mT

WoH) b.

pn + po + 1 × 1

(22)

as

T

subject to

T

,

(24) (25)

= 1,

1T Wn WTo Wn

WTn 1

WTn Wo

N

1T Wo

WTo 1

WTo Wo

, (26)

and

(13)

=

(14)

I pn × pn 0 pn × po + 1 . 0 po + 1 × pn 0 po + 1 × po + 1

(27)

The equality constraint (25) ensures that ||b||2 = 1. Thereby, a non-trivial solution of (24) is obtained. From the constrained optimization problem (24) and (25), a Lagrange function is formed

L( , ) =

T

(

T

L( , )

(15)

=(

(28)

1).

Setting the partial derivatives with respect to

) = 0.

to zero leads to (29)

Note, that (29) is a generalized eigenvalue problem. As a solution, is the generalized eigenvector corresponding to the smallest eigenvalue λ. Additionally, the equality constraint (25) can always be satisfied by appropriate scaling of the eigenvector. Parameters corresponding to noisy regressors are then obtained by

(16)

ˆGTLS, n =

1:(pn 1)

,

(30)

pn

(17)

whereas parameters corresponding to noise free regressors are obtained by

ˆGTLS, o =

(pn + 2):(pn + po + 1) pn

.

(31)

Thereby, the index notation i :j denotes the sub-vector of from the i th to the j th element. The complete parameter vector ˆGTLS is then obtained by combining and rearranging corresponding to the chosen partition (17). Note, that if the regressor X is free of error and ˜ and Wo = X , the GTLS estimate is equivalent to the LS therefore Wn = Y ˜, Y ˜ ] and Wo = , the total least squares solution (7), whereas if Wn = [X solution is obtained.

(18)

subject to

Yˆ

minimize

=

where Wn denotes the pn × N sub-matrix which corresponds to columns subjected to uncorrelated noise with zero mean and equal variance. Subsequently, Wo denotes the po × N sub-matrix of regressor elements which are noise-free. The previously introduced tilde-notation ˜ for vectors and matrices subjected to noise is omitted from here on due to redundancy. For the GTLS parameter estimate, only the noisy components of the augmented regressor are to be reconstructed:

min||Wn

WoH)T (Wn

(23)

WTn Wn

Note, that (15) can be seen as the extension of the LS Problem, with the ˜ is addition that in order to find a parameter estimate ˆ , the regressor X subject to reconstruction as well. For the GTLS algorithm [33,41], W is rearranged as

˜ = [ Wn Wo ], W

1mT

with

subject to

Yˆ

ˆ n ||2F = bT (Wn W

the minimization of (22) can be written as an equality constrained quadratic program

˜ in Thereby, ˆ = diag( 1, 2, … M , 0) contains the singular values of W descending order, with the most insignificant singular value set to zero, and U and V consist of the corresponding left and right singular vectors. The parameter vector obtained from the kernel base of the apˆ , is called the total least squares (TLS) solution of proximated matrix W the formally stated minimization problem:

˜ min||W

N × 1.

b = mT b , Hb

(11)

is given by

ˆ = U ˆ VT . W

(21)

WoH) b] bT ,

By defining a combined parameter vector

˜ is is generally of full rank. According to the considerations above, W ˆ using the Eckarttherefore approximated by a rank deficient Matrix W Young-Mirsky Theorem [40]: ˜ , that minimizes the Frobenius The best rank M approximation of W norm

ˆ ||2F , W

1mT

With 1=[1 1 … 1] Inserting (21) into (18) leads to

(10)

where Ker( ) denotes the Kernel (or null space) of the augmented matrix N × (M + 1) W . A base for Ker(W) is given via singular value decomposition (SVD) as the right singular vector belonging to the singular value that equals zero. By appropriate scaling of this base vector such that the last component equals 1, the augmented parameter vector is obtained from which can be extracted. When subject to noise, the matrix

˜ min||W

ˆ n = [(Wn W

T

where N denotes the number of recorded observations and M the number of parameters to be estimated. Since Y Image(X) , it follows that Ker(W) , e.g. W is rank deficient. Specifically,

˜ = [ X + X Y + Y ], W

(20)

HT wo (k )T ) b] bT .

m

Aggregated for all N data records, this leads to

= 0,

= Ker(W),

[(wn (k )

(19)

2.3.2. Recursive generalized total least squares While it is shown above how to obtain the GTLS parameter estimate from a batch of recorded input and output data, it is now discussed how to adapt the GTLS solution as a recursive algorithm suitable for online parameter estimation. Special consideration about the computational

For the reconstruction, the noisy components Wn are projected onto pn × 1, a unit normal a hyperplane defined by a reference point m p × 1 n vector b and a linear combination of the noise-free regressors po × pn po × 1 given by HT wo (k ) with H and wo (k ) , which is a row vector of Wo : 153

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

efficiency are thereby necessary. Assume that from the previous time step the matrix (k 1) and the corresponding generalized eigenvector (k 1) , belonging to the smallest eigenvalue, are known. At the current time instance k, new measurements of the input u (k ) and output y (k ) are available. After updating the current regressor vector xT (k ) , the augmented regressor vector is formed as

wT (k ) = [ wTn (k ) wTo (k )].

historic data. One commonly used approach is the rectangular receding window. Thereby, historic data at a time instance j = k L is discarded from the estimation, where L is the data-length of the rectangular window, which serves as tuning parameter for the inherent tradeoff between tracking performance and covariance of the parameter estimate. Analogous to the addition of new data to the estimation, the current augmented regressor matrices can be partitioned as

(32)

It can be seen that by integrating the new measurements as

Wn (k ) =

Wn (k

1)

,

wTn (k )

Wn (k ) =

Wo (k ) = Wo (k

1)

,

wTo (k )

(k

1) +

(34)

(k )

(35)

T (k ) ,

with

¯ (k ) =

wn (k ) . 1 wo (k )

(k ) =

¯

(k )

(j )

(44)

T (j ),

One possibility is the shift-and-invert power iteration [42]: By introducing the auxiliary matrix

(k ) =

1 (k )

1 (k ) T (j )

(j ) T (j ) 1 (k ) (j )

1 (k )

1

,

(45)

(37)

1

)(k )

1

which is then to be used in the auxiliary matrix (37) for the inverse power iteration. The recursive GTLS algorithm with sliding rectangular window is summarized as follows:

(36)

(k ) ) (k ) = 0.

S(k ) = ( )

(43)

and by applying the Woodbury matrix identity as to avoid the matrix inversion, the following recursive forgetting of historic data is obtained

The next step is to solve for the desired eigenvector (k ) corresponding to the smallest eigenvalue (k ) of the time-dependent generalized eigenvalue problem

( (k )

wTo (j) . ¯ o (k ) W

¯ n (k ) and W ¯ o (k ) are the current noisy and noise-free subThereby W matrices excluding the data at time instance j from the estimation. Note that all data preceding j have already been excluded prior to the current time-step k. The influence of historic data is removed from the estimation by

matrix (26) can be recursively updated as

(k ) =

(42)

and

(33)

and

Wo (k ) =

wTn (j) , ¯ n (k ) W

a sequence {qi + 1 (k )}i = 0 generated by (38)

pi + 1 (k ) = S (k ) q i (k ), q i + 1 (k ) =

pi + 1 (k ) ||pi + 1 (k )||2

,

(39)

converges to a generalized eigenvector (k ) corresponding to an eigenvalue in the vicinity of σ. Since (k ) is positive definite, convergence to the eigenvector corresponding to the smallest eigenvalue is achieved by setting = 0 . The iteration is stopped after a suitable convergence criterion is satisfied, such as

ql (k ), q l 1 (k )

1

(40)

,

(k ) with sufficient accuracy. The after which it is assumed that q l number of l iterations required depends on the chosen initial vector from which the power iteration is started. An appropriate initial vector constitutes the eigenvector from the previous time-step (k 1) . Since ˆGTLS (k ) is determined from (k ) , the recursive update of (35) and the inverse power iteration (38)–(39) already constitute a recursive GTLS algorithm. However, one major drawback is the computationally expensive matrix inversion in (37) that has to be carried out at each time-step. This can be avoided by using the Woodbury matrix identity [43] to update the inverse of (k ) directly: 1 (k )

=

1 (k

1)

1 (k

1+

1) (k ) T (k ) 1 (k 1) . T (k ) 1 (k 1) (k )

The condition for past data to be removed is deliberately left ambiguous in the algorithm above. While a rectangular window of constant length is easily obtained, it is interesting to note, that a variable window length can be implemented as-well, which may improve the convergence to time-varying parameters, see Ref. [44].

(41)

Note, that the denominator in (41) is scalar. 2.3.3. Forgetting of historic data To improve the tracking capability in the online parameter estimation of time-varying systems, it is necessary to limit the influence of

2.4. Algorithmic synthesis of equivalent circuit models When analyzing EIS data in the frequency domain, the choice of the 154

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

equivalent circuit model is of crucial importance. For exceptional EIS data quality it is possible to use equivalent circuits with increasing complexity to model many aspects of dynamic PEMFC behavior. However, one critical aspect of model selection, that is sometimes paid too little attention, is, that an increase of the model order, and therefore an increase of the available degrees of freedom, will always lead to better fit of the model to the given data set used in the estimation. But at some point, further increasing the model order will lead to an overfit, as the model tries to explain random disturbances and measurement noise present in the data. Nevertheless, the goodness-of-fit is increasing giving the false impression of increased accuracy. This circumstance is known as the bias-variance trade off [16]. A model that balances this trade-off, in the sense that there are enough parameters as to capture relevant system properties, but not too many as to produce an over-fit, is called parsimonious model. For impedance fitting in the frequency domain, a popular technique that is applied in the sense of finding a parsimonious model, is the distribution of relaxation times (DRT) analysis [45,46] to detect significant time constants. For the impedance estimation in the time domain, the equivalent circuit models are typically selected a-priori. Certain model components which are commonly found in the frequency domain fitting are simplified such that the time domain representation only consists of ordinary differential equations. For example the Warburg impedance can be approximated by a finite number of RC elements [15], and constant phase elements are replaced by a capacitor [24]. After discretizing the chosen equivalent circuit model, a difference equation (1) with a specific order n is obtained. The relation between discrete model parameter vector (2) and physical equivalent circuit parameters needs to be derived. Commonly applied time domain models can be found in Refs. [21–24,28] and can be readily applied for the online parameter estimation with the proposed recursive GTLS algorithm. As an alternative to the a-priori model selection, an inverse approach is proposed to obtain equivalent circuits from an arbitrary difference equation (1). At the first step, the model order n has to be selected. As discussed, it is of importance to find a parsimonious model order. Many methods have been proposed that can be utilized for finding a parsimonious model with respect to a given set of data, such as Akaike's information criterion [34]:

AIC = N ln(

˜ ||Y

ˆ (n)||22 Y ) + 2(2n + 1). N

sampling frequency and the excitation signal used. Therefore, the model order with a minimum AIC value does not necessarily reveal the ”true” order of the underlying system, but it gives a decision on the number of parameters that are significant in reproducing the given data. Finding a concise yet significant set of parameters is especially useful for fault detection as to avoid redundancies and ambiguities in the diagnostic signals. After embedding the discrete filter model with selected order into Algorithm 1, and obtaining numerical values for the estimated parameter vector (2), the next step is to synthesize an equivalent circuit whose impedance matches that of the difference equation. For that, (1) is transformed into the discrete frequency domain by applying the Ztransformation:

G (z ) =

b0 z n + b1 z n 1+…+ bn . z n + a1 z n 1+…+an

(47)

Since electric circuit synthesis techniques are typically formulated in the continuous frequency domain, (47) is further transformed using a suitable inverse discretization method, such as the bilinear transformation (Tustin-method) [47], leading to

Gc (s ) =

bc,0 s n + bc,1 s n 1+…+ bc, n . s n + ac ,1 s n 1+…+an

(48)

Note, that due to the transformation from discrete to continuous frequency domain, the parameters in (47) and (48) differ. A suitable equivalent circuit synthesis technique constitutes the Foster synthesis [36,48], which deals with the problem of finding a passive electrical network, corresponding to the given impedance as a transfer function (48). Applying a partial fraction expansion to (48) leads to:

ri

Gc (s ) = i

s

pi

+ k (s ).

(49)

The residuals ri and poles pi are either real valued, complex conjugates, or purely imaginary. Thus, the possible number of combinations is thereby limited. Each of these fundamental partial fractions now corresponds to a fundamental equivalent circuit, which can be seen in Fig. 1. The algebraic relation between residuals and poles with respect to the physical parameter values of the equivalent circuit are given in Table 1. The complete equivalent circuit is then formed as a serial connection of the fundamental circuits. Note that the order n of the transfer function correlates with the number of occurring fundamental equivalent circuits, however, which specific fundamental equivalent circuit is present after synthesis is solely determined by the location of the poles and zeros of the transfer function (48) on the complex plane. A vague idea about the expected poles and zeros can be inferred from a step response. When applying a direct current step, an overdamped dynamic settling of the voltage indicates the presence of purely real valued poles and zeros, leading to an equivalent circuit that

(46)

Thereby, N denotes the number of samples in a data set, n is the discrete filter order and Yˆ (n) is the predicted output obtained from an n-th order model. The idea is, that several models with increasing order, e.g. n [0, …, 10] are estimated from a data batch and the corresponding AIC values are compared. A minimum AIC gives an indication about an appropriate model order, balancing the goodness of fit and number of parameters. Note that since this test is depending on a data set, different decisions can be made for the same physical system, depending on the

Fig. 1. Fundamental partial fractions of the transfer function and equivalent Foster circuit elements. An asterisk ()* denotes the complex conjugate. 155

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

Commercial gas diffusion layers (H14C10, Freudenberg & Co. KG) were used at cathode and anode side for all tests. The gas diffusion layer compression was set to 20% by selecting a proper gasket thickness.

Table 1 Relation between residuals and poles of partial fractions and the physical parameters of corresponding equivalent circuits (cf. Fig. 1). Nr.

Residues & Poles

1) 2)

r r

,p= ,p

Physical Parameters

r

4)

r= +i p= +i

r , p

R=

C= 3)

2.7. Testing procedure

R=r

,p=0

C= ,

After a conditioning procedure [53], two different load points (1.5 Acm−2 and 2.0 Acm−2) were investigated at humid conditions (80 °C, 100% RH, 300 kPa abs ) with a flow rate of 2000 nccm (normal cubic centimeters per minute) H2 on the anode and 5000 nccm air on the cathode. After reaching steady state (5 min), electrochemical impedance spectroscopy (from 100 kHz to 1 Hz) and APRBS excitation with a sampling frequency of 1 kHz and 10 kHz were conducted. Finally, an APRBS excitation during a transient load step was investigated by analyzing a current jump from 1.5 Acm−2 to 2 Acm−2.

1 r 1 r

a0 = 2( + ) , a1 = 2 , b0 = |p|, b1 = 2 ,

a0 L, a1 a12 1 ,C= R2 = a1 b1 a0 a1 a13 L= 2 a1 b0 a0 (a1 b1 a0 )

R1 =

5) 6)

r r

,p= ,p=i

3. Results

L=r C=

1 , r

L=

2r 2

The results section is structured into two parts: At first, the statistical properties of the GTLS algorithm in the context of impedance estimation are assessed using simulation data of a simple equivalent circuit. In the second part, measurement data of a 5 cm2 active surface PEMFC are analyzed and compared to reference electrochemical impedance spectra, obtained by standard EIS. To highlight the tracking capabilities regarding transient operation, a direct current step is applied during the experiment together with the superimposed APRBS signal.

consists of fundamental circuits Nr. 1) and 2). Observing a voltage overshoot indicates complex valued poles and/or zeros leading to the presence of fundamental circuits 4) and 6). In Ref. [49], an equivalent circuit is proposed for analyzing overshoot phenomena in high temperature PEMs, that consists of a simplified randles circuit extended with an additional oscillatory circuit, which is precisely Nr. 4) in Fig. 1. Overshoot phenomena can also occur in low temperature PEMs [50], for which oscillatory circuits are proposed [51]. With the proposed methodology of parsimonious order selection and circuit synthesis, the appropriate equivalent circuit structure is automatically derived from the time-domain data.

3.1. Equivalent circuit simulation The equivalent circuit that is used to investigate the statistical properties consists of a resistor (2.304 m ) and two RC-Elements (R1 = 18.1 m , C1 = 0.587 F , R2 = 2m and C2 = 0.2 F ) in series. The applied APRBS signal is assumed to be measured without error, whereas the simulated voltage measurement is corrupted by Gaussian noise with a standard deviation of = 5 10 5 V . Exemplary time domain voltage and current data used in the estimation, including the additive Gaussian noise, can be seen in Fig. 2. The order of the discrete difference equation (1) is set to n = 2 , and for the GTLS estimation, according to (10), the first, second and last column of W are defined as the sub-matrix Wn , since only the voltage measurement is corrupted by noise. The sampling frequency is chosen to be 10 kHz .

2.5. Fuel cell test setup In order to validate the presented method, impedance spectra as well as voltage responses to APBRS current excitations were measured in a 5 cm2 single cell setup (Fuel Cell Technologies, Inc.). Custom 7 channel serpentine graphite flow fields (Poco Graphite, Inc.) were used at anode and cathode with a thermocouple in the center of the cathode block in order to measure the cell temperature. More details regarding the cell setup are shown in Ref. [52]. The cell was tested on a Greenlight Innovation G60 fuel cell test station equipped with an Agilent N3306A load bank. Electrochemical impedance spectroscopy as well as the APBRS excitation analysis were performed by a Reference 3000 potentiostat (Gamry Instruments, Inc.).

3.1.1. Analyzing the bias reduction The discrete system parameters are estimated repeatedly (5000 times) for the same batch of experimental data, but with varying Gaussian noise added to the voltage measurement. Thereby, the data set used in the estimation spans a time interval of 1 s. All 5000 realizations of the parameter estimate are inserted into the transfer function (48) and each transfer function is evaluated along a logarithmically spaced frequency grid, consisting of 130 frequency points. Thus, for each

2.6. Materials Fuel cell tests were performed with commercial catalyst coated membranes (W. L. Gore & Associates, Inc.) with electrode loadings of 0.4 mg (Pt) cm 2 on the cathode and 0.1 mg (Pt) cm 2 on the anode.

Fig. 2. Dynamic voltage response of the equivalent circuit subjected to an APRBS current excitation. The voltage measurement is corrupted by Gaussian noise. 156

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

frequency point on the impedance curve, a stochastic sample consisting of 5000 realizations is available. This now facilitates an analysis of the underlying probability distribution of real and imaginary part of the impedance at each frequency point. The assumed null hypothesis is, that the real and imaginary parts originate from a normal distribution. This hypothesis is tested by the Anderson-Daling test at a 5% significance and has not been rejected at any frequency point. Therefore, the calculation of sample mean and covariance, defined a bi-variate normal distribution of the impedance at each frequency point, is justified. From this bi-variate normal distribution, the confidence bounds [0.001, 0.1] are calculated for varying levels of confidence 1 , [54]. Each confidence bound is an ellipsis centered around its sample mean, which should ideally coincide with the actual impedance of the system. The results of the analysis above can be seen in Fig. 3. Even under the moderate noise conditions as depicted in Fig. 2, using the least squares algorithm for the impedance estimation leads to biased parameter estimates due to the fact that the noise assumption is inherently incorrect. Subsequently, the sample mean of the estimated impedance is off-center with respect the analytic impedance. Using the GTLS estimation with the correct partitioning into Wn and Wo results in a significant reduction of the estimation bias. However, the covariance is slightly increased in the low frequency domain of the impedance.

several orders of magnitude. With the algorithm settings as used above, mass transport phenomena are difficult to estimate, as the sampling frequency of 10 kHz leads to a drastic over-sampling of mass transport phenomena with time constants of below 1 s [55]. To improve the low frequency accuracy, the following options are available:

3.1.2. Influence of algorithm settings and inherent trade-offs Besides the appropriate choice of the estimation algorithm, the overall accuracy is depending on a variety of influencing factors. The information about the system dynamics in the output data is bounded by the Nyquist frequency and the sliding data-window length. Choosing the sampling frequency too low leads to aliasing and folding of the input and output spectrum, reducing the high frequency fit of the estimated model. On the other hand, oversampling typically leads to a decrease of the low frequency fit, since one-step-ahead prediction error methods emphasize high frequency model components. Drastic oversampling additionally leads to numerical problems as all the poles of the discrete transfer function (47) are converging to the same point on the unit circle when increasing the sampling frequency significantly above the eigenfrequencies of the system. Increasing the sliding data window length improves the estimation by reducing the parameter covariance and enabling the estimation of low frequency model components, but the computational burden increases and the convergence speed of the parameters to changing operating conditions decreases. Another key factor in obtaining sufficiently accurate parameter estimates is the choice of the applied excitation signal as the input spectrum directly influences the bias distribution of the estimator [16]. It is obvious, that due to these conflicting goals, certain trade-offs have to be made. In the context of time-domain PEMFC impedance estimation, they become especially evident in the presence of mass transport phenomena, as thereby the spread of characteristic time constants reaches

The experimental results are divided into two parts. At first, the initial order selection and necessary initializations for the recursive algorithm are carried out using batch data. Additionally, the equivalent circuit synthesis is illustrated using the numerical parameter estimates obtained from this initial estimate. Afterwards, the recursive GTLS algorithm is applied to obtain time resolved impedance estimates for the experimental data, which is subjected to transient operating conditions.

1. 2. 3. 4.

Increasing the data window length used in the estimation. Lowering the sampling frequency. Lowering the bandwidth of the excitation signal. Using lowpass/bandpass pre-filters for the time-domain input and output data.

However, the low-frequency accuracy is only improved by either sacrificing high frequency accuracy (due to aliasing or shifting the bias distribution) or by increasing the computational complexity which, at some point, interferes with real-time computability. Alternatively, the low frequency accuracy can be improved by modifying the algorithm, such as running a dual recursive GTLS configuration at different sampling frequencies (as proposed for LS impedance estimation [56]) or by switching experiments and imposing the obtained low frequency impedance as equality constraints on the GTLS estimation [41]. 3.2. Experimental results

3.2.1. Order selection, initialization, and equivalent circuit synthesis To select an appropriate model order for the present experimental data, Akaike's information criterion (46) is evaluated. A comparison of AIC values for increasing model orders n is shown in Fig. 4. According to (1), the model with n = 0 is a purely static relation between voltage and current defined by a single parameter b0 . A sharp drop in AIC value is observed by increasing the model order to n = 1, e.g. by considering first order dynamics, and the minimum is found at n = 2 . The parsimonious model candidate is therefore chosen with a model order of n = 2 . Afterwards, an initial parameter estimate from the batch data is calculated from (26)–(31). The results are used to start the recursive 1 (k 1) , (k 1) need to be supplied. In algorithm, as the matrices Fig. 5 this initial impedance estimate is compared to the reference impedance spectrum obtained through standard single sinusoidal EIS. Thereby, a remarkable characteristic of the time domain estimation in combination with broad-band excitation is evident: Whereas the reference impedance required 180 s to be recorded, the initial GTLS

Fig. 3. True impedance of the equivalent circuit compared to the distribution of the estimated impedance at selected frequencies. a) LS estimation leads to a strong parameter bias in the presence of noisy voltage measurements. b) Applying the GTLS algorithm significantly reduces the estimation bias due to the correct noise assumption. 157

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

Fig. 6. Equivalent circuit obtained by the Foster synthesis.

transfer function is obtained:

Gc (s ) =

0.0056s 2 + 36.6s + 2.145 10 4 . s 2 + 4955s + 6.112 105

(51)

Carrying out the partial fraction expansion leads to

Gc (s ) = 0.0056 +

8.3356 3.5147 + , s + 4.8286·103 s + 126.6

(52)

The first term in (52) corresponds to element 1) in Fig. 1, wheres the second and third term correspond to element 2) as both poles and residues are real valued. The resulting equivalent circuit can be seen in Fig. 6. The physical parameter values are obtained as Rs = 5.6 m , R1 = 1.7 m , C1 = 0.12 F, R2 = 27.8 m and C2 = 0.28 F by using the algebraic relations in Table 1.

Fig. 4. Comparison of AIC values for increasing model order. A minimum AIC value is obtained at order n = 2.

3.2.2. Discussion on the synthesized equivalent circuit The proposed method of model order selection, GTLS parameter estimation and circuit synthesis is purely exercised from a data-driven, black-box modeling point of view. No physical, first-principle considerations about the fuel cell system are thereby necessary. On the one hand, this increases the flexibility as opposed to an a-priori model selection which can easily lead to an ill-conditioned estimation problem when not chosen suitable for the present data, on the other hand, the interpretation of the obtained results and connection to physical processes of the fuel cell is ambiguous and far from trivial. However, with the proposed procedure a significant set of parameters is automatically extracted from streaming time-domain voltage and current data, suitable for a subsequent fault detection and isolation scheme. 3.2.3. Tracking of impedance during transient operation In Fig. 7 the experimental data used for the recursive estimation is shown. The voltage response due to the superimposed APRBS excitation is visible throughout the data. At 12 s, the direct current load undergoes a step-wise change from 1.5 Acm 2 to 2.0 Acm 2 . The settling to the new operating conditions requires several seconds. Since the recursive estimation is linked to the sampling frequency, the temporal resolution of the obtained parameter estimates is also 10 kHz. In the lower three plots of Fig. 7, the online estimated electrochemical impedance, is compared to reference impedance spectra at 1.5 Acm 2 and 2.0 Acm 2 at specific time instances. The reference impedance spectra were obtained through standard EIS prior to this transient experiment. At the first time instance at 7 s and third time-instance at 20 s, a good agreement between the online estimated impedance and the respective reference EIS is shown. However, at the second time-instance at 13 s, the low-frequency impedance is increased as compared to the reference impedance at 2.0 Acm 2 . This is due to the fact that a steady-state operating point has not been reached yet, which is also evident in the time-voltage response after the current jump. This kind of temporally resolved impedance information during transient operation of the fuel cell is not obtainable with standard EIS procedures and highlights the advantage of the proposed online algorithm. Throughout the estimation, the poles and zeros of the estimated transfer function (48) stay real valued. Therefore, the synthesized equivalent circuit does not change its structure during the experiment.

Fig. 5. Estimated Impedance from time domain experimental data using a 1 s data batch, compared to a reference impedance spectrum. The effects of proper choice of the sampling frequency is shown (1 kHz vs. 10 kHz). For the 1 kHz measurement, a high frequency bias can be observed.

estimate of the impedance is obtained after 1 s, while still maintaining a satisfactory level of accuracy. Note, that the phenomenon of the 45° high frequency branch, visible in the cut-out of Fig. 5, is not accurately captured. In the frequency domain, transmission-line models [57] are typically employed which consist of a high order of RC-elements in series. However, in the time-domain such an increase in model order does not comply with the necessary model parsimony in order to obtain a well-conditioned estimation problem. Also shown in Fig. 5 is the effect of an improper choice of sampling frequency: If the sampling frequency is set too low, high frequency effects can not be captured due to aliasing and distortions of the discrete system around the Nyquist frequency. Increasing the sampling frequency from 1 kHz to 10 kHz significantly improves the fit of the GTLS estimate compared to the reference impedance, but also increases the computational burden in terms of required calculations per second for the online estimation, as discussed above. The numerical values of the initial parameter estimate are used to illustrate the equivalent circuit synthesis, that is carried out at each sampling instance: With the parameter estimate

ˆGTLS=[

1.6044 0.6093 0.0056

0.0080 0.0025]T .

3.2.4. Time-resolved equivalent circuit parameters and confidence intervals The time-resolved equivalent circuit parameter estimates are shown in Fig. 8. The vertical dashed lines mark the time interval during which

(50)

inserted into (47) and applying the inverse discretization, the following 158

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

Fig. 7. Tracking of the electrochemical impedance from experimental data using the online GTLS algorithm: A direct current jump from 1.5Acm 2 to 2.0Acm 2 has been applied during the experiment. The online estimated impedance is compared to the reference impedance spectra at three time instances. At t1 and t2 , good agreement between online estimate and reference impedance is evident. At t2 , an increase of low frequency impedance is visible, as a steady-state operating point has not been reached yet.

the current jump is present in the sliding estimation data window. As the data are not compatible during this time interval, heavily distorted parameter estimates are the result which are removed from the plot. The gray areas around the parameter estimates denote the 95% confidence intervals which are determined by the following hybrid analytic/numerical approach: An analytic approximation of the parameter covariance matrix is given in Ref. [32] as

(k ) = cov(

method [58]. Knowledge of the confidence intervals is invaluable for a test decision whether a significant change in parameters has occurred, issuing an alert in a subsequent fault detection scheme. However, the hybrid analytic/numeric calculations above, which have to be carried out at each time-step are computationally costly and inadequate for an online application. One possibility to obtain confidence intervals online is the use of suitable approximation techniques for the nonlinear transformation of probability distributions such as a first order Taylor series expansion of (55) or the unscented transformation [59]. The most prominent change in physical parameters due to the change of the operating condition is visible for R2 , which corresponds to the resistance of the low frequency RC loop. The relative size of the confidence intervals also reflects the emphasis on high frequency model components for the recursive GTLS algorithm, which was also evident in Fig. 3.

GTLS (k ))

(1 + || ˆGTLS (k )||22 )

2

(k )(X˜ T (k ) X˜ (k )

N

2

(k ) IG).

(53)

Thereby IG denotes an identity matrix with the diagonal elements corresponding to noise-free parameters set to zero, and 2 (k ) is the error variance approximated by 2

(k )

(k ) , N

(54)

4. Conclusion

where (k ) is the smallest eigenvalue of the generalized eigenvalue problem (36) at the current time step k. This approximation has been found to be reasonably accurate when the amount of additive noise is moderate. However, to obtain the covariance of the physical parameters, grouped in vector pˆ (k ) , it is necessary to consider the nonlinear transformation

pˆ (k ) = f ( ˆGTLS (k )),

In this work, an online estimation algorithm of the electrochemical impedance, based on the generalized total least squares algorithm is presented. Using simulated data of an equivalent circuit, the bias-reduction in the presence of noisy measurements, as compared to the least squares algorithm, has been shown. Furthermore, the presented online impedance estimation scheme has been applied to real measurement data of a 5 cm2 PEMFC, subjected to APRBS current excitations. The appropriate model order is selected by using Akaike's information criterion, leading to a parsimonious model. Equivalent circuits are automatically generated by applying Foster's analogue circuit synthesis to the parsimonious model. A good agreement between the online estimate and the reference impedance, obtained through conventional EIS, has been observed when the fuel cell is operated under stationary conditions. For the present data, satisfactory initial estimates of the impedance are already obtained after 1 s, which is a significant time reduction as compared to the EIS procedure. Moreover, the applicability during transient operations has been highlighted using measurement data subject to a step-wise change of the operating point. With the proposed online algorithm, time resolved impedance information is obtained during transient operation.

(55)

which consists of the inverse discretization and analogue circuit synthesis described in section 2.4. Due to the complexity of this nonlinear transformation, an analytic expression of the physical parameter distribution and its confidence intervals is tedious to obtain. However, with (53), together with the assumption of an unbiased GTLS estimator, a multivariate normal distribution N ( ˆGTLS (k ), ) for the discrete system parameters is defined. From this distribution, a stochastic sample can be drawn which is then inserted into the transformation (55). The hypothesis, that the distribution of the physical parameter estimates is normal, has been rejected. Therefore the confidence bounds cannot be determined from the sample covariance and Student's t-statistic. Instead, the 95% confidence intervals are calculated numerically from the physical parameter samples, which is known as bootstrapping 159

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

Fig. 8. Time resolved physical parameters from transient experimental data. The most prominent effect due to the direct current jump is visible in the resistance R2 , corresponding to the RC element describing the low frequency impedance arc. Vertical dashed lines mark the time interval during which the presence of current jump in the sliding data window distorts the estimation. The gray area around the parameter estimates denotes the 95% confidence interval.

Acknowledgment This work was supported by the project SoH4PEM (FFG. No. 854867) in cooperation with AVL List GmbH. Moreover, the authors would like to thank Abram Krebs (Gamry Instruments, Inc.) for providing scripts and technical support to enable APRBS measurements. Furthermore the authors appreciate the support of Hubert Gasteiger for his scientific advice and the use of his laboratory.

[4] [5] [6]

References

[7]

[1] E.L. Miller, D. Papageorgopoulos, N. Stetson, K. Randolph, D. Peterson, K. CierpikGold, A. Wilson, V. Trejos, J.C. Gomez, N. Rustagi, et al., US department of energy hydrogen and fuel cells program: progress, challenges and future directions, MRS Advances 1 (42) (2016) 2839–2855. [2] V. Das, S. Padmanaban, K. Venkitusamy, R. Selvamuthukumaran, F. Blaabjerg, P. Siano, Recent advances and challenges of fuel cell based power system architectures and control–A review, Renew. Sustain. Energy Rev. 73 (2017) 10–18. [3] L. Dubau, L. Castanheira, F. Maillard, M. Chatenet, O. Lottin, G. Maranzana,

[8] [9] [10]

160

J. Dillet, A. Lamibrac, J.-C. Perrin, E. Moukheiber, et al., A review of PEM fuel cell durability: materials degradation, local heterogeneities of aging and possible mitigation strategies, Wiley Interdisciplinary Reviews: Energy Environ. 3 (6) (2014) 540–560. M. Jouin, R. Gouriveau, D. Hissel, M.-C. Péra, N. Zerhouni, Prognostics and health management of PEMFC–state of the art and remaining challenges, Int. J. Hydrogen Energy 38 (35) (2013) 15307–15317. T. Sutharssan, D. Montalvao, Y.K. Chen, W.-C. Wang, C. Pisac, H. Elemara, A review on prognostics and health monitoring of proton exchange membrane fuel cell, Renew. Sustain. Energy Rev. 75 (2017) 440–450. J. Wu, X.Z. Yuan, H. Wang, M. Blanco, J.J. Martin, J. Zhang, Diagnostic tools in PEM fuel cell research: Part I Electrochemical techniques, Int. J. Hydrogen Energy 33 (6) (2008) 1735–1746. J. R. Macdonald, E. Barsoukov, Impedance spectroscopy: theory, experiment, and applications, History 1 (8). N. Fouquet, C. Doulet, C. Nouillant, G. Dauphin-Tanguy, B. Ould-Bouamama, Model based PEM fuel cell state-of-health monitoring via ac impedance measurements, J. Power Sources 159 (2) (2006) 905–913. C. Brunetto, A. Moschetto, G. Tina, PEM fuel cell testing by electrochemical impedance spectroscopy, Elec. Power Syst. Res. 79 (1) (2009) 17–26. W. Choi, J.W. Howze, P. Enjeti, Development of an equivalent circuit model of a fuel cell to evaluate the effects of inverter ripple current, J. Power Sources 158 (2) (2006) 1324–1332.

Journal of Power Sources 405 (2018) 150–161

D. Ritzberger et al.

identification of discrete-time dynamic system, Int. J. Eng. Technol. [36] R. Ionutiu, J. Rommes, A framework for synthesis of reduced order models, (CASAreport; Vol. 0928) Eindhoven: Technische Universiteit Eindhoven. [37] A. Lasia, Electrochemical Impedance Spectroscopy and its Applications vol. 7, Springer, 2014. [38] M.U. Iftikhar, D. Riu, F. Druart, S. Rosini, Y. Bultel, N. Retière, Dynamic modeling of proton exchange membrane fuel cell using non-integer derivatives, J. Power Sources 160 (2) (2006) 1170–1182. [39] S. Alavi, C. Birkl, D. Howey, Time-domain fitting of battery electrochemical impedance models, J. Power Sources 288 (2015) 345–352. [40] C. Eckart, G. Young, The approximation of one matrix by another of lower rank, Psychometrika 1 (3) (1936) 211–218. [41] C. Hametner, S. Jakubek, Nonlinear identification with local model networks using GTLS techniques and equality constraints, IEEE Transactions on Neural Networks 22 (9) (2011) 1406–1418. [42] J. Rommes, Arnoldi and Jacobi-Davidson methods for generalized eigenvalue problems with singular, Math. Comput. 77 (262) (2008) 995–1015. [43] W.W. Hager, Updating the inverse of a matrix, SIAM Rev. 31 (2) (1989) 221–239. [44] J. Jiang, Y. Zhang, A novel variable-length sliding window blockwise least-squares algorithm for on-line estimation of time-varying parameters, Int. J. Adapt. Contr. Signal Process. 18 (6) (2004) 505–521. [45] M. Saccoccio, T.H. Wan, C. Chen, F. Ciucci, Optimal regularization in distribution of relaxation times applied to electrochemical impedance spectroscopy: ridge and lasso regression methods-a theoretical and experimental study, Electrochim. Acta 147 (2014) 470–482. [46] A. Weiß, S. Schindler, S. Galbiati, M.A. Danzer, R. Zeis, Distribution of relaxation times analysis of high-temperature pem fuel cell impedance spectra, Electrochim. Acta 230 (2017) 391–398. [47] A. Oppenheim, R. Schafer, R. John, Discrete-time Signal Processing, Prentice Hall Inc, New Jersey, 1989. [48] E. A. Guillemin, Synthesis of Passive Networks: Theory and Methods Appropriate to the Realization and Approximation Problems. [49] C. de Beer, P.S. Barendse, P. Pillay, B. Bullecks, R. Rengaswamy, Online diagnostics of HTPEM fuel cells using small amplitude transient analysis for CO poisoning, IEEE Trans. Ind. Electron. 62 (8) (2015) 5175–5186. [50] Y. Tang, W. Yuan, M. Pan, Z. Li, G. Chen, Y. Li, Experimental investigation of dynamic performance and transient responses of a kw-class pem fuel cell stack under various load changes, Appl. Energy 87 (4) (2010) 1410–1417. [51] C. Raga, A. Barrado, A. Lázaro, C. Fernández, V. Valdivia, I. Quesada, L. Gauchía, Black-box model, identification technique and frequency analysis for pem fuel cell with overshooted transient response, IEEE Trans. Power Electron. 29 (10) (2014) 5334–5346. [52] C. Simon, F. Hasché, D. Müller, H.A. Gasteiger, Influence of the gas diffusion layer compression on the oxygen mass transport in PEM fuel cells, ECS Transactions 69 (17) (2015) 1293–1302. [53] C. Simon, F. Hasché, H.A. Gasteiger, Influence of the gas diffusion layer compression on the oxygen transport in PEM fuel cells at high water saturation levels, J. Electrochem. Soc. 164 (6) (2017) F591–F599. [54] M. Bensimhoun, N-dimensional cumulative function, and other useful facts about Gaussians and normal densities, Jerusalem, Israel, Tech. Rep. (2009) 1–8. [55] S. Keller, T. Özel, A.-C. Scherzer, D. Gerteisen, U. Groos, C. Hebling, Y. Manoli, Characteristic time constants derived from the low-frequency arc of impedance spectra of fuel cell stacks, Journal of Electrochemical Energy Conversion and Storage 15 (2) (2018) 021002. [56] C. Zhang, W. Allafi, Q. Dinh, P. Ascencio, J. Marco, Online estimation of battery equivalent circuit model parameters and state of charge using decoupled least squares technique, Energy 142 (2018) 678–688. [57] M. Eikerling, A. Kornyshev, Electrochemical impedance of the cathode catalyst layer in polymer electrolyte fuel cells, J. Electroanal. Chem. 475 (2) (1999) 107–123. [58] B. Efron, R.J. Tibshirani, An Introduction to the Bootstrap, CRC press, 1994. [59] A. Honkela, Approximating nonlinear transformations of probability distributions for nonlinear independent component analysis, Neural Networks, 2004. Proceedings. 2004 IEEE International Joint Conference on, vol. 3, IEEE, 2004, pp. 2169–2174.

[11] E. Lechartier, E. Laffly, M.-C. Péra, R. Gouriveau, D. Hissel, N. Zerhouni, Proton exchange membrane fuel cell behavioral model suitable for prognostics, Int. J. Hydrogen Energy 40 (26) (2015) 8384–8397. [12] S. Asghari, A. Mokmeli, M. Samavati, Study of PEM fuel cell performance by electrochemical impedance spectroscopy, Int. J. Hydrogen Energy 35 (17) (2010) 9283–9290. [13] W.J. Wruck, R.M. Machado, T.W. Chapman, Current interruptioninstrumentation and applications, J. Electrochem. Soc. 134 (3) (1987) 539–546. [14] A. Forrai, H. Funato, Y. Yanagita, Y. Kato, Fuel-cell parameter estimation and diagnostics, IEEE Transactions on Energy Conversion 20 (3) (2005) 668–675. [15] M. Rubio, A. Urquia, S. Dormido, Diagnosis of PEM fuel cells through current interruption, J. Power Sources 171 (2) (2007) 670–677. [16] L. Ljung, System Identification: Theory for the User, second ed., Prentice Hall PTR, 2006. [17] A. Debenjak, P. Boskoski, B. Musizza, J. Petrovcic, D. Juricic, Fast measurement of proton exchange membrane fuel cell impedance based on pseudo-random binary sequence perturbation signals and continuous wavelet transform, J. Power Sources 254 (2014) 112–118. [18] Z. Zheng, R. Petrone, M.-C. Péra, D. Hissel, M. Becherif, C. Pianese, N.Y. Steiner, M. Sorrentino, A review on non-model based diagnosis methodologies for PEM fuel cell stacks and systems, Int. J. Hydrogen Energy 38 (21) (2013) 8914–8926. [19] O. Nelles, Nonlinear System Identification: from Classical Approaches to Neural Networks and Fuzzy Models, Springer Science & Business Media, 2013. [20] R. Isermann, Process fault detection based on modeling and estimation methods: a survey, Automatica 20 (4) (1984) 387–404. [21] M.A. Danzer, E.P. Hofer, Electrochemical parameter identificationan efficient method for fuel cell impedance characterisation, J. Power Sources 183 (1) (2008) 55–61. [22] M. Rubio, A. Urquia, R. Kuhn, S. Dormido, Electrochemical parameter estimation in operating proton exchange membrane fuel cells, J. Power Sources 183 (1) (2008) 118–125. [23] C. Jeppesen, S.S. Araya, S.L. Sahlin, S.J. Andreasen, S.K. Kær, An eis alternative for impedance measurement of a high temperature pem fuel cell stack based on current pulse injection, Int. J. Hydrogen Energy 42 (24) (2017) 15851–15860. [24] N. Fouquet, Real time model-based monitoring of a pem fuel cell flooding and drying out, Vehicle Power and Propulsion Conference (VPPC), IEEE, 2010, pp. 1–8 IEEE, 2010. [25] C. Fleischer, W. Waag, H.-M. Heyn, D.U. Sauer, On-line adaptive battery impedance parameter and state estimation considering physical principles in reduced order equivalent circuit battery models: Part 1. Requirements, critical review of methods and modeling, J. Power Sources 260 (2014) 276–291. [26] C. Fleischer, W. Waag, H.-M. Heyn, D.U. Sauer, On-line adaptive battery impedance parameter and state estimation considering physical principles in reduced order equivalent circuit battery models part 2. Parameter and state estimation, J. Power Sources 262 (2014) 457–482. [27] B. Xia, X. Zhao, R. De Callafon, H. Garnier, T. Nguyen, C. Mi, Accurate Lithium-ion battery parameter estimation with continuous-time system identification methods, Appl. Energy 179 (2016) 426–436. [28] B. Fridholm, T. Wik, M. Nilsson, Robust recursive impedance estimation for automotive lithium-ion batteries, J. Power Sources 304 (2016) 33–41. [29] P. Stoica, T. Söderström, Bias correction in least-squares identification, Int. J. Contr. 35 (3) (1982) 449–457. [30] S.V. Huffel, P. Lemmerling, Total Least Squares and Errors-in-variables Modeling: Analysis, Algorithms and Applications, Springer Science & Business Media, 2013. [31] C.C. Paige, M. Wei, Analysis of the generalized total least squares problemAX B when some columns ofA are free of error, Numer. Math. 65 (1) (1993) 177–202. [32] S.V. Huffel, J. Vandewalle, Analysis and properties of the generalized total least squares problem AXB when some or all columns in A are subject to error, SIAM J. Matrix Anal. Appl. 10 (3) (1989) 294–315. [33] S. Jakubek, C. Hametner, Identification of neurofuzzy models using GTLS parameter estimation, IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 39 (5) (2009) 1121–1133. [34] H. Akaike, Canonical correlation analysis of time series and the use of an information criterion, Mathematics in Science and Engineering, vol. 126, Elsevier, 1976, pp. 27–96. [35] M. F. A. Samad, A. R. M. Nasir, Comparison of information criterion on

161