Iterative Fuzzy Modeling of Hydrogen Fuel Cells by the Extended Kalman Filter

Hydrogen economy is one of the recently opened alternatives in the field of non-polluting energy. Hydrogen fuel cells show high performance, high reliability in stationary applications and minimal environmental impact. To increase the efficiency of the hydrogen fuel cell it is very important to have a good model to predict its dynamic behavior. In addition, this model must be able to adapt iteratively to the changes that occur in its performance due to operating conditions and even to the degradation through its lifespan. This paper presents the application of an iterative fuzzy modeling methodology based on the extended Kalman filter applied to a real hydrogen fuel cell. Two algorithms based on the Kalman filter will be compared with the well-known backpropagation algorithm from three different initializations: by uniform partitioning, subtractive clustering and CMeans clustering. The used data have been collected during the actual operation of a real 3.4 kW proton exchange membrane fuel cell. As the article experimentally shows, the Takagi-Sugeno type fuzzy model allows to create a very accurate nonlinear dynamic model of the fuel cell, which can be very useful to design an efficient fuel cell control system.


I. INTRODUCTION
Today's society faces the challenge of changing the economic model to one that is more sustainable and respectful with the environment. In this way, the main aim of new energy conversion technologies is the production of non-polluting energy. The so-called ''hydrogen economy'' is one of the recently opened alternatives in the field of non-polluting energy. Hydrogen can be ecologically produced from renewable sources using electrolyzers. Also it can be stored (under pressure, in the form of metal hydrides or even liquid, despite the expenditure of energy required for it) and can be converted back, in a non-polluting way, into electric power by hydrogen fuel cells (hereinafter FC both plural and singular). These three characteristics make it possible to create a circular energy cycle based on clean or green hydrogen (productionstorage-usage) [1]. FC are nowadays, despite being a technology in continuous evolution, a realistic (and practical) The associate editor coordinating the review of this manuscript and approving it for publication was Xiaochun Cheng. solution thanks to its technological development (mainly due to the automotive field), high performance, high reliability and minimum environmental impact [2].
A FC is a nonlinear complex system based on the serial connection of individual cells making up a stack, together with other systems necessary to its proper operation (the so-called balance of plant or BoP, which is governed by a control system responsible for the best performance of the FC [3]). Compared to other green technologies, for example wind or photovoltaic generation [4], FC can operate in any geographic location, without taking into account environmental constraints (availability of sun, wind, etc.). Specifically, proton exchange membrane fuel cells, also known as polymer electrolyte membrane fuel cells (PEMFC), are considered one of the best FC alternatives for transport applications as well as stationary and portable applications. Other important advantages of PEMFC are its low operating temperature (50 o C-100 o C) and the fact that the only by-product resulting its operation is water, which allows a circular cycle if it is recirculated back to the electrolyzer to produce hydrogen again. Regarding electrical grids applications, PEMFC can be connected to them [5], or installed as grid-independent generators [6]. High power FC could be used connected to the grid in power stations [5]- [9]; while smaller FC could be used in mobile power stations [10]. Additionally, as an important difference with regard to the usual batteries, a FC can supply energy continuously, whenever there is hydrogen available at its input. In transport applications, there is a lot of successful research on FC electric vehicles (scooter, car, bus or truck), with realities already present in the market [2], [8].
The output voltage of a PEMFC can vary depending on the state of its stack: temperature, hydrogen pressure, etc. [11], [12], and especially with the demand for current at its output [13], [14]. To increase applicability and efficiency of the FC it is very important to have a model to predict its dynamic behavior [15]- [18], which is particularly important in an FC due to its relatively slow response to rapid changes in demand for current at its output. However, it is not easy to find an accurate mathematical model due, mainly, to the large number of physicochemical parameters and laws that need to be taken into account, as well as to the nonlinear intrinsic behavior of the FC. In addition, a FC is subject to considerable changes in its dynamic behavior throughout its lifespan; mainly due to the stack degradation. Thus, a model-based FC controller can lose efficiency over the time, whereby would be very interesting to have different FC models through its lifespan. However, it is not easy to carry out this by mathematical models, whereby it is increasing the interest in modeling techniques based on input-output data that reflects the current state of the FC [19]. This would allow to have different FC models over its lifespan, making easier the system analysis [20]- [24] and the controller setting from its initial design, when the FC was brand new [25]- [28].
There are numerous modeling algorithms present in the literature, both for offline [29]- [31] and iterative modeling [32], [33]; however, when these algorithms are going to be used in real systems, it is necessary to take into account considerations such as the convergence speed or its robustness against the measurement noise.
The aim of this paper is to present an iterative modeling of a real PEMFC, efficient and robust to the presence of measurement noise. For this, a Takagi-Sugeno (TS) fuzzy model [34] hybridized with the Kalman's filter [35] will be used; since both, working together, allow to obtain very simple, accurate, efficient and robust nonlinear models in the presence of noise. Fuzzy models are rule-based models, where the output of the model is computed as the combination of all the outputs of each of its rules based on the degree of compliance of each of them [36].
Fuzzy modeling, particularly modeling based on TS fuzzy systems, allows to obtain accurate models from a small number of rules [37] and, therefore, it can be a proper technique to model a nonlinear plant as a PEMFC, whose parameters also change over time and use. The fuzzy modeling methodology developed by the authors is based on the extended Kalman filter (EKF) [38], and it allows iterative modeling of nonlinear systems based on input-output data in presence of measurement noise. The developed methodology can work iteratively, online with the system, even in presence of noise, and is as computationally efficient as the Kalman filter is [39]. Modeling using input-output data allows to obtain, in each time, realistic models of the system under study; which, over the lifetime of the system, will be more accurate than theoretical ones because they reflect the current state of the system. The Kalman filter (KF) [35] is an efficient recursive filter that estimates the internal states of a dynamic system from a series of noisy measurements. It uses a series of observed over time data, that is capable to lead with noise and other inaccuracies, to estimate unknown variables more accurately than techniques based on a single measurement [40]. The BoP required to operate a FC is composed of switching devices (inherently noisy), which makes it very interesting that the modeling algorithm can work properly in the presence of noise. In addition, in many cases, not only for modeling an FC, it is required that modeling algorithm works iteratively to adapt the model to the changes in plant behavior. KF is used in a wide range of engineering applications and it is an important topic in control systems theory and control engineering. If a non-linear system, such as a FC, needs to be modeled, the EKF is an interesting option to be considered [41] when the system supports linearized models around any working point.
The KF hybridized with fuzzy logic has been previously used in several applications [42]- [44]. In 2002, Simon introduced the use of KF to adjust the parameters of a fuzzy model [45]. Later, other proposals were made [39], [46]- [48] to improve and generalize the first ones.
In short, the combined use of a TS fuzzy model and a Kalman filter allows obtaining a fast and efficient modeling algorithm, which is capable of working iteratively and in the presence of noise, something very useful in real systems. Based on the previous analysis, this article presents a general methodology of FC fuzzy modeling based on the EKF, which allows to obtain accurate models iteratively. This methodology works in line with the system, in the presence of noise, and it is computationally efficient. This methodology is applied to the online modeling of a real 3.4 kW PEMFC commercial stack.
After this general introduction to the problem that is intended to be addressed, and the principle of operation of PEMFC that will be used, this article is organized as follows: First, in section II, a review of the latest advances present in the literature, as well as highlighting the novelty of this work. In the section III, fuzzy TS modeling, the Kalman filter, and its application to fuzzy TS modeling are introduced. In this section, the modeling techniques used in this paper are presented, along with the notation used. Section IV is devoted to the application of these techniques to the modeling of a real PEMFC. This section shows the characteristics of the performed tests and their results. Subsequently, in section V, a discussion of the most relevant results from the authors'  point of view is carried out. Finally, some conclusions are presented in section VI.

