A fuel cell fault diagnosis method based on double stack consistency constraints

Abstract

The application discloses a fuel cell fault diagnosis method based on double-stack consistency constraints and belongs to the technical field of fuel cells. The method comprises the following steps: S1, collecting original data in the operation process of double stacks of a fuel cell, performing pretreatment, and generating samples; S2, constructing a double-stack information mapping table, calculating double-stack difference characteristics, evaluating the state representation ability of the double-stack difference characteristics, and determining reserved double-stack difference characteristics; S3, constructing a global bias baseline, and performing bias correction on the reserved double-stack difference characteristics based on the constructed global bias baseline; S4, dividing the bias-corrected double-stack difference characteristics into working condition sensitive characteristics and non-working condition sensitive characteristics, and generating an adaptive threshold; and S5, acquiring and processing double-stack operation data of a fuel cell to be diagnosed, and performing fault determination based on the adaptive threshold. The application can realize fault diagnosis of the fuel cell under complex working condition conditions and has high accuracy.

Description

Technical Field

[0001] This invention relates to the field of fuel cell technology, and more specifically to a fuel cell fault diagnosis method based on dual-stack consistency constraints. Background Technology

[0002] Fuel cells have promising applications in vehicle power, distributed energy supply, and energy storage systems due to their cleanliness and high efficiency. However, in actual operation, fuel cells are susceptible to factors such as water management imbalance, abnormal gas supply, and thermal management fluctuations, which can lead to performance degradation or even failure.

[0003] Existing diagnostic methods mainly include data-driven diagnostic frameworks based on black-box models such as neural networks, and diagnostic frameworks based on rule-based thresholds and empirical criteria. While the former possesses strong pattern recognition capabilities, its interpretability is relatively limited, and it is quite sensitive to data quality and sample distribution. The latter, on the other hand, often relies heavily on empirical rules in threshold setting and boundary judgment, making it difficult to adapt to the diagnostic needs under complex and variable operating conditions. Furthermore, current research on fuel cell fault diagnosis largely focuses on the operation of individual cells or single-stack systems, with relatively insufficient research on dual-stack and multi-stack systems. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention proposes a fuel cell fault diagnosis method based on dual-stack consistency constraints.

[0005] The technical solution adopted in this invention is: A method for fault diagnosis of fuel cells based on dual-stack consistency constraints includes the following steps: S1. Collect raw data during the operation of the dual-stack fuel cell, preprocess the raw data, and generate samples; S2. Based on the samples generated in S1, construct a dual-pile information mapping table, calculate the dual-pile difference features, evaluate the state representation ability of the dual-pile difference features, and determine the retained dual-pile difference features. S3. Based on the retained dual-stack difference features determined in S2, a global bias baseline is constructed using samples under normal operating conditions, and the retained dual-stack difference features are biased and corrected based on the constructed global bias baseline to obtain the corrected dual-stack difference features. S4. Conduct a working condition sensitivity analysis on the corrected dual-stack difference characteristics obtained in S3, divide them into working condition sensitive features and non-working condition sensitive features, and construct dynamic threshold boundaries and global static boundaries respectively to form an adaptive threshold boundary library. S5. Acquire and process the dual-stack operation data of the fuel cell to be diagnosed, and determine the fault based on the adaptive threshold boundary library generated in S4.

[0006] The beneficial technical effects of this invention are: This invention is based on the principle that measurement points in a fuel cell should maintain basic consistency under the same operating conditions. Specifically, in a dual-stack fuel cell system, corresponding measurement points should have a stable correspondence under the same or comparable operating conditions. It constructs dual-stack difference features and further establishes a global bias baseline to correct the original difference features, thereby decoupling the inherent inconsistencies between the two stacks from actual fault disturbances. By analyzing the relationship between the corrected difference features and corresponding operating conditions, the difference features are divided into condition-sensitive and non-sensitive features. Dynamic threshold boundaries and global static boundaries that change with operating conditions are constructed respectively, achieving adaptive threshold generation that continuously changes with operating conditions. The candidate parameter combinations are evaluated using a validation sample set, reducing the subjective influence of manual experience-based settings. This invention can achieve fault diagnosis of fuel cells under complex operating conditions with high accuracy.

[0007] Specifically, the present invention has the following advantages: (1) High anti-interference capability: Through the construction of dual-stack difference features and global bias correction, the sensor zero drift and system static inconsistency are effectively eliminated, and the signal-to-noise ratio of fault signals is improved; (2) Adaptive capability of operating conditions: By utilizing the dynamic threshold boundary based on the operating condition compartment, the false alarm problem caused by baseline drift of fuel cells under wide operating conditions is solved; (3) Robustness of nonparametric statistics: The boundary is constructed by using nonparametric statistics such as quantiles, which reduces the dependence on the assumption of the characteristic distribution form and improves the applicability under complex nonlinear conditions. Attached Figure Description

[0008] Figure 1 This is a flowchart illustrating a fuel cell fault diagnosis method based on dual-stack consistency constraints according to the present invention. Figure 2 The results of sensitivity analysis ranking of different differential features in a specific application example of the present invention; Figure 3 This is a schematic diagram illustrating the construction of the dynamic threshold boundary of the working condition sensitive feature corr_diff4 in a specific application example of the present invention; Figure 4 This is a schematic diagram of window-level diagnostic result analysis in a specific application example of the present invention. Detailed Implementation

[0009] This invention is based on the difference characteristics of corresponding measurement points in a dual-stack fuel cell system. It uses a sliding window method for sample division. First, a global bias baseline is constructed using normal training samples, and the original difference characteristics are corrected. Then, the correlation between each corrected difference characteristic and the operating condition variables is analyzed to identify operating condition-sensitive and non-operating condition-sensitive characteristics, and dynamic threshold boundaries and global static boundaries that change with the operating conditions are established respectively. On this basis, the corresponding boundaries are adaptively called for the sample to be diagnosed, diagnostic features are generated, and state determination is performed, thereby realizing fuel cell fault diagnosis under complex operating conditions with high accuracy.

[0010] The following is a more detailed explanation in conjunction with the accompanying drawings.

[0011] like Figure 1 As shown, a fuel cell fault diagnosis method based on dual-stack consistency constraints includes the following steps: S1. Collect raw data during the operation of the dual-stack fuel cell, preprocess the raw data, and generate samples.

[0012] S2. Based on the samples generated in S1, construct a dual-stack information mapping table, calculate the dual-stack difference features, evaluate the state representation ability of the dual-stack difference features, and determine the retained dual-stack difference features.

[0013] S3. Based on the retained dual-stack difference features determined in S2, a global bias baseline is constructed using samples under normal conditions, and the retained dual-stack difference features are biased and corrected based on the constructed global bias baseline to obtain the corrected dual-stack difference features.

[0014] S4. Perform a working condition sensitivity analysis on the corrected dual-stack difference features obtained in S3, divide them into working condition sensitive features and non-working condition sensitive features, and construct dynamic threshold boundaries and global static boundaries respectively to form an adaptive threshold boundary library.

[0015] S5. Acquire and process the dual-stack operation data of the fuel cell to be diagnosed, and determine the fault based on the adaptive threshold boundary library generated in S4.

[0016] The above S1 includes the following steps: S11, Data Acquisition; Raw time-series data is collected during the operation of a dual-stack proton exchange membrane fuel cell system. This raw time-series data includes operating parameters at corresponding measurement points in both stacks, as well as operating condition variables reflecting changes in system operating conditions. The measurement points can be determined based on the measurement point layout of the dual-stack fuel cell system, fault diagnosis requirements, and the correlation between operating parameters and fault states. For example, the operating condition variable can preferably be current.

[0017] S12. Preprocess the raw time series data; The original time-series data undergoes data cleaning, data alignment and segmentation, and standardization to generate standardized input data required for subsequent research.

[0018] S13, Generate samples; For the standardized input data after the above processing, samples are generated using an overlapping sliding window method. Each sample includes a sequence of operating parameters for the corresponding measurement points of the two reactors within a sliding window, along with their corresponding operating condition variable sequences. For example, based on the approximately 20-second stabilization period set during fuel cell sampling, the window length is set to 20 and the step size to 5. Dynamic evolution continuity is achieved through data overlap between windows.

[0019] In this invention, the size and step size of the sliding window are selected based on the characteristics of data sampling. In actual operation, the size and step size of the sliding window can be adjusted according to the length of time it takes to stabilize after changes in working conditions.

[0020] S14. Sample partitioning; The samples are divided into training and validation sets in an 8:2 ratio. Normal training samples refer to samples assigned to the training set under normal operating conditions, and normal validation samples refer to samples assigned to the validation set under normal operating conditions. Fault training samples refer to samples assigned to the training set under fault conditions, and fault validation samples refer to samples assigned to the validation set under fault conditions.

