Proton exchange membrane fuel cell life prediction method based on comprehensive dynamic aging factor

Abstract

The invention discloses a method for predicting the service life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor. The method specifically comprises the following steps: step 1, drawing electrochemical impedance spectrograms at different moments and polarization curves at different moments; 2, constructing a proton exchange membrane fuel cell aging equivalent circuit model, and extracting a dynamic factor S1 data set based on the equivalent circuit model; 3, establishing a semi-empirical equation, associating the semi-empirical equation with the equivalent circuit model to obtain an association model, and extracting a dynamic factor S2 data set based on the association model; 4, performing normalization processing on the dynamic factors S1 and S2, and obtaining a comprehensive dynamic aging factor S data set according to a weight superposition mode; 5, estimating the voltage and the aging state of the fuel cell in real time to obtain a prediction data set of the aging state of the fuel cell; and step 6, estimating the remaining service life of the proton exchange membrane fuel cell. The problem that an existing battery life prediction method is inaccurate in prediction result is solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of fuel cell life prediction, and relates to a method for predicting the life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor. Background Art

[0002] Proton exchange membrane fuel cells are a new energy technology with great development potential because of their high energy conversion efficiency and environmentally friendly reaction products. However, due to the limitations of the internal structure and material properties of proton exchange membrane fuel cells themselves, the performance of current proton exchange membrane fuel cells decreases significantly during operation, resulting in a short service life of proton exchange membrane fuel cells and hindering the large-scale commercial development of proton exchange membrane fuel cells. If the fuel cell system can be predicted for failure in time before a fault occurs, providing a basis for formulating effective control and maintenance strategies, thereby extending the service life of the fuel cell and reducing the maintenance cost, it has great research significance for the development of proton exchange membrane fuel cell systems.

[0003] Currently, the prediction methods for the remaining service life of proton exchange membrane fuel cells usually adopt a hybrid method of model-driven and data-driven. The hybrid method combines a mechanism model, an empirical model or a semi-empirical and semi-mechanistic model with a machine learning algorithm to achieve the prediction of the remaining service life of proton exchange membrane fuel cells. The prediction effect of the hybrid method depends on the accuracy of the construction of the dynamic aging factor. However, the proton exchange membrane fuel cell system is a complex system with multiple physical domains and multiple time scales. The existing hybrid methods do not fully understand the aging mechanism of proton exchange membrane fuel cells, and there are problems with low accuracy in the constructed dynamic aging factor. Currently, there is no method for predicting the life of a proton exchange membrane fuel cell with a comprehensive dynamic aging factor considering both macro and micro aspects. By constructing a comprehensive dynamic aging factor for proton exchange membrane fuel cells to more comprehensively characterize the internal degradation of fuel cells, the historical state of the fuel cell system can be estimated more accurately and effectively, and the remaining service life of the fuel cell can be predicted.

[0004] Due to limited research, the existing learning of the aging mechanism of proton exchange membrane fuel cells is relatively one-sided, so there are great limitations in constructing a dynamic factor characterizing fuel cell aging, which affects the prediction effect of the remaining service life of proton exchange membrane fuel cells. Currently, there is no method for predicting the life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor combining micro and macro analysis. The dynamic aging factor used in fuel cell aging estimation cannot effectively and accurately predict the remaining service life. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for predicting the life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor, which solves the problem of inaccurate prediction results existing in the existing battery life prediction methods.

[0006] The technical solution adopted by the present invention is a method for predicting the life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor, which specifically includes the following steps:

[0007] Step 1, draw the electrochemical impedance spectroscopy diagrams at different times and the polarization curves at different times;

[0008] Step 2, construct an aging equivalent circuit model of the proton exchange membrane fuel cell, and extract a data set of dynamic factor S1 based on the equivalent circuit model;

[0009] Step 3, establish a semi-empirical equation, associate it with the equivalent circuit model to obtain an associated model, and extract a data set of dynamic factor S2 based on the associated model;

[0010] Step 4, normalize the dynamic factors S1 and S2, and obtain a data set of comprehensive dynamic aging factor S in the way of weighted superposition;

[0011] Step 5, estimate the voltage and aging state of the fuel cell in real time to obtain a prediction data set of the fuel cell aging state;

[0012] Step 6, realize the estimation of the remaining service life of the proton exchange membrane fuel cell.

[0013] The characteristics of the present invention also lie in:

[0014] The specific process of Step 1 is as follows: Step 1.1, obtain the impedance-frequency data set of the proton exchange membrane fuel cell stack through the impedance test of the proton exchange membrane fuel cell stack;

[0015] Step 1.2, use the moving average filtering algorithm to filter the data in the data set obtained in Step 1.1 to obtain the filtered values of the real part impedance and the imaginary part impedance after preprocessing;

[0016] Step 1.3, import the impedance data set after preprocessing in Step 1.2 and the frequency data set in Step 1.1 into MATLAB, and draw the electrochemical impedance spectroscopy diagrams of different single fuel cells, at different times, and under different currents;

[0017] Step 1.4, respectively test and obtain the proton exchange membrane fuel cell stack voltage data set, current, and working temperature data set through a voltage sensor, a current sensor, and a temperature sensor;

[0018] Step 1.5, use the moving average filtering algorithm to filter the data in the data set obtained in Step 1.4 to obtain the filtered values of the stack voltage, current, and working temperature after preprocessing;

[0019] Step 1.6: Import the preprocessed stack voltage and current data sets from Step 1.5 and the time data set from Step 1.4 into MATLAB, and plot the polarization curves at different times.

[0020] The specific process of Step 2 is as follows: Step 2.1: Connect the inductor element L Ohm and the resistor element R Ohm in series to represent the impedance Z Ohm part of the ohmic polarization loss inside the fuel cell. Connect the capacitor element C Ca and the resistor element R Ca in parallel to represent the impedance Z Ca part of the cathode polarization loss inside the fuel cell. After connecting the impedance Z Ca in series with the resistor element R Con , then connect it in parallel with the constant phase capacitor element CPE Con to represent the impedance Z Con part of the concentration polarization loss inside the fuel cell;

[0021] Step 2.2: Connect the impedance Z Ohm and the impedance Z Con constructed in Step 2.1 in series successively to obtain an equivalent circuit model representing the aging of the proton exchange membrane fuel cell;

[0022] Step 2.3: Design the equivalent circuit model constructed in Step 2.2 in ZView software, and import the impedance-frequency data set of the proton exchange membrane fuel cell stack after the filtering process in Step 1.2;

[0023] Step 2.4: Based on the equivalent circuit model designed in Step 2.3, respectively fit the electrochemical impedance spectra of different single fuel cells, at different times, and at different currents plotted in Step 1.3, and extract the characteristic parameters of the equivalent circuit model for different single fuel cells, at different times, and at different currents;

[0024] Step 2.5: At different times t, set the coupling and superposition weight of the characteristic parameters of the equivalent circuit at a constant current and different times to J, and calculate the dynamic aging factor S t data set;

[0025] Step 2.6: At the a-th different current, set the coupling and superposition weight of the characteristic parameters of the equivalent circuit under different current conditions to E, and calculate the dynamic aging factor S a data set;

[0026] Step 2.7: Through fusion calculation of the dynamic aging factors S a extracted for different operating currents in Step 2.6, obtain the final dynamic factor S1 data set based on the equivalent circuit model.

[0027] The specific process of step 2.4 is as follows: Step 2.4.1, under the conditions of equal current and consistent time, the number of single fuel cells constituting the proton exchange membrane fuel cell stack is N cell , use ZView software to check N cell The electrochemical impedance spectrum of each single fuel cell is fitted, and the characteristic parameters of the equivalent circuit model of different single fuel cells with constant current and equal time are directly extracted according to the feedback data of ZView software. wherein i represents the i-th single fuel cell, and the value range is i=1, 2, ..., N cell ; L Ohm-i and R Ohm-i represent the ohmic inductance data set and the ohmic resistance data set of the i-th single fuel cell respectively; C Ca-i and R Ca-i Respectively represent the polarization capacitance data set and polarization resistance data set of the i-th single fuel cell; R Con-i represents the concentration polarization resistance data set of the i-th single fuel cell; Q i and Respectively represent the double layer capacitance data set and the exponential parameter data set of the constant phase capacitance element of the i-th single fuel cell; [] T Represents a matrix full-time symbol; the fuel cell stack consists of N cell The characteristic parameters of the equivalent circuit model of the fuel cell stack are calculated by connecting the single fuel cells in series. Data set:

[0028]