II. RELATED WORKS
Regarding FC systems, as it can be deduced from [49] for a proper operation and accurate diagnosis for the detection and identification of faults, it is crucial to have a good model to predict the FC system behavior, not only from a static point of view but also dynamic throughout.
Despite the high number of models than can be found [50], they can be classified as shown in Table 2. Regarding this classification, FC system dynamic models based on empirical approach for the electrochemistry description (highlighted in black) are required to model important transients such as startup, shutdown and load changes. Additionally, they are also employed for fuel cell system degradation as the thermal stress associated with load and thermal cycling that may contribute to cell failure, and of course, dynamic models are prerequisites for control systems design. A control system automatically regulates the response of the system and keeps it at the desired value by manipulating some variables such as temperature, flow rate, or composition of reactant streams.
Exploring in the bib liography in the last decade, Kim et al. proposed in [51] a dynamic model that reflects the stack temperature evolution. Physical parameters for the proposed model were estimated based on the steady and transient results experimentally, obtained due to the stepwise variation parameters adjusted with an excel solver. Despite the model is accurate it does not include pressure behavior, so it is not possible to shows the reactant gas effects.
In this same research line, Qi et al. [52] proposed a model that includes the effects of stack temperature and reactant gas flow. In spite of the model is more detailed, the dynamic behavior is defined as the conventional ''charge double layer'' phenomenon by means the adjustment of theoretical parameters, and the proposed model does not include the BoP. A more complex model is proposed in [53], where authors address a PEFC dynamic model associated with the reactant flow, cooling waterflow, and water phase change. Nevertheless, the model is also based on theoretical parameters adjustment and assumes auxiliary components are well controlled, ignoring their impacts on fuel cell behavior.
Advanced modeling techniques are used by Na et al. in [54], and Chang et al. in [55], where authors propose a multiple-input and multiple-output (MIMO) nonlinear PEMFC model that can directly utilize the feedback linearization for the nonlinear control. Recently, adaptive neural network and fuzzy logic controllers were applied in [56], [57]. Benchouia et al. [56] develop the models of stack voltage and stack power, but keeps out the stack temperature and pressure reactant gas, while and Abbaspour et al. [57] include reactant gases pressure and flow rate but do not take into account the stack temperature.
Based on the above, Table 3 summarizes the main contributions of this paper and explains why it is novel. The submitted proposal develops a fuzzy model of a FC system based on data in the presence of measurement noise. Additionally, the developed methodology can work iteratively, online with the system, even in presence of noise.