[0021] The above S2 includes the following steps: S21. Construct a dual-heap information mapping table; Based on the samples generated in S13, a mapping relationship between corresponding sensor measurement points in the dual-stack fuel cell system is established. For the physical quantity to be analyzed, the corresponding measurement points in stack 1 and stack 2 are determined respectively. A unified dual-stack information mapping table, i.e., a measurement point mapping table, is constructed and used... Indicates the difference characteristics, where, Indicates the first Each dual-stacking pile corresponds to a pair of measurement points. The dual-stacking information mapping table mainly includes: Number, corresponding physical quantity, column corresponding to pile 1, column corresponding to pile 2.

[0022] S22. Calculate the difference characteristics between the two piles; Based on the dual-heap information mapping table constructed using S21, the difference between each mapping point is calculated, and a difference feature vector is constructed. For the first... There are corresponding measurement points, at time... The dual-stack difference characteristic is defined as: ; in, Indicates the first heap in heap 1 The measured values ​​at each measuring point at time t. This indicates the time at which the corresponding measurement point in pile 2 is located. The measured value, Indicates the first The dual-pillar difference characteristics at time t of each measurement point. In this invention, the first... The pair of measurement points corresponding to the double piles and the first point calculated from it. Each of the two pile differences corresponds one-to-one; therefore, in the following text... A unified representation is used to indicate the pair number of the measurement points corresponding to the two piles and the corresponding difference feature number of the two piles.

[0023] The above S3 includes the following steps: S31. Determine the global bias baseline; Dual-reactor systems have inherent biases due to factors such as manufacturing, installation, external environment, and operating time; therefore, this invention needs to determine the global bias baseline and calculate and correct the dual-reactor difference characteristics.

[0024] (1) Single-sample bias estimation; The first generated in step S13 The first sample The bias of the two-stack difference feature, namely the difference feature calculated in step S22, is defined as: ; in, Indicates the first The number of valid sample points in a sample. Indicates the first In the nth sample The two-stack difference features at time... The value of , This represents the bias corresponding to the k-th sample.

[0025] (2) Global baseline synthesis; The biases of all normal training samples are summed in a second order, and the global bias baseline of the i-th dual-stack differential feature is: ; in, Represents the total number of samples. Indicates the first Global bias baseline for dual-stack difference characteristics This represents the bias corresponding to the k-th sample.

[0026] S32. Calculate and correct the difference characteristics between the two stacks; Based on the global bias baseline determined in S31, the retained differential features obtained in S2 are subjected to bias correction, and the calculation formula is as follows: ; in, For the first In the nth sample The two-stack difference features at time... The value of , For the first Global bias baseline for dual-stack difference characteristics For the corrected first In the nth sample Two-stack difference characteristics. This formula is used to calculate and decouple static errors and dynamic anomalies in the original signal.

[0027] S41. Conduct a working condition sensitivity analysis; To further distinguish between normal fluctuations caused by changes in operating conditions and abnormal deviations caused by faults in the differential characteristics after bias correction, this step analyzes the relationship between the differential characteristics of each diagnostic corrected dual-stacking system and operating condition variables, identifying characteristics that are significantly affected by operating conditions. Specific relevant indicators are as follows: Calculate the differences between the two stacks and the operating condition variables. Correlation indicators: Calculate the Pearson linear correlation coefficient and the Spearman rank correlation coefficient to characterize the degree of linear correlation and monotonicity.

[0028] ; In the formula, Indicates the first The corrected dual-stack difference characteristics at time 1 The value of , This represents the mean of the corrected differential characteristics. Indicates time Operating condition variables, This represents the mean of the operating condition variables. The Pearson linear correlation coefficient; express Rank in the entire sample express Rank in the entire sample Indicates the first The mean of the ranks of each corrected difference feature. The mean of the rank of the operating condition variable. Indicates the first Spearman's rank correlation coefficient between the corrected differential characteristics and the operating condition variables.

[0029] Furthermore, to determine whether the corrected dual-stack difference characteristics undergo systematic drift with changes in operating condition variables, a sensitivity analysis of operating conditions was conducted based on normal training samples. Since the range of operating condition variables within a single normal training window is limited, and the overall statistics of all normal training windows are insufficient to reflect the distributional differences between different operating condition intervals, the operating condition variables corresponding to the normal training samples were divided into several operating condition sensitivity evaluation bins using a quantile binning method. .

[0030] The quantile binning method is as follows: multiple quantile points are determined according to a preset number of bins, and the operating condition variable values ​​corresponding to adjacent quantile points are used as the boundaries of the operating condition interval.

[0031] Extract the operating condition variable values ​​corresponding to the normal training samples and sort them in ascending order of value; then, according to the preset number of bins... Sure The corresponding quantiles are used; then, the operating condition variable values ​​corresponding to adjacent quantiles are used as the boundaries of the operating condition intervals to form several operating condition sensitivity evaluation bins. Each normal training sample is assigned to the corresponding bin according to its operating condition variable value. When the operating condition variable values ​​corresponding to adjacent quantiles are the same, or when the number of samples in a bin is too small to stably calculate the bin statistics, adjacent bins can be merged, or the preset number of bins can be adjusted and re-binded. The operating condition sensitivity evaluation bins are used to calculate the bin mean, bin standard deviation, and inter-group explained value of the corrected dual-stack difference characteristics within different operating condition ranges to determine whether the corrected dual-stack difference characteristics undergo systematic drift with changes in operating condition variables.

[0032] For each load condition sensitivity evaluation bin, the sample size, bin mean, and bin standard deviation of the i-th corrected difference feature within that bin are calculated. The changes in bin mean, bin standard deviation, and inter-group explanatory power are compared across different load condition sensitivity evaluation bins to determine whether the normal value variation of the corrected difference feature in different load condition intervals undergoes a systematic drift with changes in load conditions. The normal value variation is described by the sample size, bin mean, and bin standard deviation within each load condition sensitivity evaluation bin, as well as the population centrality and population standard deviation across all normal training samples. The specific calculation process is as follows: (1) Mean drift index: ; (2) Changes in fluctuation scale: ; (3) Intergroup explanatory power index: ; Where b represents the working condition compartment number; The mean shift index represents the i-th corrected difference feature, used to represent the... The maximum difference in the mean of each correction difference feature under different operating conditions and binning is the proportion of the overall fluctuation scale of the feature. This represents the bin mean of the i-th corrected difference feature within the b-th bin; This represents the maximum value of the mean across all operating conditions. This represents the minimum value of the mean across all operating conditions. For the first The overall standard deviation of each corrected difference feature across all normal training samples; It is a very small positive number, used to prevent the denominator from being zero.

[0033] For the first The fluctuation scale change index of each characteristic; Indicates the first The correction difference feature in the first Standard deviation within each operating condition compartment; This represents the maximum standard deviation within all operating conditions across all bins. This represents the minimum standard deviation within all operating conditions.

[0034] Indicates the inter-group explanatory power index. Indicates the corrected binary stack difference feature number; Indicates the total number of compartments under different operating conditions; This represents the total number of normal training samples. Indicates the first The number of normal training samples in each working condition bin; Indicates the first In the nth normal training sample The values ​​of the corrected dual-stack difference features; This represents the overall center position of the i-th corrected difference feature across all normal training samples.

[0035] Using the Pearson linear correlation coefficient, Spearman rank correlation coefficient, mean drift index, fluctuation scale change index, and between-group explanatory power index, the condition sensitivity evaluation vector for the i-th corrected difference feature is formed: ; The sensitivity thresholds include strong sensitivity thresholds and weak sensitivity thresholds. For Pearson linear correlation coefficient, Spearman rank correlation coefficient, mean drift index, fluctuation scale change index, and between-group explanatory power index, the values ​​of all corrected dual-stack difference features on the corresponding indexes are statistically analyzed and sorted according to their numerical values. The upper quartile of the sorted results is taken as the strong sensitivity threshold, and the median of the sorted results is taken as the weak sensitivity threshold. When any operating condition sensitivity index reaches the corresponding strong sensitivity threshold, or at least two operating condition sensitivity indices reach the corresponding weak sensitivity threshold, the corrected dual-stack difference feature is classified as an operating condition sensitive feature. Otherwise, it is classified as a non-operating condition sensitive feature. The strong and weak sensitivity thresholds are only used to separate operating condition sensitive features from non-operating condition sensitive features.

[0036] S42. Characteristic flow division and boundary type determination; Based on the operating condition sensitivity evaluation vector obtained from S41, the dual-stacking correction difference features are divided into operating condition sensitive features and non-operating condition sensitive features. For operating condition sensitive features, a dynamic threshold boundary that varies with the operating condition is constructed; for non-operating condition sensitive features, a global static boundary is constructed.

[0037] I. Construction of dynamic threshold boundaries for operating condition sensitive features.