[0029] in, Indicates that for i=1, 2, ..., N cell The symbol for the superposition and summation of the characteristic parameters of the equivalent circuit model of a single fuel cell; L Ohm and R Ohm Respectively represent the ohmic inductance and ohmic resistance of the fuel cell stack; C Ca and R Ca Respectively represent the polarization capacitance and polarization resistance of the fuel cell stack; R Con represents the concentration polarization resistance of the fuel cell stack; Q and They respectively represent the double layer capacitance and the exponential parameter of the constant phase capacitance element of the fuel cell stack;

[0030] Step 2.4.2, under the condition of equal current and different time, mark the characteristic parameters of the equivalent circuit model of the fuel cell stack at different times as F t, where \(t\) represents different moments, and its value range is \(t = 1, 2, \ldots, T\). According to the fuel cell electrochemical impedance spectroscopy diagrams at different moments drawn in Step 1.3, the electrochemical impedance spectra are successively fitted through ZView software, and the characteristic parameter \(F\) of the fuel cell equivalent circuit model at different times under a constant current is extracted. t , \(F\) t is consistent with the expansion of the characteristic parameter \(F\) in formula (1), \(F\) t represents \(F\) at different times;

[0031] In Step 2.4.3, under the condition of different currents, the characteristic parameter of the equivalent circuit model of the fuel cell stack with different current values is marked as \(F\) a , where \(a\) represents the number of different current values, and its value range is \(a = 1, 2, \ldots, A\). According to the fuel cell electrochemical impedance spectroscopy diagrams at different moments drawn in Step 1.3, the electrochemical impedance spectra of the fuel cell under different current conditions are successively fitted through ZView software, and the characteristic parameter \(F\) of the fuel cell equivalent circuit model operating at different current values is extracted. a , \(F\) a is consistent with the expansion of the characteristic parameter \(F\) in formula (1), \(F\) a represents \(F\) at different currents.

[0032] In Step 2.5, the coupling and superposition weight \(J\) of each characteristic parameter of the equivalent circuit at a constant current and different moments is calculated using the following formula (2):

[0033]

[0034] where \(w_1, w_2, \ldots, w_7\) are the main diagonal elements in the coupling and superposition weight \(J\), representing the weight factors of the proportions of the 7 characteristic parameters \(F\) t in the equivalent circuit characteristic parameter; the symbol “-” represents normalization processing; in formula (2) the symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (3) is:

[0035]

[0036] where \(v\) represents the sorting number of the characteristic parameter \(F\) t , and its value is \(v = 1, 2, \ldots, 7\); \(F\) t-v (t) and \(F\) t-v (t - 1) respectively represent the \(v\)-th characteristic parameter in \(F\) t at the \(t\) and \((t - 1)\) moments; \(w\) v represents the weight factor corresponding to the \(v\)-th characteristic parameter in the coupling and superposition weight \(J\).

[0037] In Step 2.6, the coupling superposition weight E of each characteristic parameter of the equivalent circuit under different current conditions is calculated using the following formula (5):

[0038]

[0039] where e1, e2, …, e7 are the main diagonal elements in the coupling superposition weight E, representing the weight factors of the proportions of the 7 characteristic parameters in the equivalent circuit characteristic parameter F a ; the symbol “-” represents normalization processing; in formula (5) the symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (6) is:[[]]

[0040]

[0041] where m represents the sorting number in the characteristic parameter F a , with the value range of m = 1, 2, …, 7; F a-m (a) and F a-m (a - 1) respectively represent the m-th characteristic parameter in F a under the current values of a and (a - 1); w m represents the weight factor corresponding to the m-th characteristic parameter in the coupling superposition weight E.

[0042] The specific process of Step 2.7 is to fuse and calculate the dynamic aging factors S a extracted in Step 2.6 under different working currents through formula (8) and formula (9) to obtain the final data set of the dynamic factor S1 based on the equivalent circuit model:

[0043]

[0044] where j1, j2, …, jA are the main diagonal elements in the coupling superposition weight M, representing the weight factors of the proportions of the A dynamic factors in the dynamic factor S a ; the symbol “-” represents normalization processing; in formula (8) the ▽ symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (9) is:[[]]

[0045]

[0046] where S a (t) and S a (t - 1) respectively represent the dynamic factors at the t-th and (t - 1)-th moments under the a-th current value; j a represents the weight factor corresponding to the dynamic factor under the a-th current value in the coupling superposition weight M.

[0047] The specific process of Step 3 is:[[]]

[0048] Step 3.1, construct a semi-empirical model of the proton exchange membrane fuel cell, as shown in the following formula (11):

[0049] U stack = N cell ·(E OCV - U Ohm - U Con - U Ca ) (11)

[0050] Among them, U stack is the stack voltage, N cell is the number of single fuel cells in the stack, E OCV is the open circuit voltage at rated pressure and rated temperature, U Ohm is the ohmic polarization voltage, U Con is the concentration polarization voltage, U Ca is the polarization voltage;

[0051] In the semi-empirical model, the expression of the ohmic polarization voltage obtained from Ohm's law is:

[0052] U Ohm = R eq ·i stack (12)

[0053] Among them, R eq is the total resistance, i stack is the load current density;

[0054] In the semi-empirical model, the expression of the concentration polarization voltage U Con is:

[0055]

[0056] Among them, λ is the charge transfer coefficient, i limit is the limiting current density, R is the ideal gas constant, T work is the stack temperature, with a value of T0 °C, n e is the number of electrons participating in the electrochemical reaction, F ara is the Faraday constant;

[0057] In the semi-empirical model, the expression of the polarization voltage U Ca is:

[0058]

[0059] Among them, i change is the varying current density, i leak is the leakage current density. In equations (13) and (14), λ is a constant. In equation (14), ileak much smaller than i stack , i leak is ignored; finally, Equation (11) is simplified to:

[0060]

[0061] where A ta is the Tafel constant and B con is the concentration constant;

[0062] Step 3.2: Combine the fuel cell semi-empirical model and the equivalent circuit model to obtain a macro-micro combined correlation model for the proton exchange membrane fuel cell, as shown in the following formula (16):

[0063] U stack = N cell ·(E OCV - M Ohm - M Ca - M Con ) (16)

[0064] where M Ohm , M Ca and M Con represent the ohmic polarization loss, activation polarization loss, and concentration polarization loss in the model that correlates the equivalent circuit characteristic parameters with the semi-empirical equation, respectively, and are specifically expressed as:

[0065]

[0066] where S area is the cross-sectional area through which the operating current passes, with the unit of A / cm 2 ; |·| is the absolute value symbol, i.e., modulus calculation; Q Ohm , Q Ca and Q Con are the correlation coefficients of the ohmic polarization loss, activation polarization loss, and concentration polarization loss in the correlation model, respectively; Z Ohm , Z Ca and Z Con are the ohmic polarization impedance, activation polarization impedance, and concentration polarization impedance in the equivalent circuit model, respectively;