III. MATERIAL AND METHODS
To collect the data of this research, a 3.4 kW FCgen R -1020ACS fuel cell stack integrated in a work bench forming a FC has been implemented [58] (Fig. 1 shows the used diagram and Fig. 2 shows the real implementation). This stack is formed by 80 BAM4G polymeric single cells [59] with a porous carbon cloth anode and cathode and a platinum based catalyst [60]. The whole stack has graphite plates between cells and aluminum end plates, all of them join by compression. The pressure of the hydrogen inlet is around 1.36 bar. The maximum output power of the stack is 3.4 kW, with 45.33 V and 75 A as typical values for voltage and current, respectively. The stack has its own air based refrigeration system and an oxidant subsystem; both built by the authors, following the manufacturer's instructions [58]. Fig. 1 shows all the subsystems that make up the BoP that, together with the commercial 3.4 kW stack, configure the FC workbench. The real laboratory equipment used  in this research is shown in Fig. 2 [61]. To perform different tests, a AMREL R PLA5K-120-1200 programmable electronic load has been used; with it, different output current profiles can be carried out. The SCADA (supervisory control and data acquisition system), described in [62] and [63], collects the measurements of the sensors (values of the variables of interest: voltage, current, temperature, pressure and hydrogen flow, etc.) and stores them.
Hereinafter in this section, a methodology already published by authors for online TS fuzzy modeling by the EKF will be briefly introduced.