[0038] For operating condition-sensitive features, their normal values ​​fluctuate with changes in operating condition variables. Therefore, we consider constructing a dynamic threshold boundary that changes with operating conditions based on normal training samples. First, we use the number of bins for the operating condition, the inner quantile parameter, and the outer quantile parameter as parameters to be optimized, and construct candidate parameter combinations. The specific steps are as follows: (1) Determine the number of bins; divide the normal training samples into several working condition bins according to the working condition variables, generate a candidate bin number set according to the number of normal training samples, and satisfy the minimum bin sample number constraint, as follows: ; in, ; This refers to the number of boxes. This represents the number of training windows for normal samples. The minimum number of samples allowed per bin. To minimize the number of boxes, This represents the maximum number of bins. To ensure statistical reliability, the number of valid samples within each bin for each operating condition should not be less than the minimum number of samples.

[0039] (2) Construction of feature sets within bins; For the i-th condition-sensitive feature and the b-th condition bin, the corrected double-stack difference feature values ​​of the condition variable values ​​falling into the bin of that condition are extracted from the normal training samples, forming a bin feature set. Based on the bin feature set, the median and interquartile range are calculated, which are used as the center position and fluctuation scale under that condition bin, respectively. The purpose of this binning step is mainly to establish dynamic threshold boundaries for features that have been identified as condition-sensitive under different conditions. Specifically: For the The first condition-sensitive characteristic, in the first Within each working condition bin, extract samples from normal training samples that meet the requirements. Feature set: ; in, Indicates the first The sensitive feature of the working condition is in the first The normal training set sample set within each working condition bin. Indicates the first In the nth normal training sample The corrected feature values ​​of each difference feature. Indicates the first The operating condition variable values ​​corresponding to each window Indicates the first Each working condition is divided into separate boxes. This represents the normal training sample set.

[0040] In each Within this range, the median and interquartile range of each feature are calculated, serving as the center position and the fluctuation scale, respectively, to describe the normal sample distribution within this operating condition range. The median is used to represent the baseline level of the feature within this operating condition range, while the interquartile range is used to represent the dispersion of the feature within this operating condition range.

[0041] (3) Construction of inner and outer boundaries; An inner and outer boundary are constructed, where the inner boundary represents the main fluctuation range of normal samples, and the outer boundary represents the boundary for determining significant anomalies. Let the inner and outer quantile parameters be: and And satisfy Then, the inner layer partition boundary of the sensitive feature of the i-th working condition under the b-th working condition is: ; in, Indicates the first The sensitive feature of the working condition is in the first The lower inner boundary of each working condition compartment Indicates the first The sensitive feature of the working condition is in the first The inner upper boundary of the compartment under each working condition. Indicates the first The sensitive feature of the working condition is in the first The set of feature values ​​of normal training samples within each working condition bin.

[0042] The outer boundary is: ; in, Indicates the first The sensitive feature of the working condition is in the first The outer lower boundary of the compartment under each working condition. Indicates the first The sensitive feature of the working condition is in the first The outer upper boundary of the compartment under each working condition.

[0043] Each operating condition bin corresponds to a set of local normal boundaries. When a sample to be diagnosed enters the diagnostic phase, the system calls the boundaries under the corresponding bin to make a judgment based on the operating condition bin to which its current operating condition variable value belongs.

[0044] (4) Dynamic threshold boundary interpolation; To address the threshold abrupt change at the boundary of adjacent operating condition bins, a linear interpolation algorithm is introduced. When the operating condition variable of the window to be diagnosed is located within an adjacent operating condition bin... and When interpolating between adjacent operating conditions, the boundary parameters corresponding to the sub-bins are interpolated. The specific formula is as follows: ; in, Indicates the first Individual working condition distribution center Indicates the first Individual working condition distribution center Indicates the first The operating condition variable values ​​for each window to be diagnosed. Indicates the first One feature is located in the center of the upper compartment. Boundary parameters at that location, Then it means the first One feature is located in the center of the upper compartment. The boundary parameters at the location include at least one of the inner lower boundary, inner upper boundary, outer lower boundary, and outer upper boundary.

[0045] (5) Boundary validity check; To ensure that the generated threshold boundaries have effective diagnostic width and numerical stability, the candidate dynamic threshold boundaries are subjected to validity checks. When a candidate boundary exhibits zero scale, overlapping inner and outer boundaries, boundary reversal, an outer boundary narrower than the inner boundary, or a boundary width lower than the minimum effective width, it indicates that the candidate boundary cannot stably represent the normal fluctuation range, and the system automatically discards the candidate parameter combination.

[0046] Wherein, the boundary width is the difference between the upper boundary and the lower boundary, and the effective boundary width should satisfy:

[0047] The minimum effective width can be determined by multiplying the overall fluctuation scale of features in the normal training samples by a preset width coefficient. The overall fluctuation scale can be one of the standard deviation, interquartile range, or range. The preset width coefficient is used to avoid zero-width boundaries or invalid boundaries caused by insufficient sample size, sensor resolution limitations, or excessively small fluctuations within bins, and is not used as the final fault diagnosis threshold itself.

[0048] This process avoids threshold boundary degradation caused by excessive sample concentration or insufficient binning samples, thereby improving the stability of subsequent fault determination and continuous deviation calculation.

[0049] (6) Evaluation of normal validation sample coverage; Its coverage was evaluated on normal validation samples. The number of normal samples was [number missing]. , No. The feature in the first The feature value on a normal verification window is Its operating condition variables are The corresponding inner and outer boundaries can be obtained by calling or interpolating based on the operating condition value.

[0050] Among them, the proportions of samples falling into the inner and outer layer boundaries in normal verification samples are as follows: ; And calculate the coverage deviation based on the difference between the actual coverage and the theoretical target coverage: ; in, This indicates the actual coverage rate of normal validation samples falling into the inner layer boundary. This indicates the actual coverage rate of normal validation samples falling into the outer boundary. To ensure a normal number of validation samples, For the first In the first normal verification sample One corrected dual-stack difference eigenvalue, Let K be the operating condition variable corresponding to the k-th normal verification sample. and These represent the i-th feature in the operating condition variables. The corresponding dynamic inner lower boundary and dynamic inner upper boundary are below. and These represent the i-th feature in the operating condition variables. The corresponding dynamic outer lower boundary and dynamic outer upper boundary are below. This indicates an indicator function that takes the value 1 if the condition within the parentheses is true, and 0 otherwise. The coverage rate of the inner layer theoretical target. For the coverage of the outer theoretical target, This represents coverage deviation.

[0051] Coverage evaluation is used to determine whether the candidate dynamic boundary can reasonably cover the normal verification samples. It is compared with a set threshold. If the coverage is too low, it means that the boundary is too narrow, which is prone to false alarms; if the coverage is too high, it means that the boundary is too wide, which may weaken the fault identification capability.

[0052] (7) Determination of target parameters for the dynamic threshold library; Based on the boundary validity check results, coverage deviation, and diagnostic performance evaluation results, a target threshold boundary parameter combination is determined from the candidate threshold boundary parameter combinations. A dynamic threshold boundary corresponding to the condition-sensitive features is then generated based on the target threshold boundary parameter combination. The candidate threshold boundary parameter combination consists of the number of candidate bins, candidate inner-layer quantile parameters, and candidate outer-layer quantile parameters. The diagnostic performance evaluation result is the diagnostic performance index calculated based on the known state labels of the validation set samples after the candidate threshold boundary parameter combination has undergone boundary invocation, feature-level judgment, and window-level aggregation processes.

[0053] The specific process is as follows: 1. Perform boundary validity checks according to step (5). Eliminate candidate solutions with insufficient binning sample size, inverted boundaries, outer boundaries that do not include inner boundaries, or boundary widths that are too small; 2. Following step (6), calculate the inner and outer layer coverage and coverage deviation of the normal validation sample; 3. Diagnostic performance is evaluated by combining the false alarm rate, false negative rate, and F1 score on the validation set. On the validation set, the candidate boundary is used for characteristic-level fault determination, and the false alarm rate is calculated. (Proportion of normal samples judged as faulty), false negative rate (The proportion of faulty samples judged as normal) and F1 score.

[0054] Finally, among the candidate parameter combinations that passed the boundary validity check, the parameter combination with coverage close to the theoretical target coverage, small coverage deviation, boundary width meeting the minimum effective width requirement, and superior diagnostic performance indicators was selected as the target threshold boundary parameter combination. Based on this target threshold boundary parameter combination, a dynamic threshold boundary corresponding to the condition-sensitive features was generated. In this method, the number of candidate bins and candidate quantile parameters only serve as the parameter search space; the final dynamic threshold boundary is jointly determined by normal training samples, normal validation samples, and diagnostic performance evaluation results.

[0055] II. Construction of global static boundaries for non-operating condition sensitive features; For features that are not sensitive to operating conditions, a global static boundary independent of operating conditions is constructed directly based on the center position, fluctuation scale, and quantile boundaries of normal training samples within the global scope. The specific steps are as follows: Let the first The set of non-operating condition sensitive features in the normal training samples is as follows: ; in, This represents the global feature set of the i-th non-condition-sensitive feature in the normal training samples. This represents the value of the i-th non-condition-sensitive feature in the k-th normal training sample. This represents the normal training sample set.