[0067] Step 3.3: Based on the correlation model, use the non-linear least squares method to fit the polarization curves at different times, and capture the characteristic parameters of the correlation model at different times H = [E OCV , Q Ca , Q Ohm , Q Con T ;

[0068] ​Step 3.4, at different times, set the coupling superposition weight of each characteristic coefficient of the correlation model at different times as D, and calculate the correlation coefficient S2 of the correlation model through the characteristic parameter H of the correlation model:

[0069]

[0070]

[0071] Among them, d1, d2, d3, and d4 are the main diagonal elements in the coupling superposition weight D, representing the weight factors of the proportion of characteristic parameters in the characteristic parameters of the correlation model; the symbol "-" represents normalization processing; the X symbol in formula (18) represents a specific decoupling calculation method, and the specific algorithm combined with formula (19) is:

[0072]

[0073] Among them, l represents the l-th element in the characteristic parameter H of the correlation model, and the value is l = 1, 2, 3, 4; d l represents the weight factor corresponding to the l-th correlation model coefficient in the coupling superposition weight H; H l (t) and H l (t - 1) respectively represent the l-th characteristic parameter in the characteristic parameter H of the correlation model at times t and (t - 1).

[0074] The specific process of Step 4 is as follows:

[0075] Step 4.1, perform conventional normalization processing on the dynamic factor S1 data group and the dynamic factor S2 data group, and map the original data group to the dimension interval of [0, 1];

[0076] Step 4.2, at different times, set the coupling superposition weight of the dynamic factor based on the equivalent circuit model and the dynamic factor based on the correlation model at different times as C, and perform weight superposition on the dynamic factor S1 data group and the dynamic factor S2 data group after the normalization processing in Step 4.1 to calculate the comprehensive dynamic aging factor S data group:

[0077]

[0078] Among them, c1 and c2 are the main diagonal elements in the coupling superposition weight C, representing the weight factors of the proportion of the dynamic factor based on the equivalent circuit model and the dynamic factor based on the correlation model; the symbol "-" represents normalization processing; in formula (21) the symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (22) is:

[0079]

[0080] Among them, k represents the sequence of dynamic factors, and the value range is k = 1, 2; c k represents the weight factor corresponding to the k-th dynamic factor in the coupling superposition weight C; S k (t) and S k (t - 1) respectively represent the k-th dynamic factor at times t and (t - 1). When k = 1, S k represents the dynamic factor based on the equivalent circuit model; when k = 2, S k represents the dynamic factor based on the correlation model.

[0081] The specific process of step 5 is as follows:

[0082] Step 5.1: Based on the extracted comprehensive dynamic aging factor S data set, using the discrete nonlinear system and the associated particle filter algorithm, taking the comprehensive dynamic aging factor as the state variable of the system state equation, and according to formula (24), complete the state estimation of the system parameter vector T including the equivalent circuit model and the joint model for 1 time step:

[0083] T i = f i (T i-1 , ω i-1 ) (24)

[0084] Among them, i is the time step; T is the system parameter vector including the equivalent circuit model and the joint model; f(·) is the degradation model describing the state change; ω i is the system noise;

[0085] Step 5.2: Taking the stack voltage data of the proton exchange membrane fuel cell as the observation variable in the observation equation, and the current and operating temperature data of the proton exchange membrane fuel cell as the input variables in the observation equation, using the system state parameters estimated in step 5.1, and according to formula (25) and formula (15), complete the observation estimation of the stack voltage of the proton exchange membrane fuel cell for 1 time step:

[0086] z t = h i ((T t-1 , x t ), v t ) (25)

[0087] Among them, z i is the observation state at time i, T i-1 is the system state at time (i - 1), x t represents the input load current and operating temperature that can be planned in advance, v t represents the observation noise, and h(·) is the observation equation;

[0088] Step 5.3, compare the actual values of the stack voltage data group with the estimated values of the stack voltage data group estimated in Step 5.2, and achieve the correction of the nonlinear system through the adjustment of importance weights. Calculate the weights according to the observation likelihood function, and the likelihood function is usually defined in the Gaussian form:

[0089]

[0090] where r is the number of particles, and the value range is r = 1, 2, …, R; T is the matrix transpose symbol; V is the observation noise covariance matrix, representing the uncertainty of the observed values; exp represents the exponential operation with the natural constant e as the base;

[0091] Step 5.4, repeat Step 5.1 to Step 5.3 until the estimation of the stack voltage data group is completed;

[0092] Step 5.5, end the estimation, and output the stack voltage data group of the observed parameters estimated in real time and the comprehensive dynamic aging factor data group of the output state parameters estimated in real time.

[0093] The beneficial effects of the present invention are as follows: When predicting the aging of a proton exchange membrane fuel cell system, the present invention uses a comprehensive dynamic aging factor as an evaluation index, and can directly obtain information on the historical or future fuel cell health state through sampling data. The present invention extracts the comprehensive dynamic aging factor through the combination of microscopic analysis and macroscopic analysis, and can intuitively give the change of the aging state of the proton exchange membrane fuel cell with accurate prediction. The estimated value of the historical comprehensive dynamic degradation factor constructed by the present invention conforms to the monotonic change trend of the experimental data, and can reflect the monotonic change of the historical comprehensive dynamic degradation factor. When the present invention estimates the remaining service life of a proton exchange membrane fuel cell in real time, it applies the particle filter algorithm to estimate the future comprehensive dynamic aging factor, and the estimation result is reasonable. The present invention considers the requirements for the fuel cell life threshold in practical applications, and conducts long-term prediction on the aging state of the fuel cell through a neural network, and can achieve a more accurate and effective prediction effect on the service life of the fuel cell. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 It is a schematic diagram of the electrochemical impedance spectrum in the proton exchange membrane fuel cell life prediction method based on the comprehensive dynamic aging factor of the present invention;

[0095] Figure 2 It is a schematic diagram of the polarization curve in the proton exchange membrane fuel cell life prediction method based on the comprehensive dynamic aging factor of the present invention;

[0096] Figures 3(a) to 3(c)Schematic diagrams of partial structures of equivalent circuits for characterizing ohmic polarization loss, polarization loss, and concentration polarization loss inside a fuel cell in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention; Fig. 3(d) is a schematic diagram of the completely constructed equivalent circuit in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention;

[0097] Figure 4 Result graph for extracting the comprehensive dynamic aging factor in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention;

[0098] Fig. 5(a) and Fig. 5(b) are result graphs of real-time estimation of the operating condition stack voltage and the comprehensive dynamic aging factor in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention;

[0099] Figure 6 Schematic diagram of the structure for constructing an improved neural network in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention;

[0100] Fig. 7(a) to (f) are result graphs of long-term prediction of the degradation trend in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention, which are result graphs for training sets of 40%, 50%, 60%, 70%, 80%, and 90% respectively;

[0101] Figure 8 Result graph for estimating the remaining service life in the proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention. Detailed implementation manners

[0102] The present invention will be described in detail below with reference to the accompanying drawings and specific implementation manners.

[0103] Example 1

[0104] The proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to the present invention specifically includes the following steps: Step 1, through the impedance (including real part and imaginary part)-frequency data set of the proton exchange membrane fuel cell stack, draw the electrochemical impedance spectroscopy diagrams of different single fuel cells, at different times, and at different currents, as Figure 1 shown. Through the current and voltage data sets of the proton exchange membrane fuel cell stack, draw the polarization curves at different times, as Figure 2 shown.

[0105] Step 2: Construct an equivalent circuit model for the aging of a proton exchange membrane fuel cell, as shown in Figure 3. Fit the electrochemical impedance spectrogram in Step 1 through software ZView, capture the characteristic parameters of the fuel cell aging equivalent circuit model, and extract the data set of the dynamic factor S1 based on the equivalent circuit model.

[0106] Step 3: Establish a semi-empirical equation and associate it with the equivalent circuit model constructed in Step 2 to obtain an associated model. Use the nonlinear least squares method to fit the polarization curve in Step 1, capture the dynamic coefficient of the associated model, and extract the data set of the dynamic factor S2 based on the associated model.

[0107] Step 4: Normalize the dynamic factor S1 extracted in Step 2 and the dynamic factor S2 extracted in Step 3, and obtain the data set of the comprehensive dynamic aging factor S by the method of weighted superposition.

[0108] Step 5: Take the stack current and operating temperature obtained in Step 1 as the inputs of the particle filter algorithm, take the stack voltage obtained in Step 1 and the aging state extracted in Step 4 as the outputs, take the semi-empirical model as the observation equation, set the initial values of the operating parameters and noise in the particle filter algorithm, and estimate the voltage and aging state of the fuel cell in real time to obtain the prediction data set of the fuel cell aging state.

[0109] Step 6: Select the fuel cell aging state prediction data set in Step 5 as the input, set the optimal parameter values of the leakage rate, spectral radius, and regularization factor of the 4 sub-reservoirs through the training results of historical data, and obtain a decoupled echo state network. Perform long-term prediction on the fuel cell aging trend through the trained neural network to realize the estimation of the remaining service life of the proton exchange membrane fuel cell.

[0110] Example 2

[0111] Specifically, Step 1 includes: Step 1-1, the impedance (including real part and imaginary part)-frequency data set of the proton exchange membrane fuel cell stack used in the present invention is obtained through impedance testing of the proton exchange membrane fuel cell stack, and the data covers the test results under two constant current conditions of 20 A and 85 A. Each set of constant current test data contains the impedance data (real part and imaginary part) of 8 single fuel cells, and each single fuel cell is tested at 4 time points of 0 h, 120 h, 240 h, and 360 h. At each time point, each single fuel cell contains 45 groups of data, and each group of data consists of the actual values of frequency, real part impedance, and imaginary part impedance.

[0112] Step 1-2: Use the sliding average filtering algorithm to filter the impedance (including real part and imaginary part) data set obtained for each single fuel cell at each time point in Step 1-1 to obtain the filtered values of the real part impedance and imaginary part impedance after preprocessing.

[0113] Steps 1-3: Import the real part impedance and imaginary part impedance data sets pre-processed in Step 1-2 and the frequency data set in Step S1-1 into MATLAB, and plot the electrochemical impedance spectra of different single fuel cells, at different times, and under different currents. The electrochemical impedance spectrum data of the 7th single fuel cell and the 8th single fuel cell connected in series at different times under the constant current condition of 50 A are as Figure 1 shown.

[0114] Steps 1-4: The proton exchange membrane fuel cell stack voltage data set, current, and operating temperature data set used in the present invention are respectively obtained by testing with a voltage sensor, a current sensor, and a temperature sensor. The stack voltage, current, and operation of the fuel cell are tested at a total of 10 time points: 0 h, 85 h, 150 h, 200 h, 250 h, 300 h, 350 h, 400 h, 450 h, and 500 h. At each time point, the fuel cell stack contains 25 groups of data, and each group of data consists of the actual values of time, voltage, current, and operating temperature.

[0115] Steps 1-5: Adopt a moving average filtering algorithm to filter the fuel cell stack voltage, current, and operating temperature obtained in Step 1-4 at each time point to obtain the filtered values of the pre-processed stack voltage, current, and operating temperature.

[0116] Steps 1-6: Import the pre-processed stack voltage and current data sets in Step 1-5 and the time data set in Step S1-4 into MATLAB, and plot the polarization curves at different times. The polarization curves of the fuel cell stack at different times are as Figure 2 shown.

[0117] Example 3

[0118] The specific process of Step 2 is as follows: Step 2-1: Connect the inductance element L Ohm and the resistance element R Ohm in series to represent the impedance Z Ohm part of the ohmic polarization loss inside the fuel cell, as shown in Fig. 3(a). Connect the capacitance element C Ca and the resistance element R Ca in parallel to represent the impedance Z Ca part of the cathode polarization loss inside the fuel cell, as shown in Fig. 3(b). After connecting the impedance Z Ca in series with the resistance element R Con , then connect it in parallel with the constant phase capacitance element CPE Con to represent the impedance Z Con part of the concentration polarization loss inside the fuel cell, as shown in Fig. 3(c).

[0119] Steps 2-2: The impedance Z constructed in Step 2-1Ohm and impedance Z Con (Including impedance Z Ca ) are connected in series in sequence to obtain an equivalent circuit model that characterizes the aging of proton exchange membrane fuel cells, as shown in Figure 3(d).

[0120] Step 2-3, design the equivalent circuit model constructed in step 2-2 in ZView software, and import the impedance (including real and imaginary parts)-frequency data set of the proton exchange membrane fuel cell stack after filtering in step S1-2.

[0121] Step 2-4, based on the equivalent circuit model designed in step 2-3, respectively fit the electrochemical impedance spectra of different single fuel cells, at different times and at different currents drawn in step 1-3, and extract the characteristic parameters of the equivalent circuit model of different single fuel cells, at different times and at different currents.

[0122] Step 2-4-1, under the conditions of equal current and consistent time, the number of single fuel cells constituting the proton exchange membrane fuel cell stack is N cell , use ZView software to check N cell According to the feedback data from ZView software, the characteristic parameters of the equivalent circuit model of different single fuel cells with constant current and equal time are directly extracted. data group.

[0123] Where i represents the i-th single fuel cell, and the value range is i=1, 2, ..., N cell ; L Ohm-i and R Ohm-i represent the ohmic inductance data set and the ohmic resistance data set of the i-th single fuel cell respectively; C Ca-i and R Ca-i Respectively represent the polarization capacitance data set and polarization resistance data set of the i-th single fuel cell; R Con-i represents the concentration polarization resistance data set of the i-th single fuel cell; Q i and Respectively represent the double layer capacitance data set and the exponential parameter data set of the constant phase capacitance element of the i-th single fuel cell; [] T Represents the matrix full-time symbol.

[0124] The fuel cell stack is composed of N cell The characteristic parameters of the equivalent circuit model of the fuel cell stack are calculated by connecting the single fuel cells in series. Data set:

[0125]

[0126] in, Denotes the summation symbol for the characteristic parameters of the equivalent circuit model of the \(i = 1, 2, \ldots, N\) cell single fuel cells; \(L\) Ohm and \(R\) Ohm respectively represent the ohmic inductance and ohmic resistance of the fuel cell stack; \(C\) Ca and \(R\) Ca respectively represent the polarization capacitance and polarization resistance of the fuel cell stack; \(R\) Con represents the concentration polarization resistance of the fuel cell stack; \(Q\) and respectively represent the double-layer capacitance and exponential parameter of the constant phase capacitance element of the fuel cell stack.

[0127] Step 2-4-2, under the condition of equal current and different times, mark the characteristic parameters of the equivalent circuit model of the fuel cell stack at different times as \(F\) t , where \(t\) represents different times, and the value range is \(t = 1, 2, \ldots, T\). According to the fuel cell electrochemical impedance spectra at different times drawn in step S1-3, use ZView software to fit the electrochemical impedance spectra in turn, and extract the characteristic parameters \(F\) of the fuel cell equivalent circuit model at different times with a constant current t . \(F\) t is consistent with the expansion of the characteristic parameter \(F\) in formula (1), and \(F\) t represents \(F\) at different times.

[0128] Step 2-4-3, under the condition of different currents, mark the characteristic parameters of the equivalent circuit model of the fuel cell stack with different current values as \(F\) a , where \(a\) represents the number of different current values, and the value range is \(a = 1, 2, \ldots, A\). According to the fuel cell electrochemical impedance spectra at different times drawn in step S1-3, use ZView software to fit the fuel cell electrochemical impedance spectra under different current conditions in turn, and extract the characteristic parameters \(F\) of the fuel cell equivalent circuit model operating at different current values a . \(F\) a is consistent with the expansion of the characteristic parameter \(F\) in formula (1), and \(F\) a represents \(F\) at different currents.

[0129] Step 2-5, among different times \(t\), set the coupling superposition weight of the characteristic parameters of the equivalent circuit at a constant current and different times as \(J\), and use the characteristic parameter \(F\) t extracted in step 2-4-2 to calculate the dynamic aging factor \(S\) when the current is equal and the time is different t Data set:

[0130]

[0131] Among them, w1, w2, …, w7 are the main diagonal elements in the coupled superposition weight J, representing the weight factors of the proportion of the seven characteristic parameters in the equivalent circuit characteristic parameter F at different times; the symbol "-" represents the normalization process; in formula (2) t The symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (3) is:

[0132]

[0133] Among them, v represents the sorting number of the characteristic parameter F t in, and the value range is v = 1, 2, …, 7; F t-v (t) and F t-v (t - 1) respectively represent the v-th characteristic parameter in the characteristic parameter F at times t and (t - 1); w t represents the weight factor corresponding to the v-th characteristic parameter in the coupled superposition weight J v

[0134] Step 2 - 6, at the a-th different current, set the coupled superposition weight of each characteristic parameter of the equivalent circuit under different current conditions to E, and calculate the dynamic aging factor S with different working currents through the characteristic parameter F a extracted by step 2 - 4 - 3 a Data set:

[0135]

[0136]

[0137] Among them, e1, e2, …, e7 are the main diagonal elements in the coupled superposition weight E, representing the weight factors of the proportion of the seven characteristic parameters in the equivalent circuit characteristic parameter F at different currents; the symbol "-" represents the normalization process; in formula (5) a The symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (6) is:

[0138]

[0139] Among them, m represents the sorting number of the characteristic parameter F a in, and the value range is m = 1, 2, …, 7; F a-m (a) and F a-m (a - 1) respectively represent the m-th characteristic parameter in the characteristic parameter F at current values a and (a - 1); w a represents the weight factor corresponding to the m-th characteristic parameter in the coupled superposition weight E m

[0140] ​​Step 2-7: The dynamic aging factors S of different working currents extracted in Step 2-6 a are fused and calculated through Formulas (8) and (9) to obtain the final data set of dynamic factor S1 based on the equivalent circuit model:

[0141]

[0142] where j1, j2, …, jA are the main diagonal elements in the coupling superposition weight M, representing the weight factors of the proportions of A dynamic factors in the dynamic factor S a ; the symbol "-" represents normalization processing; the ▽ symbol in Formula (8) represents a specific decoupling calculation method, and the specific algorithm combined with Formula (9) is:

[0143]

[0144] where S a (t) and S a (t - 1) respectively represent the dynamic factors at times t and (t - 1) under the a-th current value; j a represents the weight factor in the coupling superposition weight M corresponding to the dynamic factor under the a-th current value.

[0145] Example 4

[0146] Specifically included in Step 3: Step 3-1, constructing a semi-empirical model of a proton exchange membrane fuel cell:

[0147] U stack = N cell ·(E OCV - U Ohm - U Con - U Ca ) (11)

[0148] where U stack is the stack voltage, N cell is the number of single fuel cells in the stack, E OCV is the open circuit voltage under rated pressure and rated temperature, U Ohm is the ohmic polarization voltage, U Con is the concentration polarization voltage, U Ca is the polarization voltage. In the semi-empirical model, the expression of the ohmic polarization voltage obtained from Ohm's law is:

[0149] U Ohm = R eq ·i stack (12)

[0150] where R eq is the total resistance, i stackis the load current density. In the semi-empirical model, the concentration polarization voltage U Con has the following expression:

[0151]

[0152] where λ is the charge transfer coefficient, i limit is the limiting current density, R is the ideal gas constant with a value of 8.314 J / (mol·K), T work is the stack temperature with a value of T0 °C, n e is the number of electrons participating in the electrochemical reaction with a value of 2, and F ara is the Faraday constant with a value of 96485 C / mol.

[0153] In the semi-empirical model, the polarization voltage U Ca has the following expression:

[0154]

[0155] where i change is the varying current density, and i leak is the leakage current density, which is related to gas permeation and circuit short-circuit. In Equations (13) and (14), λ is strictly stable and is usually regarded as a constant. In Equation (14), i leak is much smaller than i stack , and i leak is usually ignored. Finally, Equation (11) can be simplified to:

[0156]

[0157] where A ta is the Tafel constant, and B con is the concentration constant.

[0158] Step 3-2: Combine the fuel cell semi-empirical model constructed in Step 3-1 and the equivalent circuit model constructed in Step 2 to obtain a macro-micro combined correlation model of the proton exchange membrane fuel cell:

[0159] U stack = N cell ·(E OCV - M Ohm - M Ca - M Con ) (16)

[0160] where M Ohm , M Ca and M Con respectively represent the ohmic polarization loss, activation polarization loss, and concentration polarization loss in the model of the correlation equivalent circuit characteristic parameters and the semi-empirical equation, and are specifically expressed as:

[0161]

[0162] Among them, S area is the cross-sectional area through which the operating current passes, with the unit of A / cm 2 ; |·| is the absolute value symbol, i.e., modulus calculation; Q Ohm , Q Ca and Q Con are the correlation coefficients of ohmic polarization loss, activation polarization loss, and concentration polarization loss in the correlation model, respectively; Z Ohm , Z Ca and Z Con are the ohmic polarization impedance, activation polarization impedance, and concentration polarization impedance in the equivalent circuit model constructed in step 2-1, respectively.

[0163] Step 3-3, based on the correlation model constructed in step 3-2, use the non-linear least squares method to fit the polarization curves at different times drawn in step 1-6, and capture the characteristic parameters H = [E OCV , Q Ca , Q Ohm , Q Con of the correlation model at different times. T .

[0164] Step 3-4, at different times, set the coupling superposition weight of each characteristic coefficient of the correlation model at different times to D, and calculate the correlation coefficient S2 of the correlation model through the characteristic parameter H of the correlation model in step 3-3:

[0165]

[0166] Among them, d1, d2, d3, and d4 are the main diagonal elements in the coupling superposition weight D, representing the weight factors of the proportion of characteristic parameters in the characteristic parameters of the correlation model; the symbol "-" represents normalization processing; the symbol X in formula (18) represents a specific decoupling calculation method, and the specific algorithm combined with formula (19) is:

[0167]

[0168] Among them, l represents the l-th element in the characteristic parameter H of the correlation model, with the value of l = 1, 2, 3, 4; d l represents the weight factor corresponding to the l-th correlation model coefficient in the coupling superposition weight H; H l (t) and H l (t - 1) represent the l-th characteristic parameter in the characteristic parameter H of the correlation model at times t and (t - 1), respectively.

[0169] Example 5

[0170] Specifically, Step 4 includes: Step 4-1, performing conventional normalization processing on the dynamic factor S1 data group extracted in Step 2 and the dynamic factor S2 data group extracted in Step 3, and mapping the original data group to the dimensional interval of [0, 1].

[0171] Step 4-2, at different times, setting the coupling superposition weight of the dynamic factor based on the equivalent circuit model and the dynamic factor based on the correlation model at different times as C, and performing weight superposition on the dynamic factor S1 data group and the dynamic factor S2 data group after the normalization processing in Step 4-1 to calculate the comprehensive dynamic aging factor S data group, as Figure 4 shown:

[0172]

[0173] where c1 and c2 are respectively the main diagonal elements in the coupling superposition weight C, representing the weight factors of the proportion of the dynamic factor based on the equivalent circuit model and the dynamic factor based on the correlation model; the symbol "-" represents the normalization processing; in formula (21) the symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (22) is:

[0174]

[0175] where k represents the sequence of the dynamic factor, and the value is k = 1, 2; c k represents the weight factor corresponding to the kth dynamic factor in the coupling superposition weight C; S k (t) and S k (t - 1) respectively represent the kth dynamic factor at times t and (t - 1). When k = 1, S k represents the dynamic factor based on the equivalent circuit model; when k = 2, S k represents the dynamic factor based on the correlation model.

[0176] Example 6

[0177] Specifically, Step 5 includes: Step 5-1, taking the comprehensive dynamic aging factor S with non-linear change extracted in Step S4 as the discrete non-linear system in the particle filter algorithm. Therefore, taking the comprehensive dynamic aging factor as the state variable of the system state equation, which is represented in the form of a particle set. Completing the state estimation of the system parameter vector T of the equivalent circuit model constructed in Step 2 and the joint model constructed in Step 3 for 1 time step, and realizing the state estimation of the comprehensive aging factor particle set: The estimation expression is:

[0178] T i = f i (T i-1 , ω i-1 ) (24)

[0179] where \(i\) is the time step; \(T\) is the system parameter vector including the equivalent circuit model constructed in step 2 and the combined model constructed in step 3; \(f(\cdot)\) is the degradation model describing the state change; \(\omega\) i is the system noise, and the noise state variables follow a normal distribution and are independent of each other.

[0180] Step 5-2: Take the stack voltage data of the proton exchange membrane fuel cell in steps 1-4 as the observation variable in the observation equation, and the current and operating temperature data of the proton exchange membrane fuel cell in steps 1-4 as the input variables in the observation equation. Using the system state parameters estimated in step 5-1, substitute formula (15) into formula (25) to calculate the observation estimate of the stack voltage of the proton exchange membrane fuel cell for 1 time step:

[0181] z t =h i ((T t-1 ,x t ),v t ) (25)

[0182] where \(z\) i is the observed state at time \(i\), that is, the stack voltage measurement value in steps 1-4; \(T\) i-1 is the system state at time \((i - 1)\), that is, the system parameter vector of the equivalent circuit model constructed in step 2 and the combined model constructed in step 3; \(x\) t represents the input load current and operating temperature that can be pre-planned, that is, the current and operating temperature data of the proton exchange membrane fuel cell in steps 1-4; \(v\) t represents the observation noise, the noise state variables follow a normal distribution and are independent of each other, and \(h(\cdot)\) is the observation equation, that is, formula (15).

[0183] Step 5-3: Compare the actual values of the stack voltage data set in steps 1-4 with the estimated values of the stack voltage data set estimated in step 5-2, and calculate the weights of each particle in the non-linear system particle set. The weight calculation uses the observation likelihood function, and the likelihood function is usually defined in a Gaussian form:

[0184]

[0185] where \(\omega\) is the weight of each particle in the non-linear system particle set; \(r\) is the number of particles, and the value range is \(r = 1, 2, \ldots, R\); \(T\) is the matrix transpose symbol; \(V\) is the observation noise covariance matrix, representing the uncertainty of the observed value; \(\exp\) represents the exponential operation with the natural constant \(e\) as the base. Normalize the weights of each particle after calculation:

[0186]

[0187] Where Nu is the total number of particles in the particle set; jr represents the jr-th particle in the particle set. Construct a cumulative weight vector (with a length of Nu) for random sampling of the particle set:

[0188]

[0189] Where Cu is the constructed cumulative weight vector. Generate Nu independent random numbers from the uniform distribution U(0,1), and then determine the source of each new particle through inverse transform sampling:

[0190]

[0191] Where x r represents the comprehensive aging factor state variable corresponding to the r-th particle; u r represents the Nu independent random numbers generated by the r-th in the uniform distribution U(0,1). Repeat formula (29) Nu times to obtain the resampled particle set and complete the state correction.

[0192] Step 5-4, repeat Step S5-1 to Step S5-3 until the fuel cell stack voltage data set is estimated.

[0193] Step 5-5, end the estimation, output the fuel cell stack voltage data set of the observed parameters estimated in real time, the estimation result is shown in Fig. 5(a), and the comprehensive dynamic aging factor data set of the output state parameters estimated in real time, the estimation result is shown in Fig. 5(b). The result in Fig. 5(a) shows that by using the particle filter algorithm, the fuel cell stack voltage can be estimated in real time through the input current and operating temperature; the result in Fig. 5(b) shows that by using the particle filter algorithm, the change of the comprehensive aging index of the fuel cell can be estimated in real time through the input current and operating temperature.

[0194] Example 7

[0195] Specifically included in Step 6: Step 6-1, divide the reservoir with Res neurons into 4 sub-reservoirs, each sub-reservoir has Res / 4 neurons, connect the 4 sub-reservoirs in parallel in turn, and add an inhibitory coupling mechanism to obtain an improved decoupled echo state network structure, as Figure 6 shown. In the improved decoupled echo state network structure, set the number of input nodes of the sub-reservoir to L, the number of state nodes of the sub-reservoir to M, and the number of output nodes to N. The corresponding update state equation of the sub-reservoir is:

[0196]

[0197] Where, and respectively represent the activation vector and the corresponding updated state of the neurons in the reservoir at time t; the activation function of the neurons is f(·), usually tanh(·); is the input weight matrix, is the reservoir recurrent weight matrix; is the input vector at time (t - 1), is the activation vector at time (t - 1); the reservoir leakage rate α ∈ (0, 1]. In the improved decoupled echo state network structure, the output state equation of each sub-reservoir is as follows:

[0198] q(t) = S out *[p(t - 1); k(t)] (31)

[0199] where, is the output vector; is the output weight matrix; the root mean square error is calculated by ridge regression to obtain S out . Equation (32) defines the training objective of the output layer weight S out , by minimizing the root mean square error between the network output and the target value, ensuring that the output approaches the expected value, which is the theoretical expression of the training objective; Equation (33) directly gives the analytical solution of the output weight, solved by ridge regression, which is the theoretical expression of the training objective.

[0200]

[0201] S out = Q target V out T (V out V out T +γI) -1 (33)

[0202] where, D is the number of all data points in the training dataset, V out is the reservoir output matrix, Q target is the target value output matrix; γ is the regularization parameter. I is the identity matrix. A suppression hidden mechanism is added to the reservoir to obtain the reservoir recurrent weight and the updated states of 4 sub-reservoirs. The reservoir recurrent weight is calculated as:

[0203]

[0204] where, is the recurrent weight matrix of the reservoir in the improved decoupled echo state network structure; is the identity matrix; in the circulant matrix, is the cyclic weight of 4 sub-reservoirs. Substituting formula (34) into formula (30), the calculation expressions for the updated state and the current state of the reservoir in the decoupled echo state network can be obtained as follows:

[0205]

[0206] Among them, the updated state of the entire reservoir at time t is composed of the updated states of 4 sub-reservoirs is the current state of the reservoir at time t; meanwhile, 4 identical inputs constitute a new 4-column input matrix which is the input matrix of the entire reservoir; is the input matrix of 4 sub-reservoirs; is the state matrix of the entire reservoir; is the state matrix of 4 sub-reservoirs; is the identity matrix; the entire reservoir leakage rate parameter matrix constituted by the 4 sub-reservoir leakage rate parameter matrices Γ is defined as a prediction operation symbol used to calculate the associated state of the sub-reservoir at time t:

[0207] Γ·k i (t - 1) = f(S in(i) *p(t) + S i *k i (t - 1)) (37)

[0208] Substitute formula (37) into formula (34), and then combine with formula (36) to calculate the entire reservoir cyclic weight matrix S R , and then combine with formula (31), formula (32) and formula (33) to calculate the output matrix of 4 sub-reservoirs Then calculate the outputs of 4 sub-reservoirs:

[0209] q i (t) = S out(i) *[p(t - 1); k i (t)] (38)

[0210] Among them, is the output of 4 sub-reservoirs. According to the output weights, the outputs of 4 sub-reservoirs are superimposed to obtain the final output of the decoupled echo state network structure:

[0211]

[0212] Among them, represents the superimposed value of the inputs of sub-reservoir 1 and sub-reservoir 2; Represents the input superposition value of sub-reservoirs 3 and 4; Represents And The superposition value of; w1, w2, and w3 respectively represent the weight ratios of sub-reservoir 1, sub-reservoir 3, and the superposition value q 12 (t), obtaining the calculation framework for running a decoupled echo state network.

[0213] Step 6-2, set the parameter range of the decoupled echo state network structure constructed in Step 6-1, use the comprehensive dynamic aging factor data group estimated in Step 5 as the input of the neural network structure, and then use the particle swarm algorithm to optimize the parameters of the decoupled echo state network structure constructed in Step 6-1.

[0214] Step 6-2-1, set the parameter range of the decoupled echo state network structure constructed in Step 6-1: leakage rate α m ∈(0, 1], spectral radius ρ m ∈(0, 1.5], regularization coefficient γ m ∈(0, 1] (the number of the m-th sub-reservoir m = 1, 2, 3, 4), superposition weights w1, w2, w3 ∈(0, 1], a total of 15 parameters to be optimized.

[0215] Step 6-2-2, first, use the comprehensive dynamic aging factor data group estimated in Step 5 as the input of the neural network structure, set the number of particles as Num, input the range of 15 parameters to be optimized, and randomly initialize the position and velocity of each particle.

[0216] Step 6-2-3, then, calculate the fitness value of each particle, and set the fitness function as the root mean square error to evaluate the optimization result:

[0217]

[0218] Among them, Fitness RMSE Is the fitness function, and the calculation form is the root mean square error, which is calculated through the comprehensive dynamic aging factor data group in Step S5 and the data group of long-term prediction of the neural network constructed in Step S6-1; f(t) is the actual value of the comprehensive aging index; f’(t) is the predicted value of the comprehensive aging index; F is the number of prediction results.

[0219] Step 6-2-4, continuously iterate and optimize. In each iteration, update the velocity and position of the particle, and evaluate the new fitness value according to formula (40). The velocity update expression is:

[0220] v i t+1 = w·v i t+ c1r1(pbest i - x i t ) + c2r2(gbest - x i t ) (41)

[0221] where, at time t, x i is the position of the particle, the coordinate in the search space, i.e., the candidate solution; v i is the velocity of the particle; pbest i is the individual historical best particle; gbest i is the global best particle of the group; c1 and c2 are the learning factors of individual cognition and social learning respectively; w is the inertia weight, which can be linearly decreased to enhance the late local search; r1, r2 ~ U(0,1) are random numbers to increase the exploration randomness. The update expression is:

[0222] x i t+1 = x i t + v i t+1 (42)

[0223] Step 6 - 2 - 5, finally, when the maximum number of iterations T max or the improvement of gbest is less than the preset threshold, the iteration is completed, the global optimal parameters are output, the optimal parameters are obtained, and substituted into the decoupled echo state network structure constructed in step S6 - 1 to obtain the optimized decoupled echo state network structure.

[0224] Step 6 - 3, divide the comprehensive dynamic aging factor data set estimated in step 5 into a training set and a test set. Set the training set to 40%, 50%, 60%, 70%, 80% and 90% of the total duration of the comprehensive dynamic aging factor data set respectively, and the test set is the data set except the training set, a total of 6 groups.

[0225] Step 6 - 4, import the training set in step 6 - 3 into the improved structure network optimized in step 6 - 2 for training and prediction, and the prediction results are shown in Figure 7. Figures 7(a), (b), (c), (d), (e) and (f) are the results obtained by training and prediction under 40% training set, 50% training set, 60% training set, 70% training set, 80% training set and 90% training set respectively. The long - term prediction comprehensive aging index value fits the actual comprehensive aging index, and the results all show that the improved decoupled echo state network structure can achieve long - term accurate prediction of the fuel cell aging trend.

[0226] Step 6-5: Based on the comprehensive dynamic aging factor data set estimated in Step 5, take the end time of the total duration of the data set as the termination time of fuel cell failure. The time difference between the start time of prediction as the starting time and the termination time of fuel cell failure as the ending time is taken as the remaining service life. The specific calculation expression is as follows:

[0227] Life = Li(t start ) + Li(t end ) (43)

[0228] where Life is the remaining service life of the fuel cell; t start is the starting time of prediction as the starting moment; t end is the termination time of fuel cell failure; Li(t start ) is the remaining service life of the fuel cell at time t start ; Li(t end ) is the remaining service life of the fuel cell at time t end . Substitute the data set in Step 5 into formula (43) to calculate the actual remaining service life; substitute the long-term prediction data set in Step 6-5 into formula (43) to calculate the estimated remaining service life, and compare the actual and estimated remaining service lives. The comparison results are as Figure 8 shown.

[0229] From the experimental results in Figure 7, Figure 7(a), (b), (c), (d), (e) and (f) are the results obtained by training and prediction under 40% training set, 50% training set, 60% training set, 70% training set, 80% training set and 90% training set respectively. The long-term prediction values of the comprehensive aging index fit the actual comprehensive aging index. The results all show that the improved decoupled echo state network structure can achieve long-term accurate prediction of the fuel cell aging trend. The method proposed in the present invention can achieve long-term accurate prediction of the degradation performance of proton exchange membrane fuel cells under dynamic conditions. From Figure 8 the experimental results, it can be obtained that the method proposed in the present invention can accurately estimate the remaining service life of proton exchange membrane fuel cells.

Claims

1. A method for predicting the lifespan of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor, characterized in that: Specifically, it includes the following steps: Step 1, draw the electrochemical impedance spectroscopy diagrams at different times and the polarization curves at different times; Step 2, construct an aging equivalent circuit model of the proton exchange membrane fuel cell, and extract a data set of dynamic factor S1 based on the equivalent circuit model; Step 3, establish a semi-empirical equation, associate it with the equivalent circuit model to obtain an associated model, and extract a data set of dynamic factor S2 based on the associated model; Step 4, perform normalization processing on the dynamic factors S1 and S2, and obtain a data set of comprehensive dynamic aging factor S in the way of weighted superposition; Step 5, estimate the voltage and aging state of the fuel cell in real time to obtain a prediction data set of the fuel cell aging state; Step 6, realize the estimation of the remaining service life of the proton exchange membrane fuel cell.

2. The method for predicting the service life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor according to claim 1, wherein: The specific process of the said Step 1 is as follows: Step 1.1, obtain the impedance-frequency data set of the proton exchange membrane fuel cell stack through the impedance test of the proton exchange membrane fuel cell stack; Step 1.2, use the sliding average filtering algorithm to filter the data in the data set obtained in Step 1.1 to obtain the filtered values of the real part impedance and the imaginary part impedance after preprocessing; Step 1.3, import the impedance data set after preprocessing in Step 1.2 and the frequency data set in Step 1.1 in MATLAB, and draw the electrochemical impedance spectroscopy diagrams of different single fuel cells, at different times, and under different currents; Step 1.4, respectively test and obtain the data set of the fuel cell stack voltage, current, and working temperature through a voltage sensor, a current sensor, and a temperature sensor; Step 1.5, use the sliding average filtering algorithm to filter the data in the data set obtained in Step 1.4 to obtain the filtered values of the fuel cell stack voltage, current, and working temperature after preprocessing; Step 1.6, import the fuel cell stack voltage and current data sets after preprocessing in Step 1.5 and the time data set in Step 1.4 in MATLAB, and draw the polarization curves at different times.

3. The proton exchange membrane fuel cell life prediction method based on a comprehensive dynamic aging factor according to claim 2, characterized in that: The specific process of the said Step 2 is as follows: Step 2.1, connect the inductance element L Ohm and the resistance element R Ohm in series to represent the impedance Z Ohm of the ohmic polarization loss inside the fuel cell. Connect the capacitance element C Ca and the resistance element R Ca in parallel to represent the impedance Z Ca of the cathode polarization loss inside the fuel cell. Connect the impedance Z Ca in series with the resistance element R Con , and then connect it in parallel with the constant phase capacitance element CPE Con to represent the impedance Z Con of the concentration polarization loss inside the fuel cell; Step 2.2, connect the impedance Z Ohm constructed in Step 2.1 and the impedance Z Con in series successively to obtain an equivalent circuit model representing the aging of the proton exchange membrane fuel cell; Step 2.3, design the equivalent circuit model constructed in Step 2.2 in ZView software, and import the impedance-frequency data set of the proton exchange membrane fuel cell stack after filtering in Step 1.2; Step 2.4, on the basis of the equivalent circuit model designed in Step 2.3, respectively fit the electrochemical impedance spectroscopy diagrams of different single fuel cells, at different times, and under different currents drawn in Step 1.3, and extract the characteristic parameters of the equivalent circuit model of different single fuel cells, at different times, and under different currents; Step 2.5, at different times t, set the coupling superposition weight J of the constant current and the characteristic parameters of each equivalent circuit at different times, and calculate the dynamic aging factor S with equal current and different times t Data set; Step 2.6, at the a-th different current, set the coupling superposition weight of each characteristic parameter of the equivalent circuit under different current conditions to E, and calculate the dynamic aging factor S with different working currents a Data set; Step 2.7: The dynamic aging factors S of different working currents extracted in Step 2.6 a Through fusion calculation, a final data set of dynamic factor S1 based on the equivalent circuit model is obtained.

4. The method for predicting the lifespan of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor according to claim 3, wherein: The specific process of the said Step 2.4 is as follows: Step 2.4.1, under the conditions of equal current and consistent time, the number of single fuel cells that make up the proton exchange membrane fuel cell stack is N cell , use ZView software to fit the electrochemical impedance spectra of N cell single fuel cells in turn. According to the data feedback by ZView software, directly extract the characteristic parameters of the equivalent circuit model of single fuel cells with equal constant current and different times data sets; Where i represents the i-th single fuel cell, and the value range is i=1, 2, ..., N cell ; L Ohm-i and R Ohm-i represent the ohmic inductance data set and the ohmic resistance data set of the i-th single fuel cell respectively; C Ca-i and R Ca-i Respectively represent the polarization capacitance data set and polarization resistance data set of the i-th single fuel cell; R Con-i represents the concentration polarization resistance data set of the i-th single fuel cell; Q i and Respectively represent the double layer capacitance data set and the exponential parameter data set of the constant phase capacitance element of the i-th single fuel cell; [] T Represents the matrix full-time symbol; The fuel cell stack is obtained by connecting N cell single fuel cells in series, and calculate the characteristic parameters of the equivalent circuit model of the fuel cell stack data set: Among them, represents the summation symbol for the characteristic parameters of the equivalent circuit model of the i = 1, 2, …, N cell single fuel cells; L Ohm and R Ohm respectively represent the ohmic inductance and ohmic resistance of the fuel cell stack; C Ca and R Ca respectively represent the polarization capacitance and polarization resistance of the fuel cell stack; R Con represents the concentration polarization resistance of the fuel cell stack; Q and respectively represent the double-layer capacitance and exponential parameter of the constant phase capacitance element of the fuel cell stack; Step 2.4.2, under the condition of equal current and different times, mark the characteristic parameters of the equivalent circuit model of the fuel cell stack at different moments as F t , where t represents different moments, and its value range is t = 1, 2, …, T. According to the electrochemical impedance spectroscopy diagrams of the fuel cell at different moments drawn in Step 1.3, use the ZView software to fit the electrochemical impedance spectroscopy in turn, and extract the characteristic parameter F of the equivalent circuit model of the fuel cell with different times under a constant current t , F t is consistent with the expansion of the characteristic parameter F in formula (1), and F t represents F at different times; Step 2.4.3, under the condition of different currents, mark the characteristic parameters of the equivalent circuit model of the fuel cell stack with different current values as F a , where a represents the number of different current values, and the value range is a = 1, 2,..., A. According to the fuel cell electrochemical impedance spectrograms at different times drawn in Step 1.3, use ZView software to fit the fuel cell electrochemical impedance spectra under different current conditions in turn, and extract the characteristic parameter F of the equivalent circuit model of the fuel cell operating at different current values a , F a is consistent with the expansion of the characteristic parameter F in formula (1), and F a represents F under different currents.

5. The method for predicting the lifespan of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor according to claim 4, wherein: In the said Step 2.5, the coupling superposition weight J of each characteristic parameter of the equivalent circuit under a constant current and at different times is calculated by the following formula (2): Among them, w1, w2, …, w7 are the main diagonal elements in the coupling superposition weight J, representing the weight factors of the proportions of the seven characteristic parameters in the equivalent circuit characteristic parameter F at different times; the symbol "-" represents the normalization process; in formula (2) t the seven characteristic parameter proportion weight factors in; the symbol "-" represents the normalization process; in formula (2) the symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (3) is as follows: Where v represents the characteristic parameter F t The number of sorts in F is v = 1, 2, ..., 7; t-v (t) and F t-v (t-1) respectively represents the characteristic parameter F at time t and (t-1). t The vth characteristic parameter in; w v Represents the weight factor corresponding to the vth characteristic parameter in the coupling superposition weight J.

6. The proton exchange membrane fuel cell life prediction method based on the comprehensive dynamic aging factor according to claim 5, characterized in that: In the said Step 2.6, the coupling superposition weight E of each characteristic parameter of the equivalent circuit under different current conditions is calculated by the following formula (5): Among them, e1, e2, …, e7 are the main diagonal elements in the coupled superposition weight E, representing the weight factors of the proportions of the seven characteristic parameters in the equivalent circuit characteristic parameter F under different currents; the symbol "-" represents normalization processing; in formula (5) a the symbol represents the proportion of the seven characteristic parameters; the symbol "-" represents normalization processing; in formula (5) the symbol represents a specific decoupling calculation method, and the specific algorithm combined with formula (6) is as follows: Among them, m represents the sorting number in feature parameter F a , and the value range of m is m = 1, 2,..., 7; F a-m (a) and F a-m (a - 1) respectively represent the m-th feature parameter in feature parameter F a under the current values of a and (a - 1); w m represents the weight factor corresponding to the m-th feature parameter in the coupling superposition weight E.

7. The proton exchange membrane fuel cell life prediction method based on the comprehensive dynamic aging factor according to claim 6, characterized in that: The specific process of step 2.7 is as follows: the dynamic aging factors S corresponding to different working currents extracted in step 2.6 a are fused and calculated through formulas (8) and (9) to obtain the final data set of dynamic factor S1 based on the equivalent circuit model: Among them, j1, j2, …, jA are the main diagonal elements in the coupling superposition weight M, representing the weight factors of the proportions of A dynamic factors in the dynamic factor S a in; the symbol "-" represents normalization processing; the ▽ symbol in formula (8) represents a decoupled calculation method, and the specific algorithm combined with formula (9) is as follows: Among them, S a (t) and S a (t - 1) respectively represent the dynamic factors at the t-th and (t - 1)-th moments under the a-th current value; j a represents the weight factor in the coupling superposition weight M corresponding to the dynamic factor under the a-th current value.

8. The proton exchange membrane fuel cell life prediction method based on the comprehensive dynamic aging factor according to claim 7, characterized in that: The specific process of the said Step 3 is as follows: Step 3.1, construct a semi-empirical model of the proton exchange membrane fuel cell, as shown in the following formula (11): U stack = N cell ·(E OCV - U Ohm - U Con - U Ca ) (11) Among them, U stack is the stack voltage, N cell is the number of single fuel cells in the stack, E OCV is the open-circuit voltage at the rated pressure and rated temperature, U Ohm is the ohmic polarization voltage, U Con is the concentration polarization voltage, U Ca is the polarization voltage; In the semi-empirical model, the expression of the ohmic polarization voltage obtained from Ohm's law is: U Ohm = R eq · i stack (12) Among them, R eq is the total resistance, and i stack is the load current density; In the semi-empirical model, the concentration polarization voltage U Con is expressed as follows: Among them, λ is the charge transfer coefficient, i limit is the limiting current density, R is the ideal gas constant, T work is the stack temperature, with a value of T0 °C, n e is the number of electrons participating in the electrochemical reaction, F ara is the Faraday constant; In the semi-empirical model, the polarization voltage U Ca has the following expression: where, i change is the varying current density, and i leak is the leakage current density. In Equations (13) and (14), λ is a constant. In Equation (14), i leak is much smaller than i stack , and i leak can be ignored; finally, Equation (11) is simplified to: Among them, A ta is the Tafel constant, and B con is the concentration constant; Step 3.2: Combine the semi-empirical model of the fuel cell and the equivalent circuit model to obtain a macro-micro combined correlation model of the proton exchange membrane fuel cell, as shown in the following formula (16): U stack = N cell ·(E OCV - M Ohm - M Ca - M Con ) (16) Among them, M Ohm , M Ca and M Con respectively represent ohmic polarization loss, activation polarization loss, and concentration polarization loss in the model of the associated equivalent circuit characteristic parameters and the semi-empirical equation, and are specifically expressed as: Among them, S area is the cross-sectional area through which the operating current passes, with the unit of A / cm 2 ; |·| is the absolute value symbol, i.e., modulus calculation; Q Ohm , Q Ca and Q Con are the correlation coefficients of ohmic polarization loss, activation polarization loss, and concentration polarization loss in the associated model, respectively; Z Ohm , Z Ca and Z Con are the ohmic polarization impedance, activation polarization impedance, and concentration polarization impedance in the equivalent circuit model, respectively. Step 3.3, based on the correlation model, use the non-linear least squares method to fit the polarization curves at different times, and capture the characteristic parameters H = [E OCV , Q Ca , Q Ohm , Q Con of the correlation model at different times T ; Step 3.4: At different times, set the coupling superposition weight of each characteristic coefficient of the correlation model at different times to D, and calculate the correlation coefficient S2 of the correlation model through the characteristic parameter H of the correlation model: where d1, d2, d3, and d4 are the main diagonal elements in the coupling superposition weight D, representing the weight factors of the proportion of characteristic parameters in the characteristic parameters of the correlation model; the symbol "-" represents normalization processing; the X symbol in formula (18) represents a decoupled calculation method, and the specific algorithm combined with formula (19) is: Among them, \(l\) represents the \(l\)-th element in the characteristic parameter \(H\) of the correlation model, and the value range is \(l = 1, 2, 3, 4\); \(d\) l represents the weight factor corresponding to the \(l\)-th correlation model coefficient in the coupled superposition weight \(H\); \(H\) l \((t)\) and \(H\) l (t - 1) respectively represent the \(l\)-th characteristic parameter in the characteristic parameter \(H\) of the correlation model at times \(t\) and \((t - 1)\).

9. The method for predicting the service life of a proton exchange membrane fuel cell based on a comprehensive dynamic aging factor according to claim 8, characterized in that: The specific process of the said Step 4 is as follows: Step 4.1: Perform conventional normalization processing on the dynamic factor S1 data group and the dynamic factor S2 data group, and map the original data group to the dimension interval of [0, 1]; Step 4.2: At different times, set the coupling superposition weight of the dynamic factor based on the equivalent circuit model and the dynamic factor based on the correlation model at different times to C, and perform weight superposition on the dynamic factor S1 data group and the dynamic factor S2 data group after the normalization processing in Step 4.1 to calculate the comprehensive dynamic aging factor S data group: Among them, c1 and c2 are the main diagonal elements in the coupling superposition weight C, representing the weight factors of the dynamic factor based on the equivalent circuit model and the proportion of the dynamic factor based on the correlation model; the symbol "-" represents the normalization process; in formula (21) The symbol represents a decoupled calculation method, and the specific algorithm combined with formula (22) is as follows: Among them, k represents the sequence of dynamic factors, and the value range is k = 1, 2; c k represents the weight factor corresponding to the k-th dynamic factor in the coupling superposition weight C; S k (t) and S k (t - 1) respectively represent the k-th dynamic factor at times t and (t - 1). When k = 1, S k represents the dynamic factor based on the equivalent circuit model; when k = 2, S k represents the dynamic factor based on the correlation model.

10. The proton exchange membrane fuel cell life prediction method based on the comprehensive dynamic aging factor according to claim 9, characterized in that: The specific process of the said Step 5 is as follows: Step 5.1: Based on the extracted comprehensive dynamic aging factor S data group, adopt a discrete nonlinear system and the associated particle filter algorithm. Take the comprehensive dynamic aging factor as the state variable of the system state equation, and according to formula (24), complete the state estimation of the system parameter vector T including the equivalent circuit model and the joint model for 1 time step: T i = f i (T i-1 , ω i-1 ) (24) where, i is the time step; T is the system parameter vector including the equivalent circuit model and the joint model; f(·) is the degradation model describing the state change; ω i is the system noise; Step 5.2: Take the stack voltage data of the proton exchange membrane fuel cell as the observation variable in the observation equation, and the current and operating temperature data of the proton exchange membrane fuel cell as the input variables in the observation equation. Use the system state parameters estimated in Step 5.1, and according to formula (25) and formula (15), complete the observation estimation of the stack voltage of the proton exchange membrane fuel cell for 1 time step: z t = h i ((T t-1 , x t ), v t ) (25) where z i is the observation state at time i, T i-1 is the system state at time (i - 1), x t represents the input load current and operating temperature that can be planned in advance, v t represents the observation noise, and h(·) is the observation equation; Step 5.3: Compare the actual value of the stack voltage data group with the estimated value of the stack voltage data group estimated in Step 5.2, and realize the correction of the nonlinear system through the adjustment of the importance weight. Calculate the weight according to the observation likelihood function, and the likelihood function is usually defined in a Gaussian form: where r is the number of particles, and the value range is r = 1, 2,..., R; T is the matrix transpose symbol; V is the observation noise covariance matrix, representing the uncertainty of the observed value; exp represents the exponential operation with the natural constant e as the base; Step 5.4: Repeat Step 5.1 to Step 5.3 until the estimation of the stack voltage data group is completed; Step 5.5: End the estimation, and output the stack voltage data group of the observed parameters estimated in real time and the comprehensive dynamic aging factor data group of the output state parameters estimated in real time.

Images