A. TAKAGI-SUGENO FUZZY MODELS
For the formal implementation of control strategies on a PEMFC, it is very important to have an accurate and efficient model of the system [17], [64]. The dynamic behavior of a FC is nonlinear; this is why in this work it is proposed to use TS fuzzy models, since this type of models has demonstrated to be very accurate in any type of systems [65], [66].
Fuzzy logic is a multivalued logic that determine whether or not an element belongs to a set, but the quality of membership is not only ''yes'' or ''not'' (traditional bi-valued logic), but it is evaluated in a full range from non-membership to full membership. Fuzzy models are built based on rules of the IF-THEN form, where IF the premise of the rule is fulfilled, THEN the consequence of the rule is the desired value or the action to be carried out. There are different types of fuzzy models, among them, TS models, in which the antecedent of the rule is built with the model inputs in the form of fuzzy sets, and the consequent by a linear polynomial with all the model inputs plus an affine term (as can be seen in (1)). Although a fuzzy model has a linguistic form, its output can be expressed analytically, allowing linguistic interpretability to be combined with mathematical precision. This way of constructing the consequent allows TS models to be universal approximators, and achieve high accuracy with a small number of rules [67].
These characteristics of fuzzy logic and fuzzy models allow approximate reasoning, dealing innately with uncertainties, which gives it great skill for accurate and robust modeling [36], [68]. On the other hand, is relatively easy to convert them into nonlinear state models [65], which allows the formal analysis required in control engineering. However, it is also known that the number of rules in TS fuzzy models is increased as a lower approximation error is desired [69]. This implies that the modeling process is very important, both for analysis [70]- [72] and the design of fuzzy control systems [73], [74].
Let n be the number of input variables, x j , and m the number of output variables, y i , of a completely general system to model; its discrete TS fuzzy model can be represented by the following set of rules [34], [75]:   If the input vector x = (x 1 , x 2 , . . . , x n ) T extends in a coordinate [65], [73] byx 0 = 1, the extended vectorx takes the form:x Using the previous expression, the output y i can be calculated by [37]: being a ji (x) variables coefficients [76] defined by where w l i (x) represents the degree of activation of the fuzzy model rules: where µ l ji (x j (k), σ l ji ) represents the j-th membership function of the l rule for the i-th model output, which determines the fuzzy set A l ji . The σ l ji elements represent the set of adaptive parameters of the j-th membership functions; so these values along with a l ji must be calculated using the chosen estimation algorithm to obtain a suitable TS model.

B. EXTENDED KALMAN FILTER
KF allows to construct an optimal observer for linear systems in presence of zero mean Gaussian white noise, both in model and measurements [35], [38]. To cover nonlinear systems the KF has been adapted through the EKF [41]. EKF requires that the system supports linearized models around any working point. Although the EKF is not optimal, since it is based on a linear approximation, it is a powerful tool for estimation in noisy environments.
Considering a nonlinear discrete system: where x(k) is the state vector, u(k) the input vector, and v(k) and e(k) are white noise vectors. The Jacobian matrices of the system are: and Although both the KF and EKF can be written as a single equation, both can be conceptually considered as the sum of two distinct phases: predict and update. The predict phase uses the estimated state from the previous time-step, (k − 1), to produce an estimation of the state at the current time-step, k, a priori estimation. In the update phase, this prediction is combined with current observation information to refine the estimate state, a posteriori estimation. The EKF can be solved by iterative application of the following set of equations: Predict: Update: x(k|k) =x(k|k − 1) + K(k) y(k) −ŷ(k) (13) wherex(k|k −1) andx(k|k) represent, respectively, the priori and posteriori estimations of the true state x. P(k|k − 1) is the a priori estimate covariance matrix and P(k|k) is the a posteriori one. K(k) is the Kalman gain,ŷ(k) are the estimated outputs and, R v , R ve and R e , are the noise covariance matrices.
The iterative process starts with an initial estimate of state vector,x(0), and an initial value of the covariance matrix, P(0). From this initial estimate, the algorithm run iteratively with the system, obtaining a solution that minimizes both, estimation error and its covariance matrix for the linearization obtained at each time-step.