[0056] Compute on a set: (1) Central position : This is used to characterize the global baseline level of this non-operating condition sensitive feature in normal training samples. This represents the median function.

[0057] (2) Fluctuation scale : This is used to characterize the global dispersion of the non-operating condition sensitive feature in normal training samples, and is used for boundary validity evaluation and diagnostic deviation calculation. This represents the standard deviation function.

[0058] (3) Constructing inner and outer boundaries based on quantiles: Let the inner quantile parameters and outer quantiles be respectively: and And satisfy Then the first The global static inner layer quantile boundary for each non-operating condition sensitive feature is: ; The outer boundary is: ; in, and Let represent the lower boundary and upper boundary of the global static inner layer of the i-th non-operation condition sensitive feature, respectively; and Let $\frac{i}{i}$ represent the lower and upper boundaries of the global static outer layer for the $i$-th non-operational condition-sensitive feature, respectively. For candidate global static boundaries, normal validation samples are used to check boundary validity and evaluate coverage. Combined with the diagnostic performance evaluation results of the validation set, the target quantile parameter combination corresponding to the global static boundary is determined.

[0059] S43. Generate an adaptive threshold boundary library; The dynamic threshold boundaries corresponding to the working condition sensitive features obtained in step S42 and the global static boundaries corresponding to the non-working condition sensitive features are uniformly organized to form an adaptive threshold boundary library.

[0060] The adaptive threshold boundary library includes feature number, boundary type, boundary value, and applicable operating condition information. The boundary type is used to distinguish between dynamic threshold boundaries and global static boundaries; the boundary value includes at least one of the following: inner lower boundary, inner upper boundary, outer lower boundary, and outer upper boundary.

[0061] For dynamic threshold boundaries, the adaptive threshold boundary library can also record the number of bins used to generate the boundary, the inner quantile parameters, the outer quantile parameters, as well as the boundary validity check results and the normal validation sample coverage evaluation results. This information characterizes the generation basis and validity status of the dynamic threshold boundary, facilitating subsequent boundary calls, parameter tracing, and boundary updates. For global static boundaries, the adaptive threshold boundary library records the corresponding global boundary values ​​and boundary validity status.

[0062] The above S5 includes the following steps: S51, Adaptive Boundary Calling and Diagnostic Feature Generation; For the dual-stack operation data of the fuel cell to be diagnosed, the preprocessing and sliding window method in S1 is used to generate the sample to be diagnosed. The dual-stack difference characteristics of the sample to be diagnosed are calculated according to S2. The bias correction is performed based on the global bias baseline according to S3 to obtain the corrected dual-stack difference characteristics of the sample to be diagnosed. Based on the operating condition sensitivity analysis results and adaptive threshold boundary library obtained in S4, dynamic threshold boundaries are called for the operating condition sensitive features according to their corresponding operating condition variables, and global static boundaries are called for the non-operating condition sensitive features to generate window-level diagnostic features for fault determination.

[0063] S52, Feature-level determination; Based on the diagnostic feature vector generated by S51 and combined with the adaptive threshold boundary library, a primary discrete state determination is performed on the operating state of the fuel cell system, as follows: Let the diagnostic value of the i-th differential feature in the k-th diagnostic window be... Its corresponding inner dynamic threshold boundary is The outer dynamic threshold boundary is The inner boundary characterizes the normal fluctuation zone, while the outer boundary serves as the basis for determining whether a fault has occurred. The threshold boundary is either a dynamic threshold boundary obtained from S51 or a global static boundary.

[0064] when If so, the feature is determined to be in a normal state; when or If so, the feature is determined to be in a faulty state.

[0065] S53, Window-level aggregation rules; Based on feature-level determination, multidimensional differential features within the same sample are aggregated to form a diagnostic result from feature level to window level. Let there be P differential features in the k-th sample, of which the number of features identified as fault states is... The percentage of abnormal features is: ; in, This represents the proportion of outlier features in the k-th sample. Let the threshold for the number of outlier features be... The threshold for the proportion of abnormal features is ,when or If the condition is met, the k-th sample is classified as a faulty sample; otherwise, the k-th sample is classified as a normal sample. The threshold for the number of abnormal features... and the threshold for the proportion of abnormal features The parameters are used as window-level aggregation rules and determined through validation set evaluation. Ultimately, the window-level decision result is used as the basis for determining whether a fault has occurred.

[0066] S54. Validation set evaluation; To further illustrate the diagnostic performance evaluation method of the candidate threshold boundary parameter combinations in S42, and to determine the window-level aggregation rule parameters in S53, candidate parameter combinations are constructed, and the diagnostic effect of each candidate parameter combination is calculated based on the validation set samples. The candidate parameter combinations include the candidate threshold boundary parameters used in S42 to generate adaptive threshold boundaries, and the candidate aggregation rule parameters used in S53 for window-level aggregation determination; wherein, the threshold boundary parameters include at least one of dynamic threshold boundary parameters and global static boundary parameters, and the window-level aggregation rule parameters include at least one of anomaly feature quantity threshold and anomaly feature proportion threshold.

[0067] Specifically, each set of candidate parameter combinations is substituted into the diagnostic process from S51 to S53 to obtain the window-level prediction result corresponding to the validation set sample. The diagnostic performance index corresponding to the candidate parameter combination is calculated with reference to the known state label of the validation set sample. The diagnostic performance index includes at least one of accuracy, F1 score, false alarm rate and false negative rate.

[0068] Accuracy rate is used to represent the overall degree of correctness in judgments; The F1 score is used to comprehensively reflect the fault identification capability. False alarm rate is used to represent the proportion of normal samples that are mistakenly identified as faulty. The false negative rate is used to represent the proportion of faulty samples that are mistakenly identified as normal. Based on the diagnostic performance metrics, a target parameter combination is determined from the candidate parameter combinations. As a preferred approach, the accuracy and F1 score of the candidate parameter combinations on the validation set are maximized, while the false positive rate and false negative rate are minimized. Pareto non-dominated sorting is used to screen Pareto front solutions that are not dominated by other candidate solutions. When multiple Pareto front solutions exist, the target parameter combination is determined according to the priority order of lower false positive rate, lower false negative rate, and higher F1 score. Alternatively, a weighted approach can be used; specifically, for each set of candidate parameters, its diagnostic performance on the validation set is calculated.

[0069] The target parameter combination is used to determine the target threshold boundary parameter combination corresponding to the adaptive threshold boundary in S42, and the abnormal feature quantity threshold and abnormal feature proportion threshold corresponding to the window-level aggregation rule in S53.

[0070] As a further design of the present invention, the method also includes a stability verification step for the bias result: Based on the window bias results extracted from S31, the discrete measure of each differential feature channel in the full sample space is calculated to assess its reliability as a diagnostic baseline, i.e., to evaluate... To assess reliability, a global static bias baseline is used to characterize the long-term, stable, inherent differences between the corresponding measurement points of the two reactors under normal conditions. The specific steps are as follows: S311. Construct stability statistics; For the window bias results of each normal training sample obtained in S31, the range, absolute median and standard deviation are calculated on all normal training samples as statistics for evaluating the stability of the dual-stack difference characteristics under normal conditions.

[0071] S312, Calculation of coefficient of variation correction; Calculate the corrected coefficient of variation This is used to characterize the relative dispersion of the window bias result within the normal training sample space: ; in, Indicates the first Global bias baseline for dual-stack difference characteristics Indicates its standard deviation, A very small positive number is introduced to prevent the denominator from being zero.

[0072] Step 313: Automatically assign stability index weights; To reduce the subjectivity of manually assigned weights, the stability metrics were ranked on both normal training and normal validation samples, as follows: (1) Sort each indicator separately. For each type of indicator, the smaller the indicator value, the higher the ranking. (2) Calculate the ranking on the normal training set and the normal validation set respectively; (3) For each indicator, compare whether its ranking in the training set and the validation set is consistent. The smaller the ranking difference, the more stable and reliable the indicator is, and the higher its weight is. The specific calculation method is as follows: ; (4) Generate weights based on sort consistency; ; (5) Final score; Calculate the overall ranking score based on the statistical results of the training set: ; in, It is the first The dual-stack difference feature in the first Ranking under each indicator This is the weight of the indicator.

[0073] Step S314: Stability evaluation; A lower overall score indicates greater overall stability and a higher ranking for the differential feature. The stability assessment results are used to determine the reliability of the global bias baseline for each dual-reactor differential feature. For differential features with high stability, their inherent bias fluctuations under normal conditions are small, making bias correction using the global bias baseline highly reliable. For differential features with low stability, their bias fluctuations under normal conditions are large and may be affected by changes in operating conditions, random disturbances, or other dynamic factors. Therefore, further analysis of operating condition sensitivity is needed to determine whether the fluctuations are systematically related to operating condition variables, and based on this, to decide whether to use a dynamic threshold boundary or a global static boundary.