C. APPLICATION OF THE EKF TO TS FUZZY MODELING
The EKF can be applied to estimate the parameters of TS fuzzy systems. Firstly, it is necessary to model the system in accordance with the framework of the KF, i.e., to build a system whose states depend directly on the parameters to be estimated [45]. In (15) p(k) is the vector of parameters to be estimated and y(k) the output vector of the TS fuzzy system; e(k) is a zero mean Gaussian white noise that represents the uncertainty of the system output measurements (sensors accuracy), whose covariance is R e . Then, the expressions from (10) to (14) can be recursively applied to obtain an estimation of the true vector p(k). p(k + 1) = p(k) y(k) = h(x(k), p(k)) + e(k).
Applying (7), (8) and (9) on (15): and beingp(k) the current estimation of the parameters vector of the TS fuzzy model. Sop(k) represents both antecedents, σ l ji , and consequents, a l ji , parameters of the TS fuzzy model. Therefore, h(x(k), p(k)) corresponds to (3) which is nonlinear; regarding (18), C(p(k)) must be obtained from the derivative of h(·) with respect to each of adaptive parameters of the TS fuzzy model. Operating, it is possible to obtain that (19) represents the derivative of the TS fuzzy model with respect to the consequents, and (20) respect to antecedents [48]: being L, J and I determine each particular parameter of the set of the consequent and antecedent parameter vectors. The algorithms to obtain these derivatives are implemented in the open source library Fuzzy Logic Tools (FLT) [77]. Depending on the adjustment needs of the model, it is possible to use only one filter to determine the consequent ones, or sequentially execute two filters to adjust consequents and antecedents [48]. This article will test the response of both options to verify both, their accuracy in modeling, and the robustness of these two algorithms against the modification of their parameters. The algorithm that modifies only the consequent ones will be called Kalman c , and Kalman c+a that modifies both the antecedents and consequents of the fuzzy model. These algorithms are shown in Fig. 3, where the adjustment of the antecedents being applicable only for the Kalman c+a algorithm. About the initialization of the algorithms, if there is no information about the system it is possible to initialize all consequents to 0, or use other initialization procedures [78]. Antecedents can be initialized using an uniform distribution VOLUME 8, 2020 or by any clustering algorithm [79]- [81]. The antecedents and consequents covariance matrices can be initialized with a diagonal matrix multiplied by a positive integer that indicates the degree of confidence in these initial parameters. These initialization parameters will be called α for the consequents and β for the antecedents. In this way, the Kalman filter matrix P(0) is created as αI for the consequents, and as βI for the antecedents, being I an identity matrix of the appropriate size. A small value of α or β indicates to the Kalman filter that the initial parameters may be close to their real value, while a high value would indicate the opposite. From a practical point of view, a small value will make the filter behave more conservatively at startup (it may be slower but there is less risk of non-convergence), and a higher value will make the adjustment more aggressive (it may increase its speed, but also increasing the risk of non-convergence of the algorithm).

IV. INTERACTIVE FUZZY MODELING OF A REAL HYDROGEN FUEL CELL
In this section, the modeling methodology previously described in section III-C will be used to obtain the interest TS fuzzy models (as the described in the section III-A) of a real 3.4 kW PEMFC (Figures 1 and 2) [58]. These TS fuzzy models, adjusting their consequents (subscript 'c' in tables and figures from now on) or their antecedents and consequents (subscript 'c+a' in tables and figures from now on), in both cases by the EFK, will be compared to the adjusted by the well known backpropagation algorithm [82]. The tests will be carried out on two sets of data obtained from the PEMFC. Three different initializations methods will be used: uniform partitioning [83], subtractive clustering [84] and CMeans clustering [80]. These three forms of initialization have been chosen because they are widely used in the literature. Uniform partitioning allows you to make an initial model without knowledge of the system, simply by dividing the parameter space into equal spaces. Subtractive clustering is one of the most widely used algorithms for obtaining initial TS fuzzy models, since it allows generating a good initial model based on a reduced set of rules. Finally, the CMeans algorithm is a more traditional alternative that has been used to check the behavior of the modeling against different starting point. This will allow checking the behavior of the models based on different algorithms, with different degrees of initial error and different antecedents distributions. Although the developed algorithm supports the use of different types of membership functions, even mixed together, Gaussian type membership functions have been used in all cases since these were the ones that showed the best results in preliminary tests. The distribution of these membership functions will be obtained with the different clustering algorithms used.
The data shown in Fig. 4 will be used for modeling, and those shown in Fig. 5 for validation (following Fig. 1: output current is I S , output voltage is V S , stack temperature is T S , and hydrogen inlet pressure is P H 2 ). The time-step (sampling time) in both data sets is 0.5 s. The output of the  models will be the prediction in the next sampling time of the PEMFC output voltage. The output current, the hydrogen inlet pressure and the stack temperature, as well as delays of these signals and the output itself, will be used as input signals. The dynamics of a system can be represented by using time-delayed signals, including output delays. Since these signals are already known from previous iterations, they can be used to improve model predictions. The number of delays used as inputs to the dynamic model is a parameter that must also be adjusted to obtain a good model. A small number will not correctly capture the order of the system; while an excessive number is even worse since it will increase the number of parameters, complicating the model, and will not allow to capture the dynamics of the system adequately when using a higher order than the real system. The number of delays will be a parameter to change in the tests in order to find the most optimal values, and also to check the robustness of the algorithms to their modification. Multiple tests will also be carried out by changing the parameters setting of the different algorithms in order to obtain the best possible performance, as well as analyze their robustness against the parameter change. The fact of using a real FC to perform the tests makes the measurements have some noise. The connections of the sensors wires to the control unit can be seen in Fig. 2, and the existence of noise in the signals in Fig. 4 and 5, particularly in voltage and pressure measurements. One of the advantages of using the KF should be its robustness against noise, so it will be checked if this is the case in each test. To obtain the covariance matrices of the noise measurements, R e (see (15)) for each input, each signal has been filtered by means of a mean filter [85] with 10 samples. Then, the covariance of the resulting noise has been calculated.
The Mean Absolute Error (MAE) will be used to compare the results of the different models, as it gives an interpretable idea of their actual error: being N the number of samples.

A. UNIFORM PARTITIONING INITIAL MODEL
First, a completely unknown model is assumed to which an iterative adjustment will be made. No prior learning will be done to this model. This test will allow to determine the convergence speed of the algorithms in the iterative adjustment starting from a model very far from the real system. The initial model was created by a uniform partitioning of the discourse universes of the input variables [83], using 2 antecedents for each input. Note that a rule will be created for each possible combination of the antecedents of the inputs, so if a larger number is used the model would be excessively complex, i.e., with too many rules and parameters. On the other hand, it was assumed that the initial consequents were all zero. This methodology can be interesting when there are no initial data to create a model using clustering techniques. It has the great disadvantage that it generates many rules and, therefore, many parameters that must be adjusted. However, this excess of parameters can be useful to check the execution speed of the tested algorithms. The results of the tests performed for different parameters can be seen in Table 4 (note that α and β are the coefficients to initialize the consequents and antecedents covariance matrices respectively). Table 4 has been ordered according to the MAE of the Kalman c algorithm, which is the one with the least error. The models used to obtain Figures 6 to 9 have been the best models obtained by each of the algorithms evaluated using the validation data, that is, row 1 of Table 4 for backpropagation and Kalman c algorithms, and row 4 for Kalman c+a algorithm. Fig. 6 shows the output voltage, V S , of the real system and the best models obtained with each algorithm, Fig. 7 shows the errors both during the modeling and validation phases, and Figures 8 and 9 show the initial and final time of the modeling phase, respectively.