[0074] Therefore, stability assessment is mainly used to evaluate the reliability of the global bias baseline, while operational sensitivity analysis is used to determine the boundary modeling method for the difference characteristics after bias correction. Both provide a basis for selecting the subsequent threshold boundary modeling method. Based on the stability analysis of the bias results in S31, the reliability of the global bias baseline for the difference characteristics of each diagnostic dual-stack can be evaluated.

[0075] Furthermore, S51 above includes the following steps: S511. Based on the feature condition sensitivity analysis results obtained from S4 and the constructed baseline library, different boundary invocation strategies are executed for different types of differential features. For condition-sensitive features, the dynamic threshold boundary library is invoked according to the corresponding condition variables; for non-condition-sensitive differential features, the global static boundary library is directly invoked.

[0076] S512, Generate window-level diagnostic features and statistics; This step is performed separately in the offline modeling phase and the online diagnostic phase, as detailed below: (1) In the offline modeling stage, samples are generated according to the preprocessing and sliding window partitioning method of S1 for normal training samples, normal verification samples and fault verification samples used for diagnostic performance evaluation. The double-stack difference features in each sample are calculated according to S2, and the bias correction is performed according to the global bias baseline obtained in S3 to obtain the corrected double-stack difference feature sequence.

[0077] Within each sliding window, window-level statistics are extracted for each corrected dual-stack difference feature. These window-level statistics include mean, standard deviation, median, maximum, minimum, quantile, range, linear change slope, and deviation from the corresponding threshold boundary. These statistics are used to construct a window-level diagnostic feature matrix for subsequent feature-level and window-level aggregation determinations.

[0078] (2) For the real-time collected fuel cell dual-stack operation data, first form a sample to be diagnosed according to the same data processing and sliding window generation method in S1; then calculate the dual-stack difference characteristics according to S2, and perform bias correction according to S3; finally, call the adaptive threshold boundary library generated in S4 according to the operating condition variables corresponding to the current sample to generate the diagnostic features and statistics of the sample to be diagnosed, which are used for subsequent feature-level judgment and window-level aggregation judgment.

[0079] S513, Adaptive Boundary Call; Based on the above-described traffic splitting mechanism, the adaptive boundary is generated as follows: ; ; The dynamic threshold boundary library established in step 42, The global static boundary library is established for step 42.

[0080] The invention will be further explained below with reference to specific application examples.

[0081] S1. Data acquisition and preprocessing; S11, Data Acquisition; The collected data includes, but is not limited to, the dual-stack information in Table 1.

[0082] S12. Preprocess the raw data; The raw data is cleaned, aligned, and segmented, and then standardized to generate standardized input data required for subsequent research.

[0083] Step S13: Generate samples; The above data is processed using a sliding window with a window length of 20 and a step size of 5.

[0084] Step S14: Sample partitioning; The samples are divided into training and validation sets in an 8:2 ratio. Normal training samples are those assigned to the training set under normal conditions, normal validation samples are those assigned to the validation set under normal conditions, fault training samples are those assigned to the training set under fault conditions, and fault validation samples are those assigned to the validation set under fault conditions.

[0085] S2, Construction of dual-stack differential features; S21. Generate a dual-heap information mapping table; The generated specific information mapping table is shown in Table 1: Table 1

[0086] S22, Calculation of difference characteristics; The difference feature is calculated based on the dual-stack information mapping table in S21. The specific calculation formula is as follows: ; S3, Global bias baseline construction and consistency correction; S31. Determine the global bias baseline; This embodiment calculates the window bias for the differential features within each window based on normal training samples, and summarizes the bias results for all normal training samples to obtain the global bias baseline for each differential feature. An example of the calculation results is shown in Table 2 below, which is the global bias baseline table.

[0087] Table 2

[0088] Stability verification of bias results: To evaluate the reliability of the bias results of each differential feature as a baseline, stability indices, including standard deviation, range, coefficient of variation, and median absolute deviation, were calculated on normal training and validation samples. A smaller stability index indicates less fluctuation in the bias of that differential feature, and thus higher reliability as a global bias baseline.

[0089] To reduce the subjectivity of manually assigning indicator weights, this embodiment further employs automatic allocation of stability indicator weights based on training-validation set ranking consistency. Specifically, for each stability indicator, the differential features of each binary set are sorted in ascending order on both the normal training and validation sets. A smaller indicator value indicates less window bias fluctuation corresponding to that feature, and a higher ranking. For the... For each stability metric, the average ranking difference between the metric and the normal training and validation sets is calculated. The smaller the ranking difference, the more consistent the evaluation results of the metric are across different sample sets, and the higher the reliability of the stability evaluation. Further, the ranking consistency of the metric is calculated based on the average ranking difference, and the ranking consistency of each stability metric is normalized to obtain the automatic weights for each stability metric.

[0090] After obtaining the weights of each stability metric, a comprehensive stability score is calculated for each dual-hump differential feature based on its metric ranking on the normal training set. The comprehensive stability score is the weighted sum of the ranking of each stability metric and its corresponding automatic weight. The smaller the score, the higher the comprehensive ranking of the dual-hump differential feature under multiple stability metrics, and the more stable and reliable its global bias baseline is. This comprehensive stability ranking result is only used to evaluate the reliability of the global bias baseline of each dual-hump differential feature and is not used as the basis for individually removing features. In this embodiment, the comprehensive stability ranking result is shown in Table 3 below.

[0091] Table 3

[0092] The ranking is only used to evaluate the reliability of the bias baseline and is not used as a sole basis for exclusion.

[0093] S32. Calculate the correction difference characteristics; Based on the global bias baseline obtained from S31, bias subtraction is performed on the retained original difference features to obtain the corrected difference features: ; S4, Adaptive threshold generation; The threshold boundary library is used to characterize the allowable fluctuation range of various differential features under different operating conditions in normal state. Therefore, its construction is mainly based on normal training samples, and boundary validity checks, coverage evaluations, and coverage deviation calculations are performed using normal validation samples to screen reasonable candidate threshold boundary parameter combinations. Fault samples are not included in the construction of the normal baseline library to avoid fault disturbances contaminating the normal boundaries.

[0094] S41. Operating condition sensitivity analysis; After obtaining the correction difference characteristics, the relationship between each correction difference characteristic and the operating condition variable is analyzed using current as the operating condition variable.

[0095] In this embodiment, the load condition sensitivity evaluation bins are only used to determine whether the differential characteristics of the corrected dual-stack model undergo systematic drift with changes in load condition variables. Considering that this step is a preliminary load condition sensitivity evaluation step, and the number of bins will be determined separately in the subsequent dynamic threshold boundary database construction stage, this embodiment divides the load condition variables corresponding to the normal training samples into 4 load condition sensitivity evaluation bins according to the quantile binning method, that is, using the 0%, 25%, 50%, 75%, and 100% quantiles as the load condition interval boundaries.

[0096] If adjacent quantiles correspond to the same operating condition variable values, or if the sample size in a bin is insufficient to stably calculate the bin mean and bin standard deviation, then the bin should be merged with the adjacent bins and the bin statistics should be recalculated.

[0097] The analytical indicators include Pearson correlation coefficient, Spearman rank correlation coefficient, mean drift index, fluctuation scale change index, and between-group explanatory power index. These indicators are used to determine whether the corrected variance characteristics have undergone systematic drift with changes in operating conditions.

[0098] The results of the analysis of the relationship between the corrected difference characteristics and the operating condition variables are shown in Table 4.

[0099] Table 4

[0100] Sensitivity analysis ranking results of different differential characteristics are as follows: Figure 2 As shown.

[0101] Based on the above results, and There is a clear coupling relationship between it and the operating condition variables, so a dynamic threshold boundary is adopted; , and The system is less sensitive to operating conditions and uses a global static boundary.

[0102] S42. Characteristic flow division and boundary type determination; Based on the operating condition sensitivity evaluation vector obtained in S41, the operating condition sensitivity of each corrected dual-reactor difference characteristic is determined. The operating condition sensitivity evaluation vector includes the absolute value of Spearman correlation coefficient, the absolute value of Pearson correlation coefficient, the inter-group explanatory power, the mean drift index, and the fluctuation scale change index.

[0103] To avoid the subjectivity of setting thresholds based on human experience, the values ​​used for the initial determination of operating condition sensitivity are determined adaptively using quantiles, based on the statistical distribution characteristics of all differential features in the normal training samples. Taking this embodiment as an example: 1. Calculate the sensitivity indices for all features. For each corrected bi-group differential feature, calculate its Spearman correlation coefficient absolute value, Pearson correlation coefficient absolute value, between-group explained value η², mean shift ratio, and fluctuation scale change ratio according to method S41.

[0104] 2. Determine the candidate thresholds for each indicator. Sort all the above features by their values ​​on a specific indicator (such as the absolute value of the Pearson correlation coefficient) from smallest to largest. Select the upper quartile of the indicator sequence as the strong sensitive candidate threshold, and the median as the weak sensitive candidate threshold.

[0105] 3. A comprehensive set of judgment criteria is established. The strong sensitivity threshold for the absolute value of the Spearman correlation coefficient is 0.35, and the weak sensitivity threshold is 0.20; the strong sensitivity threshold for the absolute value of the Pearson correlation coefficient is 0.30, and the weak sensitivity threshold is 0.15; the strong sensitivity threshold for the inter-group explanatory power is 0.14, and the weak sensitivity threshold is 0.06; the strong and weak sensitivity thresholds for the mean drift index and the fluctuation scale change index are also determined based on the upper quartile and median of the corresponding index series. This forms the judgment conditions in S42.

[0106] In this embodiment, the minimum number of normal training samples allowed for each work condition bin is preset to 30, ensuring that there are enough samples in each work condition bin to calculate statistics such as quantile boundaries, median, and interquartile range. The minimum effective width coefficient is preset to 0.01. The minimum effective width is obtained by multiplying the overall feature fluctuation scale in the normal training samples by this coefficient. It is only used to remove obvious invalid boundaries, but is not used as the final diagnostic threshold itself.

[0107] The traffic splitting results are shown in Table 5 below: Table 5

[0108] I. Construction of dynamic threshold boundaries for operating condition sensitive features; Sensitive characteristics of operating conditions identified in S42 and In this embodiment, normal training samples are used to construct a dynamic threshold boundary that varies with operating condition variables, and normal validation samples are used to automatically select the candidate boundary parameters.

[0109] (1) Determine the number of boxes; First, using current as the operating condition variable, the normal training samples are divided into several operating condition bins. A candidate bin size set is generated based on the number of normal training samples, satisfying the minimum bin size constraint, as follows: .

[0110] (2) Construction of feature sets within bins; Secondly, within each operating condition bin, the distribution of normal training samples for operating condition sensitive features is statistically analyzed, and the center position, fluctuation scale, inner quantile boundary, and outer quantile boundary are calculated.

[0111] The center position can be the median (or the average in a specific scheme), and the fluctuation scale can be the interquartile range (or the absolute deviation of the median in a specific scheme). The inner boundary is used to characterize the main fluctuation range of normal samples, and the outer boundary is used for subsequent binary classification of normal and faulty samples.

[0112] (3) Construction of inner and outer boundaries; This embodiment constructs candidate dynamic threshold boundaries based on normal training samples. The candidate threshold boundary parameter combination consists of the number of candidate bins, candidate inner layer quantile parameters, and candidate outer layer quantile parameters.

[0113] Then, regarding the number of candidate bins and candidate quantile parameters, a set of candidate dynamic threshold boundaries is generated on normal training samples for each set of candidate parameters, and the inner and outer layer coverage, coverage deviation, and boundary validity are calculated on normal validation samples. Boundary validity is a necessary constraint; when a candidate boundary exhibits overlapping inner and outer layer boundaries, boundary reversal, outer layer boundary not including inner layer boundary, boundary width lower than the minimum effective width, or insufficient bin samples, the candidate parameter combination is discarded.

[0114] ; ; To avoid abrupt threshold changes at the boundaries of adjacent compartments, this embodiment further performs linear interpolation on the baseline center, inner boundary, and outer boundary corresponding to the center positions of compartments under different operating conditions, generating a dynamic threshold curve that continuously changes with the operating condition variables.

[0115] Based on the boundary validity check results, the normal verification sample coverage evaluation results, and the diagnostic performance evaluation results, the target threshold boundary parameter combination is determined from the candidate threshold boundary parameter combinations, and the dynamic threshold boundary corresponding to the working condition sensitive feature is generated based on the target threshold boundary parameter combination.

[0116] The optimal parameters for the dynamic threshold library are shown in Table 6.

[0117] Table 6

[0118] Finally, among the candidate threshold boundary parameter combinations, candidate combinations with boundary inversion, abnormal inner and outer boundary inclusion relationships, or boundary widths lower than the minimum effective width are first eliminated based on the boundary validity check results. Then, combining the normal validation sample coverage evaluation results and the validation set diagnostic performance evaluation results, the target threshold boundary parameter combination is determined from the remaining candidate combinations. The target threshold boundary parameter combination is used to determine the number of bins, inner quantile parameters, and outer quantile parameters corresponding to the operating condition sensitive features, and accordingly generates a dynamic threshold boundary that changes with the operating condition variables.

[0119] In this manner, the number of candidate bins, candidate inner quantile parameters, and candidate outer quantile parameters serve only as candidate search ranges during the dynamic threshold boundary construction process, rather than being manually specified final threshold boundary parameters. The final number of bins and quantile parameters are determined jointly by boundary validity checks, normal validation sample coverage evaluation, and diagnostic performance evaluation. The determined dynamic threshold boundary, along with its corresponding operating condition bin information, quantile parameters, and boundary type, is written into the adaptive threshold boundary library described in S43 for subsequent boundary calls and fault determination of samples to be diagnosed.

[0120] To further illustrate how the dynamic boundary changes with operating conditions, For example, the dynamic boundary results of partial working conditions for the sub-bins are shown in Table 7 below.

[0121] Table 7

[0122] Operating condition sensitive characteristics A schematic diagram of the construction of dynamic threshold boundaries is shown below. Figure 3 As shown.

[0123] II. Construction of global static boundaries for non-operating condition sensitive features; for , and Because of its weak coupling with operating condition variables, this embodiment constructs a static boundary globally based on normal training samples. Specifically, the center position, fluctuation scale, inner quantile boundary, and outer quantile boundary of each non-operating condition sensitive feature in the normal training samples are statistically analyzed to form candidate global static boundaries. Then, normal validation samples are used to check the validity of the candidate global static boundaries and evaluate their coverage effect. The inner layer coverage rate, outer layer coverage rate, and corresponding coverage rate deviation are calculated, and the target quantile parameter combination corresponding to the global static boundary is determined by combining the validation set diagnostic performance evaluation results. The results of the global static boundary database construction parameters are shown in Table 8 below.

[0124] Table 8

[0125] S43. Generate an adaptive threshold boundary library; In this embodiment, the dynamic threshold boundaries and global static boundaries obtained in S42 are uniformly organized to form an adaptive threshold boundary library. A summary table of the unified boundary library is shown in Table 9. This boundary library is invoked in S5 for feature-level state determination under different characteristics and operating conditions.

[0126] Table 9

[0127] S5. Fault determination and parameter optimization based on adaptive threshold boundaries; S51, Adaptive Boundary Calling and Diagnostic Feature Generation; For the window to be diagnosed, the diagnostic values ​​of each corrected difference feature within the window are first calculated. For operating condition-sensitive features, the dynamic threshold boundary library is called based on the operating condition variable values ​​of the current window; for non-operating condition-sensitive features, the global static boundary library is called. For operating condition-sensitive features, the dynamic threshold boundary is called; for non-operating condition-sensitive features, the global statistical results are used as the reference boundary.

[0128] S52, Feature-level binary classification determination; This embodiment uses the outer boundary as the basis for binary classification between normal and faulty conditions. Let the diagnostic value of the i-th differential feature in the k-th window to be diagnosed be... Its corresponding outer boundary is .when When the time is right, the characteristic is considered normal; when or When this occurs, the feature is determined to be abnormal.

[0129] The inner boundary is only used to characterize the degree of slight deviation and continuous deviation, and is not used as the final fault trigger condition.

[0130] S53, Optimization of window-level aggregation rules; Based on the feature-level determination results, the abnormal states of multiple differential features within the same window are aggregated. In this embodiment, the final window-level aggregation rule can be expressed as shown in Table 10: Table 10

[0131] S54 validation set evaluation; To further illustrate the diagnostic performance evaluation method of the candidate threshold boundary parameter combination in S42 and to determine the window-level aggregation rule parameters in S53, this embodiment constructs candidate parameter combinations and calculates the diagnostic performance of each candidate parameter combination based on the validation set samples. The candidate parameter combinations include the candidate threshold boundary parameters used in S42 to generate adaptive threshold boundaries, and the candidate aggregation rule parameters used in S53 for window-level aggregation determination; wherein the candidate threshold boundary parameters include at least one of dynamic threshold boundary parameters and global static boundary parameters, and the candidate aggregation rule parameters include at least one of anomaly feature quantity threshold and anomaly feature proportion threshold.

[0132] Specifically, each candidate parameter combination is substituted into the diagnostic process from S51 to S53 to obtain the window-level prediction result corresponding to the validation set samples. Using the known state labels of the validation set samples as a reference, the diagnostic performance index corresponding to that candidate parameter combination is calculated. The diagnostic performance index includes at least one of accuracy, F1 score, false positive rate, and false negative rate. Accuracy represents the overall correctness of the judgment, F1 score comprehensively reflects the fault identification capability, false positive rate represents the proportion of normal samples misclassified as faulty, and false negative rate represents the proportion of faulty samples misclassified as normal.