B. SUBTRACTIVE CLUSTERING INITIAL MODEL
The next test use modeling data shown in Fig. 4 to create an initial model based on the Chiu subtractive clustering VOLUME 8, 2020   algorithm [84]. The initial model will be adjusted using each of the studied algorithms (backpropagation, Kalman adjusting only the consequents (Kalman c ), and Kalman adjusting antecedents and consequents (Kalman c+a ) using the same modeling data. Finally, as in the previous test, the final models will be checked with the validation data shown in Fig. 5 without making adjustments on them. As in the previous case, the α and β parameters of the algorithms, as well as the input and output delays, will vary between tests to obtain their most appropriate values, (see Table 5). In this case, the Cluster center's range of influence parameter (RADII) of the clustering algorithm will also be changed between tests. From a practical point of view, this parameter inversely influences the number of rules that will be created (a higher value implies that the algorithm will create fewer rules). The best results for the backpropagation and Kalman c algorithms are shown in row 4, while the best result of the Kalman c+a algorithm is shown in row 1. Table 5 has been ordered according to the MAE of the Kalman c+a algorithm, which is the one with the least error.
The most interesting results of these tests are shown in Table 5 (of the 206 tests performed, only the most relevant results are shown in order to present a smaller and clearer table). The models used to obtain Figures 11 to 13 have been the best models obtained by each of the algorithms evaluated using the validation data, that is, row 1 of Table 5 for Kalman c+a algorithm, and row 4 for backpropagation and Kalman c algorithms. Fig. 11 shows the errors made both during the modeling and validation phases. Figures 12  and 13 show the initial and final time of the modeling phase, respectively.

C. CMeans CLUSTERING INITIAL MODEL
Finally, the last test use modeling data to create an initial model based on the CMeans clustering algorithm [80]. This test follows the same idea as the previous test, but using another clustering method to verify its effect on the modeling algorithms. Since the qualitative results are similar to the previous ones, only Table 6 is shown with the most interesting results of the performed tests. Table 6 has been ordered according to the MAE of the Kalman c+a algorithm, which is the one with the least error.

V. DISCUSSION
Next, main results obtained from the three initial model (uniform partitioning, subtractive clustering and cmeans clustering) are discussed in detail. Figure 6 shows that all models adequately predict the PEMFC output voltage, although those adjusted by the EFK show the best performance (in this case the Kalman c algorithm is slightly better). There are some modeling errors, shown more clearly in Fig. 7, which would be smaller with the validation data if the iterative adjustment had been made (obviously not performed during validation). Fig. 8 shows that the initial error is very high in all models, something logical due to the initial models are completely unknown. Kalman c and Kalman c+a algorithms quickly converge to small errors, while the backpropagation algorithm takes considerably longer. Fig. 9 shows the end of the test, where the models are fitted with the latest available data. Thus, the models are already very close to the final models because the latest test data is being processed; therefore at the end all respond in a similar way. It can be seen that, as expected, the backpropagation algorithm is more sensitive to noise than the algorithms based on the KF. The error peak is due to the sudden big change of the source current, I S , where it affect in a similar way to the 3 algorithms. Table 4 shows that the best results are obtained by the algorithm that adjust only the consequents using the KF (Kalman c ). Since uniform modeling covers the entire universe of discourse with antecedents of all possible options, a change in these antecedents proves a worsening of the model. This can also be verified by seeing that a lower β value improves the results. α and β are not parameters of the backpropagation algorithm, and β is not of the Kalman c algorithm, so there are lines for these algorithms where the results are repeated.

A. UNIFORM PARTITIONING INITIAL MODEL
Observing the obtained results it is verified that, for this initialization method, it is better to use only a delay in the output signal, V S , to use as input. If 2 delays are used in V S as inputs to the model, the results are worse (this result can be verified in lines 4 and 11 of Table 4). It can also be easily ruled out that the PEMFC is not a static system by observing the results of the last test. When a static model is assumed, without delays, the error obtained is the worst with the 3 algorithms as can be seen in the last row of Table 4. It is clear that the PEMFC is a dynamic system and cannot be properly modeled with a static model, therefore, this case will not be considered again in the other experiments. Keep in mind that this model is not worse for having fewer rules, the number of rules is less because fewer inputs are used. Since there are only 3 inputs Seeing how the change of each parameter affects the result of the algorithms, it can be checked that the Kalman c+a algorithm is more sensitive to the variations of β than those of α. On the other hand, the Kalman c algorithm seems quite robust to the variations of α parameter.
The execution times of each algorithm have not been included because they can depend heavily on the implementation of each of them, as well as the equipment used in the tests. However, to add some data in this regard, for the first test the average execution time of each algorithm has been: 54.74 ms for the backpropagation algorithm implemented by MATLAB R , 0.48 ms for Kalman c , and 2.23 ms for Kalman c+a , both implemented by the FLT library [77] in MATLAB R . VOLUME 8, 2020

B. SUBTRACTIVE CLUSTERING INITIAL MODEL
In this case the models adapt better to the PEMFC system since the initial model is more suitable than in the previous case. This behavior is logical because the clustering process performs a first modeling. Fig. 10 shows the outputs of the models and the real PEMFC, and Fig. 11 shows the modeling and validation errors. In the modeling initialization is perfectly observed that all the algorithms start from a small error due to the previous clustering process, as can be seen in Fig. 12. Note that this figure shows that the Kalman c+a algorithm begins an instant after the others. This is because the best model obtained with this algorithm uses two delays in the output (see V S delay in row 1 of Table 5), so it must wait another time-step before starting to iterate.
Again it can be seen that the backpropagation algorithm is affected by the noise of the measurements to a greater extent than the algorithms based on the KF, as Fig. 13 shows clearly.
Observing the most interesting results of shown in Table 5, it can be verified that, as in the previous case, adding delayed inputs does not improve the performance of the algorithms when modeling the PEMFC. Having an initial model with a better positioned antecedents, the Kalman c+a algorithm that can improve the results adjusting the antecedents, obtaining the best model of this test. Increasing the delays used in the output to 2 does improves the Kalman c+a results, but not the other two, which get worse again. However, the improvement is not very significant, 0.286 V versus 0.442 V of the best model with only 1 delay with a signal amplitude about 80 V. So this could be due to the increase in the number of rules due to the additional input, rather than the fact of increasing the delays applied as inputs.
Regarding the robustness against the change of the parameters, all the algorithms are sensitive to the change of the RADII parameter, but this is logical since this parameter influences in the number of initial rules that the clustering  algorithm generates. Although the Kalman c+a algorithm has obtained the best result, it is again very sensitive to the modification of the α and β parameters, especially this second one. This fact can be seen, for example, in rows 4 to 17 of Table 5. The Kalman c algorithm is more consistent in its results compared to the change of the α parameter, in addition to not being worse than the previous one. So, it would be the most recommended in general terms.

C. CMeans CLUSTERING INITIAL MODEL
Looking at the test results Table 6 with the CMeans algorithm, the best performance of the algorithms based on the KF can be verified. The parameters that get the best results from the evaluated algorithms are shown in row 1 for the Kalman c+a algorithm, in row 16 for the backpropagation algorithm, and in row 24 for the Kalman c algorithm. Although the best result is obtained in this test using a delay in I S , that is, using I S (k) and I S (k − 1) as inputs. Eliminating this delay does not make the result worse (0.345 V vs 0.365 V with other parameters, or 0.409 V with the same configuration for a signal amplitude about 80 V). Adding delays to other inputs only makes the results worse. Regarding the output, it is again confirmed that with a single delay it is sufficient to correctly model this PEMFC system.

VI. CONCLUSION
In this research, the application of an iterative fuzzy modeling methodology based on the extended Kalman filter has been applied to a real medium-power hydrogen fuel cell. Of course, this method is applicable not only to PEM FC technology, but any other, adjusting the algorithm parameters appropriately. Two algorithms based on the Kalman filter have been compared with the well-known backpropagation algorithm with three different initializations, by uniform partitioning, subtractive clustering and CMeans clustering. Results have shown that the algorithms based on the Kalman filter have a faster convergence than the backpropagation algorithm, and are more robust to the intrinsic noises that the real signals present. The algorithm based on the Kalman filter that uniquely adjusts the consequent ones, Kalman c , is more robust than the one that adjusts both antecedents and consequents ones, Kalman c+a , especially if the antecedents have not been properly created through clustering. However, the Kalman c+a algorithm well adjusted has the ability to generate better models.
To make Hydrogen economy a reality is crucial that fuel cells are reliable and have good performance. An accurate model of the fuel cell that predicts the fuel cell dynamic behavior, allows increase its efficiency. The fuzzy modeling methodology proposed in this paper is able to be adapted iteratively to the changes that occur in fuel cell performance due to operating conditions and even to the degradation through its lifespan. She has been working at the Department of Electronic Engineering, Computer Systems, and Automatic Control, University of Huelva, since 2003. She is coauthor of more than 20 articles published in JCR journals, one book chapter, and one book. She has participated in more than 30 research projects regional, national, and internationals.
Prof He is currently a Full Professor of systems engineering and automatic control with the University of Huelva. He has more than 300 publications, among them more than 100 articles published in indexed journals in the ISI Journal Citation Reports. Particularly, he has 55 publications included in the first quartile in 24 different journals; most of these are among the top ten in their categories, and several are number one. He holds 15 international patents. His main research interests include control engineering, renewable energy systems, and engineering education. Throughout his professional life, he has received 28 awards and academic honors. He has supervised 12 Doctoral Theses with eight awards. His H-index (SCOPUS) is 29, and the cites of his publications are near to 5000 (Google Scholar). He has led or co-led more than 60 research projects funded by public institutions and companies.