[0133] Based on the diagnostic performance metrics, the target parameter combination is determined from the candidate parameter combinations. As a preferred approach, the accuracy and F1 score of the candidate parameter combinations on the validation set are maximized, while the false positive rate and false negative rate are minimized. Pareto non-dominated ranking is used to screen Pareto front solutions that are not dominated by other candidate solutions. When multiple Pareto front solutions exist, the target parameter combination is determined according to the priority order of lower false positive rate, lower false negative rate, and higher F1 score. Alternatively, a weighted comprehensive score or validation set performance ranking method can also be used to determine the target parameter combination.

[0134] The target parameter combination is used to determine the target threshold boundary parameter combination corresponding to the adaptive threshold boundary in S42, and the abnormal feature quantity threshold and abnormal feature proportion threshold corresponding to the window-level aggregation rule in S53. The above target parameter combination is solidified into the final diagnostic configuration for automatic diagnosis of subsequent samples to be diagnosed.

[0135] Based on the above final diagnostic configuration, window-level binary classification diagnosis was performed on the validation set samples to obtain the window-level diagnostic results, as shown in Table 11.

[0136] Table 11

[0137] Window-level diagnostic results analysis, such as Figure 4 As shown.

[0138] For any parts not mentioned above, existing technologies can be adopted or referenced.

[0139] Of course, the above description is only a preferred embodiment of the present invention. The present invention is not limited to the above-described embodiments. It should be noted that any equivalent substitutions or obvious modifications made by those skilled in the art under the guidance of this specification fall within the scope of this specification and should be protected by the present invention.

Claims

1. A method for fault diagnosis of fuel cells based on dual-stack consistency constraints, characterized in that... Includes the following steps: S1. Collect raw data during the operation of the dual-stack fuel cell, preprocess the raw data, and generate samples; S2. Based on the samples generated in S1, construct a dual-pile information mapping table, calculate the dual-pile difference features, evaluate the state representation ability of the dual-pile difference features, and determine the retained dual-pile difference features. S3. Based on the retained dual-stack difference features determined in S2, a global bias baseline is constructed using samples under normal operating conditions, and the retained dual-stack difference features are biased and corrected based on the constructed global bias baseline to obtain the corrected dual-stack difference features. S4. Conduct a working condition sensitivity analysis on the corrected dual-stack difference characteristics obtained in S3, divide them into working condition sensitive features and non-working condition sensitive features, and construct dynamic threshold boundaries and global static boundaries respectively to form an adaptive threshold boundary library. S5. Acquire and process the dual-stack operation data of the fuel cell to be diagnosed, and determine the fault based on the adaptive threshold boundary library generated in S4.

2. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 1, characterized in that, S1 includes the following steps: S11, Data Acquisition; The raw time-series data of the dual-stack fuel cell system during operation are collected. The raw time-series data includes the operating parameters of the corresponding measurement points of the dual stacks and the operating condition variables that reflect the changes in the system operating conditions. S12. Preprocess the raw time series data; The raw time-series data is cleaned, aligned, segmented, and standardized to generate standardized input data. S13, Generate samples; For standardized input data, multiple samples are generated using the overlapping sliding window method. Each sample includes the sequence of operating parameters of the corresponding measuring points in the two piles within a sliding window and the corresponding sequence of operating condition variables. S14. Sample partitioning; The samples are divided into training and validation sets.

3. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 2, characterized in that, S2 includes the following steps: S21. Construct a dual-heap information mapping table; Based on the samples, a mapping relationship between corresponding measurement points in a dual-stack fuel cell system is established, a dual-stack information mapping table is constructed, and then... Indicates the difference characteristics between the two piles. Indicates the first Each pair of double-stack corresponding measurement points; S22. Calculate the difference characteristics between the two piles; Based on the constructed dual-stack information mapping table, the difference between corresponding measurement points is calculated, and a difference feature vector is constructed; for the th There are corresponding measurement points, at time... The dual-stack difference characteristic is defined as: ； in, Indicates the first heap in heap 1 The measured values ​​at each measuring point at time t. This indicates the time at which the corresponding measurement point in pile 2 is located. The measured value, Indicates the first The dual-pillar difference characteristics at each measuring point at time t; S23. Evaluate the state characterization capability of the dual-pillar difference characteristics and determine the retained dual-pillar difference characteristics; After obtaining multiple dual-stack difference characteristics, a preliminary judgment is made on the correlation between each dual-stack difference characteristic and the changes in fuel cell operating status. Dual-stack difference characteristics that do not show any difference changes in normal and fault states are eliminated, and the retained dual-stack difference characteristics are obtained.

4. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 3, characterized in that, S3 includes the following steps: S31. Determine the global bias baseline; No. The first sample The bias is defined as follows: ； in, Indicates the first The number of valid sample points in a sample. Indicates the first In the nth sample The two-stack difference features at time... The value of , Indicates the first The first sample The offset of each pair of measurement points corresponding to the double stack; The biases of all normal training samples are summed in a second order. The global bias baseline for the i-th dual-stack differential feature is: ； in, Represents the total number of samples. Indicates the first Global bias baseline for dual-stack difference characteristics; S32. Calculate and correct the difference characteristics between the two stacks; Based on the global bias baseline determined in S31, the retained dual-stack difference features obtained in step S2 are subjected to bias correction, and the calculation formula is as follows: ； in, For the first In the nth sample The two-stack difference features at time... The value of , For the first Global bias baseline for dual-stack difference characteristics To correct for the difference characteristics between the two stacks.

5. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 4, characterized in that, S4 includes the following steps: S41. Conduct a working condition sensitivity analysis; The correlation index between the corrected dual-reactor difference characteristics and the operating condition variables is calculated and binned to obtain the operating condition sensitivity evaluation vector. Based on the operating condition sensitivity evaluation vector, the corrected dual-reactor difference characteristics are divided into operating condition sensitive characteristics and non-operating condition sensitive characteristics. S42. Characteristic flow division and boundary type determination; For operating condition-sensitive features, a dynamic threshold boundary is constructed that varies with the operating condition; for non-operating condition-sensitive features, a global static boundary is constructed. S43. Generate an adaptive threshold boundary library; The adaptive threshold boundary library is formed by storing the dynamic threshold boundary corresponding to the working condition sensitive feature and its corresponding working condition binning information and quantile parameters, the global static boundary corresponding to the non-working condition sensitive feature and its corresponding quantile parameters, as well as the feature number and boundary type corresponding to the dynamic threshold boundary and the global static boundary.

6. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 5, characterized in that, S41 includes the following steps: Calculate the Pearson linear correlation coefficient and the Spearman rank correlation coefficient to characterize the degree of linear correlation and monotonicity. ； In the formula, This represents the mean of the corrected dual-stack difference characteristics. Indicates time Operating condition variables, This represents the mean of the operating condition variables. The Pearson linear correlation coefficient; This indicates the rank of the corrected bi-hash differential feature across all samples. Indicates time The rank of the working condition variable in all samples. Indicates the first The mean of the ranks of the differential features of each corrected double-stack. The mean of the rank of the operating condition variable. This represents the Spearman rank correlation coefficient; Binning was performed, dividing the operating condition variables corresponding to the normal training samples into several operating condition sensitivity evaluation bins according to the quantile binning method. ; For each load condition sensitivity evaluation bin, the number of samples, bin mean, and bin standard deviation of the i-th corrected dual-stack difference feature within that bin are calculated. The changes in bin mean, bin standard deviation, and inter-group explanatory power are compared between different load condition sensitivity evaluation bins to determine whether the changes in the corrected dual-stack difference feature across different load condition intervals exhibit systematic drift with changes in load conditions. The changes in the corrected dual-stack difference feature across different load condition intervals are described by the number of samples, bin mean, and bin standard deviation within each load condition sensitivity evaluation bin, as well as the population center position and population standard deviation on all normal training samples. The specific calculation process is as follows: (1) Mean drift index: ； (2) Changes in fluctuation scale: ； (3) Intergroup explanatory power index: ； Where b represents the working condition compartment number. This indicates the mean shift indicator; This represents the bin mean of the i-th corrected dual-stack difference feature within the b-th bin; This represents the maximum value among all the sub-bins under all operating conditions; This represents the minimum value of the average across all operating conditions in the binning area; Indicates the first The overall standard deviation of each corrected dual-stack difference feature across all normal training samples; It represents a very small positive number and is used to prevent the denominator from being zero; Indicators representing changes in fluctuation scale; Indicates the first The corrected dual-stack difference feature in the first Standard deviation within each operating condition compartment; This represents the maximum standard deviation within all operating conditions across all bins. This represents the minimum standard deviation within all operating conditions. Indicates the total number of compartments under different operating conditions; This represents the total number of normal training samples. Indicates the first The number of normal training samples in each working condition bin; Indicates the first In the nth normal training sample The values ​​of the corrected dual-stack difference characteristics; This represents the overall center position of the i-th corrected bi-hash differential feature across all normal training samples; Indicates the inter-group explanatory power index; Using Pearson linear correlation coefficient, Spearman rank correlation coefficient, mean drift index, fluctuation scale change index, and between-group explanatory power index, the working condition sensitivity evaluation vector for the i-th corrected dual-stack difference characteristics is formed: ； Based on the comparison results of each index in the operating condition sensitivity evaluation vector with the corresponding sensitivity judgment threshold, the corrected dual-stack difference characteristics are divided into operating condition sensitive characteristics and non-operating condition sensitive characteristics. For each operating condition sensitivity index, the value of the corrected dual-stack difference characteristic on that index is statistically analyzed and sorted according to the numerical value. The upper quartile of the sorting result is taken as the strong sensitivity judgment threshold, and the median of the sorting result is taken as the weak sensitivity judgment threshold. When any operating condition sensitivity index reaches the corresponding strong sensitivity judgment threshold, or at least two operating condition sensitivity indices reach the corresponding weak sensitivity judgment threshold, the corresponding corrected dual-stack difference characteristic is classified as an operating condition sensitive characteristic; otherwise, it is classified as a non-operating condition sensitive characteristic.

7. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 6, characterized in that, Constructing a dynamic threshold boundary that varies with operating conditions in S42 includes the following steps: (1) Determine the number of bins; divide the normal training samples into several working condition bins according to the working condition variables, generate a candidate bin number set according to the number of normal training samples, and satisfy the minimum bin sample number constraint, as follows: ； in, ; This refers to the number of boxes. This is the normal training sample window size. The minimum number of samples allowed per bin. To minimize the number of boxes, This represents the maximum number of boxes. (2) Construct the feature set within the bins; For the i-th sensitive feature of the working condition and the b-th bin of the working condition, the corrected double-stack difference feature values ​​of the working condition variable values ​​in the normal training samples that fall into the bin of the working condition are extracted to form a bin feature set; the median and interquartile range are calculated based on the bin feature set and are used as the center position and fluctuation scale under the bin of the working condition, respectively. (3) Construct the inner quantile boundary and the outer quantile boundary; Construct inner and outer quantile boundaries; let the inner and outer quantile parameters be: and And satisfy Then, the inner layer partition boundary of the sensitive feature of the i-th working condition under the b-th working condition is: ； in, Indicates the first The sensitive feature of the working condition is in the first The lower inner boundary of each working condition compartment Indicates the first The sensitive feature of the working condition is in the first The inner upper boundary of the compartment under each working condition. Indicates the first The sensitive feature of the working condition is in the first The set of feature values ​​of normal training samples within each working condition bin; The outer boundary is: ； in, Indicates the first The sensitive feature of the working condition is in the first The outer lower boundary of the compartment under each working condition. Indicates the first The sensitive feature of the working condition is in the first The outer upper boundary of the sub-compartment under each working condition; Each working condition bin corresponds to a set of local normal boundaries. When a sample enters the diagnostic stage, the boundary under the corresponding bin is called to make a judgment based on the working condition bin to which its current working condition variable value belongs. (4) Dynamic threshold boundary interpolation; When the operating condition variable in the window to be diagnosed is located in an adjacent operating condition bin and When there is a gap, interpolate the boundary parameters corresponding to the adjacent working conditions of the sub-bins; ； in, Indicates the first Individual working condition distribution center Indicates the first Individual working condition distribution center Indicates the first The operating condition variable values ​​for each window to be diagnosed. Indicates the first The boundary parameters of each working condition sensitive feature at the center of the adjacent lower compartment. Then it means the first The boundary parameters of a working condition sensitive feature at the center of an adjacent upper compartment; the boundary parameters include at least one of the inner lower boundary, inner upper boundary, outer lower boundary, and outer upper boundary.

8. The fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 7, characterized in that, The construction of dynamic threshold boundaries that vary with operating conditions in S42 also includes the following steps: (5) Boundary validity check; The validity of candidate dynamic threshold boundaries is checked, and candidate threshold boundary parameter combinations that have boundary inversion, abnormal inner and outer boundary inclusion relationship, or boundary width is lower than the minimum effective width are eliminated; wherein, the boundary width is the difference between the upper boundary and the lower boundary, and the minimum effective width is determined by multiplying the feature overall fluctuation scale in normal training samples by a preset width coefficient. (6) Evaluation of normal validation sample coverage; For candidate dynamic threshold boundaries that pass the validity check, coverage is evaluated based on normal validation samples. The actual proportions of normal validation samples falling into the inner and outer boundaries are as follows: ； And calculate the coverage deviation based on the difference between the actual coverage and the theoretical target coverage: ； in, This indicates the actual coverage rate of normal validation samples falling into the inner layer boundary. This indicates the actual coverage rate of normal validation samples falling into the outer boundary. To ensure a normal number of validation samples, For the first In the first normal verification sample One corrected dual-stack difference eigenvalue, Let K be the operating condition variable corresponding to the k-th normal verification sample. and These represent the i-th feature in the operating condition variables. The corresponding dynamic inner lower boundary and dynamic inner upper boundary are below. and These represent the i-th feature in the operating condition variables. The corresponding dynamic outer lower boundary and dynamic outer upper boundary are below. This indicates an indicator function that takes the value 1 if the condition within the parentheses is true, and 0 otherwise. The coverage rate of the inner layer theoretical target. , For the coverage of the outer theoretical target, , and These are the upper and lower quantile parameters of the inner layer, respectively. and These are the upper and lower boundary parameters of the outer layer, respectively. This is due to coverage deviation; (7) Determination of target parameters for the dynamic threshold library; Based on the boundary validity check results, coverage deviation, and diagnostic performance evaluation results, the target threshold boundary parameter combination is determined from the candidate threshold boundary parameter combinations, and the dynamic threshold boundary corresponding to the working condition sensitive feature is generated based on the target threshold boundary parameter combination. The candidate threshold boundary parameter combination consists of the number of candidate bins, the candidate inner layer quantile parameter, and the candidate outer layer quantile parameter; the diagnostic performance evaluation result is the diagnostic performance index calculated based on the known state labels of the validation set samples after the candidate threshold boundary parameter combination is applied to the validation set samples through boundary invocation and feature-level judgment.

9. A fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 8, characterized in that, Constructing a global static boundary in S42 Includes the following steps: Let the first The set of non-operating condition sensitive features in the normal training samples is as follows: ； in, This represents the global feature set of the i-th non-condition-sensitive feature in the normal training samples. This represents the value of the i-th non-condition-sensitive feature in the k-th normal training sample. This represents the normal training sample set; Compute on a set: Central position : This is used to characterize the global baseline level of this non-operating condition sensitive feature in normal training samples; Represents the median function; Fluctuation Scale : This is used to characterize the global dispersion of the non-operating condition sensitive feature in normal training samples; Represents the standard deviation function; Constructing inner and outer boundaries based on quantiles: Let the inner quantile parameters and outer quantile parameters be respectively: and And satisfy Then the first The global static inner layer quantile boundary of each feature is: ； The outer boundary is: ； in, and Let represent the lower boundary and upper boundary of the global static inner layer of the i-th non-operation condition sensitive feature, respectively; and Let $i$ represent the lower boundary and upper boundary of the global static outer layer of the i-th non-operational condition sensitive feature, respectively. For candidate global static boundaries, normal validation samples are used to check the boundary validity and evaluate the coverage effect. Based on the check and evaluation results, the target quantile parameter combination corresponding to the global static boundary is determined. Finally, the target quantile parameter combination corresponding to the global static boundary is written into the adaptive threshold boundary library.

10. A fuel cell fault diagnosis method based on dual-stack consistency constraints according to claim 9, characterized in that, S5 includes the following steps: S51, Adaptive Boundary Calling and Diagnostic Feature Generation; For the dual-stack operation data of the fuel cell to be diagnosed, the preprocessing and sliding window method in S1 is used to generate the sample to be diagnosed. The dual-stack difference characteristics of the sample to be diagnosed are calculated according to S2. The bias correction is performed based on the global bias baseline according to S3 to obtain the corrected dual-stack difference characteristics of the sample to be diagnosed. Based on the operating condition sensitivity analysis results and adaptive threshold boundary library obtained in S4, dynamic threshold boundaries are called for the operating condition sensitive features according to their corresponding operating condition variables, and global static boundaries are called for the non-operating condition sensitive features to generate window-level diagnostic features for fault determination. S52, Feature-level determination; Based on the window-level diagnostic features generated by S51, a primary discrete state determination is performed on the operating state of the fuel cell system, as follows: Let the diagnostic value of the i-th corrected double-stack difference feature in the k-th window to be diagnosed be... Its corresponding inner threshold boundary is The outer threshold boundary is The inner threshold boundary is used to characterize the normal state fluctuation area, and the outer threshold boundary is used as the basis for judging whether there is a fault. The inner threshold boundary and the outer threshold boundary are dynamic threshold boundaries or global static boundaries obtained according to the S51 call. when If so, the feature is determined to be in a normal state; when or If so, the feature is determined to be in a faulty state.

Images