Method, system and device for predicting remaining service life of proton exchange membrane fuel cell

Abstract

The invention provides a method, a system and a device for predicting the residual service life of a proton exchange membrane fuel cell. The method comprises the following steps: step 1, collecting monitoring data of a charge-discharge process of the fuel cell; 2, carrying out feature selection on the monitoring data of the charging and discharging process; step 3, using robust local mean decomposition based on dynamic time warping improvement to decompose the monitoring data of the charging and discharging process of the fuel cell; 4, reconstructing a sub-sequence by adopting an entropy control principal component polymerization method; 5, establishing a dynamic adaptive weighted extreme learning machine of the life prediction model and the residual service life prediction model; step 6, improving a projection iterative optimization algorithm; and 7, optimizing the dynamic adaptive weighted extreme learning machine of the residual service life prediction model to obtain a final prediction result. The problem that the residual service life of the fuel cell cannot be accurately predicted is solved, and the prediction accuracy of the residual service life of the cell is improved.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and specifically to a method, system, and apparatus for predicting the remaining service life of a proton exchange membrane fuel cell. Background Technology

[0002] Proton exchange membrane fuel cells (PEMFCs) have become a core candidate technology for new energy vehicles, distributed power stations, and portable power supplies due to their advantages such as zero emissions, high efficiency, and rapid low-temperature start-up. However, the electrochemical reactions inside the stack are complex, and key components such as membrane electrodes, catalysts, and bipolar plates undergo irreversible degradation during long-term operation, leading to a gradual decline in output performance and ultimately system failure. Currently, the industry generally adopts a "scheduled maintenance" or "replacement after failure" strategy, which not only has high operation and maintenance costs but also may cause safety accidents due to sudden failures. Therefore, achieving accurate prediction of the remaining useful life (RUL) of PEMFCs and subsequently developing proactive operation and maintenance and health management strategies has become a bottleneck that must be overcome for the commercialization of fuel cells.

[0003] Existing RUL prediction methods can be broadly categorized into three types: model-based, data-driven, and hybrid. Model-based methods rely on simplified electrochemical-thermal-mechanical coupling equations, making it difficult to characterize the multi-scale coupling mechanism between material aging and performance degradation. Traditional data-driven methods, while capable of building statistical or shallow machine learning models using monitoring data, are limited by weak feature extraction capabilities, sensitivity to non-stationary operating conditions, and insufficient capture of long-sequence dependencies. In recent years, deep learning methods have been introduced, but most directly apply general-purpose networks such as CNNs and LSTMs without fully considering the nonlinear, non-stationary, and multimodal characteristics of fuel cell signals. Furthermore, hyperparameters rely on empirical tuning, resulting in insufficient model generalization and interpretability. Therefore, a dedicated prediction framework for PEMFCs is urgently needed, capable of adaptively mining degradation features under complex operating conditions, reducing computational complexity, and improving prediction accuracy. Summary of the Invention

[0004] Purpose of the invention: The technical problem to be solved by the present invention is to provide a method, system and device for predicting the remaining service life of a proton exchange membrane fuel cell, which can effectively improve the accuracy of the prediction of the remaining service life of the battery, in order to overcome the shortcomings of the prior art.

[0005] The method includes the following steps: Step 1: Deploy micro-sensors within the fuel cell stack to continuously collect monitoring data on the charging and discharging process of the fuel cell; Step 2: The consensus feature selection method (CFS) is used to select features from the charging and discharging process monitoring data. Feature variables with high relevance are selected as input variables. The consensus feature selection method (CFS) includes: first, constructing shadow features from the charging and discharging process monitoring data using the Boruta feature selection algorithm for preliminary importance assessment; then, performing secondary screening using the Monte Carlo Tree Search-eXtreme Gradient Boosting (MCTS-XGBoost) algorithm, which improves upon Monte Carlo Tree Search, and retaining only features that are recognized by both secondary screenings. Step 3: Decompose the fuel cell charging and discharging process monitoring data into two or more subsequences using the robust local mean decomposition based on dynamic time warping (DTW). Step 4: Reconstruct the subsequence using the entropy-controlled principal component aggregation method to reduce computational complexity; Step 5: Establish the lifespan prediction model TST-former, use the lifespan prediction model TST-former to extract battery aging features in stages, and output multi-scale health factors as the first layer prediction results, i.e. equipment degradation feature data; establish the remaining lifespan prediction model Dynamic Adaptive Weighted Extreme Learning Machine DAVW-ELM, combine all the prediction values ​​of TST-former to reconstruct a dynamic evolution feature set, and use it as the input of the remaining lifespan prediction model Dynamic Adaptive Weighted Extreme Learning Machine DAVW-ELM; Step 6: Improve the projection-iterative-methods-based Optimizer (PIMO) algorithm: use Logistic chaotic mapping to initialize the population of the projection-iterative optimization algorithm; introduce a spiral flight search strategy in the algorithm position update stage, and ensure that the solution is always within the feasible region through the projection operator, thus obtaining the improved projection-iterative optimization algorithm TPIMO; Step 7: The weights and biases of the dynamic adaptive weighted extreme learning machine DAVW-ELM for predicting the remaining useful life are optimized using the improved projection iterative optimization algorithm TPIMO to obtain the optimal weights and biases. The optimized dynamic adaptive weighted extreme learning machine DAVW-ELM is then used to fuse and predict the test data samples of the proton exchange membrane fuel cell to obtain the final prediction result.

[0006] Step 1 includes: A 3×3 grid-like micro-sensor array is deployed in different areas of the fuel cell stack (inlet / outlet, single-cell series node) to simultaneously collect 12 data points: voltage, temperature, humidity, and gas pressure. A humidity sensor, gas flow meter, and membrane impedance monitoring module are also deployed to continuously collect data on the fuel cell's charging and discharging process. A dynamic adaptive sampling strategy with triggered sampling is employed, and a threshold is set (e.g., voltage fluctuation > 5%, temperature gradient > 5%). Set to high-frequency sampling (1kHz), and normally maintain low-power mode (10Hz). Ultimately, new data on fuel cell charging and discharging voltage, current, current density, temperature, capacity, internal resistance, cycle count, aging time, and output power are obtained.

[0007] Step 2, the Boruta feature selection algorithm includes the following steps: Step 2-a1: Boruta first copies the original data feature matrix sample set R with m rows and n columns to obtain a copy sample. Then, it performs a random row transformation on each column of the copy sample to obtain a shadow feature sample L with m rows and n columns. Finally, it concatenates R and L horizontally to obtain a new matrix N=[R,L] with dimensions m×2n, which serves as the input for the feature importance evaluation of the subsequent random forest. Finally, it outputs the optimal feature subset S. Step 2-a2: Run a random forest classifier on matrix N, mask individual features in sequence and record the average loss value of model accuracy Meanlmp; calculate the average importance MeanImp of all features using SHAP value theory, mark features with importance lower than Meanlmp as unimportant, and permanently delete them from matrix N. Step 2-a3: Extract the maximum value of the shadow feature importance, Maxlmp, and compare it with the importance score of the true feature. Among the true features, those with an importance greater than Maxlmp are marked as important, and other features are marked as tentative. Step 2-a4: Select the importance score of the optimal shadow feature as the initial screening threshold, gradually compare the importance of each true feature with the optimal shadow feature, retain the true features with higher scores than the optimal shadow feature, repeatedly remove features with lower scores than the optimal shadow feature, and continuously optimize the feature set. Steps 2-a5 continue until there are no more provisional features, or the preset number of iterations for the random forest is reached, at which point the algorithm stops.

[0008] Step 2, the MCTS-XGBoost algorithm based on Monte Carlo tree search improvement includes the following steps: Step 2-b1: Initialize a search tree, with the root node representing the optimal feature subset S; set the parameters of the Monte Carlo Tree Search (MCTS), including the number of searches N, the maximum depth of the tree D, the number of samples per node K, and the balance parameter C between search and utilization. Step 2-b2: XGBoost utilizes the feature selection capability of its tree structure to perform a secondary evaluation and ranking of the optimal feature subset S. The optimal feature subset S is used as input, and the feature weights are optimized through the built-in feature selection mechanism of XGBoost's tree structure. In each iteration, XGBoost adds a new regression tree to fit the residual of the previous prediction. Let the predicted value of the ensemble model with t-1 trees be... The newly built regression tree in round t The goal is to fit the current residual, and the iterative process is as follows: (1), in The whole represents the predicted value of the j-th sample in round 0 (initial state); Step 2-b3: Using XGBoost's feature importance evaluation mechanism, the importance of each feature in the optimal feature subset S is evaluated. When splitting a node in each tree, the optimal feature and the split node are selected based on the feature gain (Gain). The essence of the gain formula is to measure the degree to which the model's loss function decreases before and after the split. The gain formula is as follows: (2), in This indicates the gain or signal parameter associated with the left channel; Indicates the gain or signal parameter associated with the right channel; For regularization parameters; These are the weighting coefficients; α is a very small positive number used to ensure numerical stability; α and β are adjustable power parameters. and These are the feature parameters associated with the left channel and the feature parameters associated with the right channel, respectively. Step 2-b4: Based on the feature importance evaluation results, divide the optimal feature subset S into two or more subspaces, and use each subspace as a new node in the search tree; Step 2-b5: Use a Bayesian optimized surrogate model to partition the space, evaluate the node value based on the upper bound of the confidence interval UCB algorithm, calculate the upper bound of the confidence interval UCB value, and select the node with the highest upper bound of the confidence interval UCB value as the next exploration target. Step 2-b6: Within the feature subspace corresponding to the selected node, random sampling is performed again to generate a new feature set. By iteratively calculating the upper bound of the confidence interval (UCB) for each node, the search path is dynamically determined to obtain the optimal parameter configuration. The definition of the upper bound of the confidence interval (UCB) is: (3), Where K represents the degree of exploration during the search, and when C=0 it is a greedy strategy, it is easy to enter a local optimum; The number of sampling points for the parent node; This represents the number of sampling points for the current node. This represents the average value currently observed in the i-th feature subspace of node j; This represents the number of times that node j has been visited or sampled in the i-th feature subspace; Step 2-b7: Select the leaf node with the highest UCB value of the upper bound of the confidence interval from the search tree, and use it as the final feature set. Output the feature set corresponding to the leaf node with the highest UCB value of the upper bound of the confidence interval as the final result of XGBoost feature selection.

[0009] Step 3 includes the following steps: Step 3-1, Robust Local Mean Decomposition (RLMD) is a signal processing method based on Local Mean Decomposition (LMD) that optimizes boundary conditions, envelope estimation, and termination criteria. It preprocesses the original signal and extracts the solution signal. : (4), in It is a discrete signal; It is a local average function; It is a local envelope function; Robust Local Mean Decomposition (RLMD) is applied to the original signal. The target signal that is easier to decompose after preprocessing; Step 3-2, during the LMD iteration process, if the frequency of the frequency-modulated signal in two consecutive iterations satisfies w 1k+2 >w 1k+1 And w 1k+3 >w 1k+2 When the iteration stops, the result of the nth iteration is output, and the pure frequency modulated signal Y2(a) and the envelope signal h2(a) are obtained. The product function Z is obtained by multiplying the pure frequency modulated signal Y2(a) and the envelope signal h2(a). PF (a) The discrete signal X(a) is decomposed into f product function PF components and 1 residual signal. ;w 1k+2The frequency value of the frequency-modulated signal is calculated with an iteration step size of 1k+2. Step 3-3: Given two sets of PF sequences, which are sequences... and sequence Construct a dimension as Cost matrix Define the distance matrix ; in, express and of Norm (Euclidean distance); Represents a sequence The m-th element in Represents a sequence The nth element in; Represents the i-th element of sequence V eigenvectors; Represents the i-th element of sequence h eigenvectors; Steps 3-4: Construct a cumulative cost matrix , used in the cost matrix Calculate the minimum cumulative distance between two PF sequences h and V from the starting point (1,1) to the ending point (m,n), and define the weighted sum of the local cost measures as r(i, j); Steps 3-5, Cost Matrix The value of the bottom right element r(m,n) is the DTW distance between PF sequences. The smaller the distance, the greater the similarity between the sequences. Steps 3-6: Use DTW distance to screen out PF sequence clusters with high similarity; reconstruct a group of PF components from the sequences within the cluster, and then iteratively peel them off: if the DTW distance between a component and the other components is less than the set threshold range [40,60], then record the component as a new independent component PF1 and extract it. Repeat the iterative stripping operation on the remaining components until all PF components are separated into unrelated subsequences.

[0010] Step 4 includes the following steps: Step 4-1, the entropy-controlled principal component aggregation method introduces an augmented matrix improved kernel principal component analysis (KPCA) to process the high-dimensional correlation of the trajectory matrix, and then fuses singular spectral entropy to screen periodic or deterministic subsequences; Step 4-2, construct the trajectory matrix, assuming the original time series is... The length is N; where Let x represent the Nth element in the original time series x; select a window length L, and transform the one-dimensional sequence into a trajectory matrix X of L◊K using a sliding window; further aggregate the data through singular value decomposition and entropy analysis, remove noise and non-primary structures, and obtain representative subsequences: (5); Step 4-3: Assemble the sample points x into an input trajectory matrix X, where X contains N samples, and use a nonlinear mapping... Mapping sample point x to a high-dimensional space A new matrix is ​​obtained. ; Step 4-4, In higher dimensions We then perform dimensionality reduction using augmented matrix improved kernel principal component analysis (KPCA) to obtain the high-dimensional space. The covariance matrix is ​​calculated using a kernel function. Then calculate the matrix. The eigenvectors α corresponding to the larger eigenvalues ​​(calculated first, then compared) are obtained, and then the corresponding weight vectors are obtained. They are then arranged in descending order of the proportion of the eigenvalues, with the larger the proportion, the more important they are. The principal components are found in order of importance. Steps 4-5 involve constructing an augmented matrix that introduces nonlinear changes to expand the feature space of the original data, thereby enhancing the feature representation capability of KPCA. The augmented matrix... Represented as: (6), Where q represents the time delay length; Steps 4-6: Construct a segmented trajectory matrix, divide the dimensionality-reduced matrix into local sub-blocks, and process them separately; Steps 4-7: Perform singular value decomposition on the trajectory matrix X; Steps 4-8 involve normalizing the singular values ​​and calculating the proportion of each singular value to the sum of all singular values. This yields the proportion of each singular value to the total energy, which serves as the probability distribution. ; Steps 4-9: Calculate the singular spectral entropy. The larger the entropy value, the more uniform the distribution of singular values ​​and the higher the signal complexity (such as the presence of noise or nonlinear components). The smaller the entropy value, the more concentrated the signal energy is in a few singular values, indicating the existence of a significant periodic or deterministic structure. Based on the results of the entropy value analysis, reconstruct two or more representative subsequences. Steps 4-10: For each low-entropy submatrix Before selection Principal singular value reconstruction: (7), in, Represents the original low-entropy submatrix Approximate or reconstructed value; Let j be the singular value of the i-th submatrix; This represents the j-th left singular vector of the i-th submatrix; Let represent the j-th right singular vector of the i-th submatrix.

[0011] Step 5, establishing the lifetime prediction model TST-former, includes the following steps: Step 5-a1: Perform a Fast Fourier Transform on the aggregated and reconstructed data (load data). By utilizing the periodicity and symmetry of the signal to reduce the amount of computation, the time-domain signal is converted to the frequency domain. According to Passevar's theorem, calculate the energy spectral density in the frequency domain. After excluding the DC component, find the frequency index corresponding to the maximum energy. Convert the frequency index into a period to obtain the main period M of the load data.

[0012] Step 5-a2: Aggregate and reconstruct the sequence (load time series). The data is transformed into a two-dimensional trajectory matrix H; a window length W is selected, and the load data is intercepted and segmented into two or more segmented signals according to the window length W, and the segmented signals are embedded into the matrix to obtain the trajectory matrix H: (8), Wherein, the number of data segments K = l - W + 1 and K > W; This represents the k-th segment signal vector, which is a subsequence of length W extracted from the original sequence by sliding window. This represents the l-th data point in the original sequence; Step 5-a3: Perform singular value decomposition on the trajectory matrix H, decomposing it into the product of three matrices; calculate the contribution rate r of the singular spectral components, and extract noise components to reduce noise interference. The trajectory matrix X is further optimized after this processing. The contribution rate r and the filtered trajectory matrix... The formula is: (9), (10) in, It is a singular value. For the filtered first A grouped matrix; d is the upper limit of the effective rank; Step 5-a4, group each matrix Transform into a feature subsequence of length W ,in The w-th element in the representation; each feature subsequence represents a certain feature of the load data, such as long-term trend, seasonal trend, etc. Step 5-a5: Embed the time series, decomposing the time series into trend A and seasonal B components along the channel dimension; embed the seasonal components into E. B The trend is embedded into E by feeding the attention module and the feedforward network FFN through a subtraction mechanism. A Feed to interactive module A which monitors seasonal information. s and F s Finally, the seasonality and trend are combined to obtain the final prediction (focus module, feedforward network FFN, interaction module A). s and F s It is a component of the TST-former model. Step 5-a6, decompose the sequence data X A X B The sequence is mapped from the original space to a new space (the sequence reconstructed by entropy-controlled principal component aggregation is input into the TST-former model for FFT-SSD decomposition) to obtain X. A X B X A ∈R,X B ∈R, where R represents the set of real numbers; using direct linear layers and dropout layers to create a global covariate X. mark Trends and seasonal embeddings: (11), (12) Concat represents a connection operation. Indicates trend item embedding, Indicates seasonal item embedding; This represents a linear layer (fully connected layer) used to map the concatenated high-dimensional vector to the embedding space; Concat() uses a separate linear layer to contain... It is used on datasets containing information and different components.

[0013] Step 5-a7: Since the seasonal component without a moving average kernel can better highlight the intrinsic characteristics of time series data, the seasonal embedding E is added. s Feeding is provided to the attention and feedforward network FFN; Step 5-a8: Considering that untrainable moving average kernels can lead to unreliable trend patterns, seasonal information is fused through the Interactive Module (IM) to assist in learning trend branches (the Interactive Module is a component of the TST-former model); focus on weighted A S It effectively reflects the dependencies between multiple variables and can be transformed into a coefficient matrix to update the trend embedding; signals discarded by the seasonal branch are considered meaningful information to guide the representation of the trend embedding. (13) (14); in, This indicates that the trend embedding obtained in the previous step is the input representation of the trend branch; It is an activation function that compresses the output to the [0,1] interval, representing the weight or gating signal; This represents a 1×1 convolution, used for dimensionality reduction or channel transformation, mapping the output of multiple convolutions to a suitable dimension. Indicates multi-level convolution operations; This represents the fused trend representation, which combines seasonal information, trend embedding, and enhanced weights. Step 5-a9: At the end of the seasonal and trend output blocks, a gate mechanism σ is designed for the two streams to autonomously regulate the speed of information transmission; the gate mechanism for the seasonal and trend streams is represented as follows: (15) (16); in, Indicates the seasonal item output; Indicates the trend term output; , , , These represent different 1×1 convolutional layers used to extract or transform features; The intermediate feature representation of the seasonal branch (which may be the result of seasonal embedding after a certain layer of processing); Step 5-a10: Use the output of the previous TwinsBlock design overlay module as the input of the next TwinsBlock, overlay N TwinsBlocks to learn seasonal and trend representations, and then add the N TwinsBlocks together through linear projection to obtain the final prediction result, i.e., equipment degradation feature data.

[0014] Step 5, establishing the dynamic adaptive weighted extreme learning machine (DAVW-ELM) for predicting remaining useful life, includes the following steps: Step 5-b1: Combine all prediction results of TST-former to reconstruct a dynamic evolution feature set, which serves as the input to the dynamic adaptive weighted extreme learning machine DAVW-ELM for the remaining useful life prediction model. The remaining useful life prediction model, Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM), employs a dynamic adjustment strategy. It adjusts the corresponding input layer weights based on the output of each hidden layer neuron. This includes: calculating the activation level for each hidden neuron i and each training sample (combining all prediction results from TST-former to reconstruct a dynamically evolving feature set as the training sample); calculating the variance of the activation level for each hidden neuron; a larger variance indicates a stronger ability of the neuron to distinguish between different samples; and adjusting the input layer weights based on the activation level variance. Step 5-b2: Different weight update strategies have different impacts on model performance. The Dynamic Adaptive Weighted Extreme Learning Machine (DAWW-ELM) for predicting remaining lifetime uses an adaptive weight update strategy. It dynamically selects a suitable update strategy based on the training state of the DAWW-ELM and monitors the trend of training error changes. If the error decreases slowly, it indicates that the current strategy may be inefficient. Multiple weight update strategies can be selected, and a suitable weight update strategy can be dynamically chosen based on the trend of training error changes. Step 5-b3: The remaining useful life prediction model, Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM), introduces a variable weighting mechanism to differentiate between samples based on their importance. Important samples are given greater weights to ensure that the model can learn the features of these samples more accurately. Step 6 includes the following steps: Step 6-1: First, set the population size, dimension, number of iterations, and upper and lower limits of the search space for the Projection Iterative Optimization (PIMO) algorithm. Step 6-2: The Projective Iterative Optimization (PIMO) algorithm first initiates the optimization process by improving the initial population using a Logistic chaotic mapping. Given the population size S and the variable dimension G, each projective agent is considered a search agent for the algorithm, and the position of each agent is: (17) in This represents the k-th position of the i-th individual in the j-th dimension; and These represent the upper and lower bounds of the search space in the j-th dimension, respectively. This represents the initial value of the Logistic regression corresponding to the i-th individual, the j-th dimension, and the k-th position. Step 6-3: Calculate the fitness value of all search agents; Step 6-4: Determine if the fitness value meets the termination objective (if the absolute error between the current optimal search agent's fitness value and the theoretical optimal fitness value for this optimization problem is less than 10%). If the termination objective is met, the current optimal solution is output. If the objective is met, the current optimal search agent is returned as the final solution. If the objective is not met, the subsequent optimization steps are continued until the optimal fitness value is reached, and the optimal parameter configuration is output. Step 6-5, Residual Guided Projection (RGP): The algorithm repeatedly projects in the solution space, and randomly selects two individuals with the lowest residual fitness as guiding agents X. v1 and X v2 And according to X v1 and X v2 Calculate the gradient. The calculation formula is: (18) in, This is the best position at present; and It is a random number between [0,1]; C is a parameter used to adjust the weights of different parts when calculating the gradient, which can control the balance between exploration and utilization in the search process of the algorithm; Step 6-6, the Residual Guided Projection (RGP) stage introduces the Jacobian matrix J, where the element in the a-th row and b-th column of the Jacobian matrix J... Represent the components of the objective function Relative to variables The partial derivatives are then used; subsequently, the gradient and the Jacobian matrix J are incorporated into the projection step to complete the position correction. Steps 6-7, Double Random Projection Process (DRP): The algorithm first randomly selects two indices V3 and V4 from the particle swarm as reference points for projection updates. When randomly selecting gradient updates, the algorithm calculates the update direction vector k1 based on the difference between the current position and the current optimal position. (19) If we choose to use Jacobian matrix projection for updating, the update direction vector k2 is calculated as follows: (20) (twenty one), in, It is the new position after the DRP update. This indicates the step size for each iteration, controlling the "movement range" for each iteration; This represents the original position (current position) of the i-th individual before the update. Steps 6-8, Weighted Random Projection Update WRPU: Particle position updates are guided by random factors r7 and r8 and adaptive factor α, and random weighted projection and adaptive correction projection paths are dynamically selected; Steps 6-9, Lévy Flight Guided Projection (LFGP) process: This process achieves efficient global search and local optimization by simulating the random motion of organisms in nature. The trigger probability P of the Lévy flight projection is determined by the following formula: (twenty two); Where t represents the current iteration number and U represents the maximum iteration number; Steps 6-10 introduce a spiral flight search strategy into the Levy flight position update formula, combining the optimal solution. (The first three steps are Residual Guided Projection (RGP), Double Random Projection (DRP), and Weighted Random Projection Update (WRPU) to obtain the best particle) and the current solution. And dynamically adjusted using random weight d: (twenty three), in, It is a random number in the range [0,1]; the dynamic weight z= a and b are parameters for the spiral flight search; e represents the natural constant; z represents a random number or random variable, usually used to introduce randomness. This represents the updated position of an individual after the projection operation at iteration t+1; The random coefficient is usually a random number in the interval [0,1]. This dynamic adjustment mechanism enables LFGP to achieve a smooth transition from global exploration to local convergence. Steps 6-11: During the optimization process, based on the dynamic changes in fitness values ​​and the progress of the optimization objective, the residual guided projection (RGP), dual random projection (DRP), weighted random projection update (WRPU), and Levy flight guided projection (LFGP) are executed iteratively until the termination condition of the algorithm is met (stopping when the objective function value or parameter change is less than the specified precision). During the iteration, the optimal solution is gradually approached by dynamically adjusting the random factor, gradient information, and Jacobian matrix, and finally the improved projection iterative optimization algorithm TPIMO is obtained. Step 7 includes the following steps: Step 7-1: Divide the device degradation feature data obtained in Step 5 into a learning sample set and a validation set in a 2:1 ratio, and use the dynamic evolution feature set as the test set. Step 7-2: Select a portion of data from the learning sample set to initialize the DAVW-ELM dynamic adaptive weighted extreme learning model for remaining useful life prediction. Use the remaining data from the learning sample set to drive the DAVW-ELM dynamic adaptive weighted extreme learning model for remaining useful life prediction to update and learn. After learning is complete, use the validation set to train the DAVW-ELM model for remaining useful life prediction, and obtain the trained DAVW-ELM dynamic adaptive weighted extreme learning model for remaining useful life prediction, and obtain the preliminary prediction results. Step 7-3: Initialize the parameters of the improved projection iterative optimization algorithm TPIMO, including population size, individual position, dimension, maximum number of iterations, and search space boundary. Step 7-4: The prediction bias is used as the objective function by the improved projection iterative optimization algorithm TPIMO, and the individual fitness value is calculated according to the objective function; the weights and biases of the dynamic adaptive weighted extreme learning machine DAVW-ELM for the remaining useful life prediction model are optimized. Step 7-5: Determine if the maximum number of iterations has been reached. If it has, output the optimal weights and biases; otherwise, return to step 7-4. Step 7-6: Input the test set data into the Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM) model, which is optimized by the improved projection iterative optimization algorithm TPIMO, to make predictions and obtain the final remaining useful life prediction results. The device degradation characteristic data include single cell voltage, voltage fluctuation amplitude and frequency, recovery voltage magnitude and recovery time, current density, limiting current, cell internal temperature, temperature gradient, reaction gas pressure, electrochemical impedance spectroscopy, membrane performance parameters, catalyst activity, and hydrogen purity.

[0015] The present invention also provides a system for implementing the method, comprising: The data acquisition module is used to continuously collect monitoring data of the charging and discharging process of the fuel cell; The feature selection module is used to perform step 2; The data decomposition module is used to execute step 3; The data reconstruction module is used to execute step 4; The module builds the module used to execute step 5; An improved PIMO (Projection Iterative Optimization) module is used to execute step 6; The prediction module is used to perform step 7.

[0016] The present invention also provides an apparatus including a processor and a memory, the memory storing program code that, when executed by the processor, causes the processor to perform the steps of the method.

[0017] This invention has the following beneficial effects: 1. It uses the Consensus Feature Selection (CFS) method to select features from multivariate data, employs the Boruta secondary "shadow feature" mechanism to eliminate redundancy, significantly reduces the input dimensionality, and introduces the MCTS-XGBoost algorithm, which is an improvement based on Monte Carlo tree search, for secondary selection, ensuring that the selected variables have both high interpretability and high discriminative power; 2. It uses the DTW-RLMD robust local mean decomposition method, which is an improvement based on DTW, and utilizes the DTW distance to measure the similarity between submodals, achieving "high similarity clustering and low similarity stripping," preserving the degenerate main trend while compressing the data volume; 3. It uses the entropy-controlled principal component aggregation method for aggregation and reconstruction, first using the KPCA method improved by introducing an augmented matrix to filter nonlinear noise, and then combining the singular spectrum entropy to measure "information content" to automatically determine the weight of each principal component, achieving secondary noise reduction and compression, effectively suppressing error accumulation in long sequence prediction, and providing high signal-to-noise ratio input for the deep learning stage; 4. The TST-former lifetime prediction model uses a "seasonal-trend" dual-channel + The gating mechanism is specifically designed to capture the multi-scale periodic components and long-term drift of fuel cell performance degradation. It automatically estimates the main period using FFT, avoiding manual window setting, improving model versatility, and outputting multi-scale health factors to provide a highly discriminative dynamic evolutionary feature set for the final RUL regression. 5. The Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM) introduces sample importance weighting to enhance the fitting of key degradation stages and significantly reduce later prediction drift. 6. An improved projection iterative optimization algorithm (TPIMO) incorporates Logistic chaotic mapping and spiral flight search initialization to improve population diversity and global exploration efficiency. The four projection strategies—residual guidance, double randomization, weighted randomization, and Lévy flight—work synergistically to improve the ability to escape local extrema in high-dimensional non-convex spaces. The projection operator consistently pulls the solution back into the feasible region, ensuring stable and rapid convergence of the optimization process, ultimately obtaining the optimal weights and biases of DAVW-ELM. Attached Figure Description

[0018] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments, and the advantages of the present invention in the above and / or other aspects will become clearer.

[0019] Figure 1 A schematic diagram of the TST-former structure, which is the remaining useful life prediction model provided by this invention.

[0020] Figure 2 A flowchart illustrating the improved projection iterative optimization algorithm provided by this invention.

[0021] Figure 3 The diagram illustrates the weighting process of the improved projection iterative optimization algorithm provided by this invention for optimizing the dynamic adaptive weighted extreme learning machine.

[0022] Figure 4 This is a schematic diagram of the overall process of the present invention. Detailed Implementation

[0023] like Figure 1 , Figure 2 , Figure 3 and Figure 4 As shown, this invention provides a method for predicting the remaining service life of a proton exchange membrane fuel cell. First, micro-sensors are deployed within the fuel cell stack to continuously collect multi-channel dynamic data from the fuel cell, constructing a multi-dimensional time-series dataset. A consensus feature selection (CFS) method is used to select features from the multivariate data. Second, the original fuel cell data is decomposed using a robust local mean decomposition (DW-RLMD) improved by DTW to obtain sub-modes, reducing the non-stationarity of the sequence. Entropy-controlled principal component aggregation is used to reconstruct multiple sub-sequences to reduce computational complexity. The life prediction model TST-former is used to extract battery aging features in stages, outputting multi-scale health factors as preliminary prediction results. Finally, the TST-former model is combined with other methods to predict the remaining service life of the fuel cell. The predicted values ​​are reconstructed into a dynamic evolutionary feature set, which serves as the input to the Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM) for predicting remaining battery life. A Logistic chaotic mapping is used to initialize the population of the projection iterative optimization algorithm. A spiral flight search strategy is introduced during the algorithm's position update phase, and the projection operator ensures that the solution always lies within the feasible region, resulting in an improved projection iterative optimization algorithm, TPIMO. Simultaneously, TPIMO is used to optimize the weights and biases of the DAVW-ELM, and the optimal weights, biases, and test data samples are input into the DAVW-ELM prediction model for fusion prediction. This invention solves the problem of inaccurate prediction of the remaining battery life and improves the prediction accuracy of battery remaining life. The specific steps are as follows: Step 2: Boruta Algorithm in the Consensus Feature Filtering Method (CFS): (2.1): Boruta first copies the original data feature matrix sample set R with m rows and n columns to obtain a copy sample. Then, it performs random row transformation on each column of the copy sample to obtain a shadow feature sample L with m rows and n columns. Then, it concatenates R and L horizontally to obtain a new matrix N=[R,L] with dimension m×2n, which is used as the input for the feature importance evaluation of the subsequent random forest. Finally, it outputs the optimal feature subset S. (2.2): Run a random forest classifier on the mixed samples N, mask individual features sequentially and record the average loss value (Meanlmp) of the model accuracy. Calculate the average importance (MeanImp) of all features using SHAP value theory, mark features with importance lower than Meanlmp as "unimportant", and permanently remove them from the feature set; (2.3): Extract the maximum value of the importance of the shadow feature Maxlmp, compare it with the importance score of the true feature. The true features with an importance greater than Maxlmp are marked as "important", and other features are marked as "tentative".

[0024] (2.4): Select the importance score of the best shadow feature as the initial screening threshold, gradually compare the importance of each true feature with the best shadow feature, retain the true feature with the higher score, that is, whether the feature has a higher score than the best shadow feature, repeatedly remove unimportant features, and continuously optimize the feature set. (2.5): The algorithm stops when there are no more “provisional” features or when the preset number of iterations for the random forest is reached.

[0025] Step 2: Using the MCTS-XGBoost algorithm, an improvement on Monte Carlo tree search in the consensus feature selection method (CFS), select highly relevant feature variables as input variables: (2.6): Initialize a search tree, with the root node representing the entire candidate feature subset S. Set the parameters of the Monte Carlo Tree Search (MCTS), such as the number of searches N, the maximum depth of the tree D, the number of samples per node K, and the balance parameter C between search and utilization; (2.7): XGBoost utilizes the feature selection capability of its tree structure to perform secondary evaluation and ranking of the candidate feature subset S (S∈R) output by Boruta. Using the feature subset S selected by Boruta as input, the feature weights are optimized through the built-in feature selection mechanism of the XGBoost tree structure. In each iteration, XGBoost adds a new regression tree to fit the residual of the previous prediction. Let the predicted value of the ensemble model with t-1 existing trees be... The newly built regression tree in round t The goal is to fit the current residual, and the iterative process is as follows: (1), in The whole represents the predicted value of the j-th sample in round 0 (initial state).

[0026] (2.8): The importance of each feature in the current feature set S is evaluated using the XGBoost feature importance evaluation mechanism. When splitting a node in each tree, the optimal feature and split node are selected based on the feature gain. The essence of the gain formula is to measure the degree of reduction in the model loss function before and after splitting. The gain formula is as follows: (2), in and These represent the gain or signal parameters associated with the left and right channels, respectively. For regularization parameters; These are the weighting coefficients; α is a very small positive number used to ensure numerical stability; α and β are adjustable power parameters. and These are the feature parameters associated with the left and right channels, respectively.

[0027] (2.9): Based on the feature importance evaluation results, the feature subset S is divided into multiple subspaces, and each subspace is used as a new node in the search tree; (2.10): A Bayesian optimized surrogate model is used for spatial partitioning, and the node value is evaluated based on the upper bound of confidence interval (UCB) algorithm. The upper confidence bound UCB value is calculated, and the node with the highest UCB value is selected as the next exploration target. (2.11): Within the feature subspace corresponding to the selected node, random sampling is performed again to generate a new feature set. The search path is dynamically determined by iteratively calculating the UCB value of each node to obtain the optimal parameter configuration. The definition of UCB is: (3), K represents the degree of exploration during the search. When C=0, it is a greedy strategy, which is prone to entering local optima. The number of sampling points for the parent node; This represents the number of sampling points for the current node. This represents the average value currently observed in the i-th feature subspace of node j; This represents the number of times that node j has been visited or sampled in the i-th feature subspace.

[0028] (2.12): Select the leaf node with the highest UCB value from the search tree as the final feature set, and output the feature set corresponding to the node as the final result of XGBoost feature selection.

[0029] Step 3: Decompose the fuel cell charge and discharge process monitoring data using DTW-based robust local mean decomposition (DW-RLMD): (3.1): The RLMD (Robust Local Mean Decomposition) method is a signal processing method based on the Local Mean Decomposition (LMD) method, which optimizes boundary conditions, envelope estimation, and termination criteria. It preprocesses the original signal and extracts the solution signal. : (4), in It is a discrete signal; It is a local average function; It is a local envelope function; It is RLMD for the original signal Preprocessing operations are performed to extract target signals that are easier to decompose. ; (3.2): During the LMD iteration process, if the frequency of the frequency modulation signal in two consecutive iterations satisfies w 1k+2 >w 1k+1 And w 1k+3 >w 1k+2 When the iteration stops, the result of the nth iteration is output, obtaining the pure frequency modulated signal Y2(a) and the envelope signal h2(a), which are multiplied to obtain the product function Z. PF (a). The discrete signal X(a) is decomposed into f PF components and one residual signal. ; Among them, w 1k+2 The frequency value of the frequency-modulated signal is calculated with an iteration step size of 1k+2.

[0030] (3.3): Given two sets of PF sequences, respectively and Construct a dimension as Cost matrix Define the distance matrix ; in, express and of Norm (Euclidean distance); Represents a sequence The m-th element in Represents a sequence The nth element in; Represents the i-th element of sequence V eigenvectors; Represents the i-th element of sequence h eigenvectors.

[0031] (3.4): Construct a cumulative cost matrix , used in the cost matrix Calculate the minimum cumulative distance between two sets of PF sequences from the starting point (1,1) to the ending point (m,n), and define the weighted sum of the local cost measures as r(i, j); (3.5): Cost Matrix The value of the bottom right element r(m,n) is the DTW distance between PF sequences. The smaller the distance, the greater the similarity between the sequences. (3.6): First, use DTW distance to screen out PF sequence clusters with high similarity; reconstruct a group of PF components from the sequences within the cluster. Then, iteratively peel off: as long as the DTW distance between a component and the other components is less than the set threshold range [40,60], it is recorded as a new independent component PF1 and extracted. Then repeat the operation in the remaining components until all PF components are separated into unrelated subsequences.

[0032] Step 4: Reconstruct multiple subsequences using entropy-controlled principal component aggregation: (4.1): The entropy-controlled principal component aggregation method introduces an augmented matrix improved kernel principal component analysis (KPCA) to process the high-dimensional correlation of the trajectory matrix, and then integrates singular spectral entropy to screen periodic or deterministic subsequences; (4.2): Construct the trajectory matrix, assuming the original time series is... The length is N. A window length L is selected, and the one-dimensional sequence is transformed into a trajectory matrix X of L◊K using a sliding window. Through singular value decomposition and entropy analysis, the data is further aggregated, noise and non-primary structures are removed, and representative subsequences are obtained. (5), (4.3): The sample points x form the input trajectory matrix X, which contains N samples. A nonlinear mapping is used. Mapping sample point x to a high-dimensional space A new matrix is ​​obtained. ; (4.4): will exist KPCA dimensionality reduction is performed to obtain the feature space. The covariance matrix. The matrix is ​​calculated using a kernel function. Then, calculate the eigenvector α corresponding to the larger eigenvalue, and then obtain the corresponding weight vector. Sort them in descending order according to the proportion of the eigenvalues. The larger the proportion, the more important it is. Find the principal components in order of importance. (4.5): Construct an augmented matrix to introduce nonlinear changes, thereby expanding the feature space of the original data to improve the feature representation capability of KPCA. as follows: (6), In the formula: Let represent the augmented matrix, and q represent the time delay length.

[0033] (4.6): Construct a segmented trajectory matrix, divide the dimensionality-reduced matrix into local sub-blocks, and process them separately; (4.7): Perform singular value decomposition on the trajectory matrix X; (4.8): Normalize the singular values, calculate the proportion of each singular value to the sum of all singular values, and obtain the proportion of each singular value to the total energy, which serves as the probability distribution. ; (4.9): Calculate the singular spectral entropy. The larger the entropy value, the more uniform the distribution of singular values, and the higher the signal complexity (such as the presence of noise or nonlinear components); the smaller the entropy value, the more concentrated the signal energy is in a few singular values, indicating the existence of significant periodicity or deterministic structure. Based on the results of the entropy analysis, several representative subsequences are reconstructed.

[0034] (4.10): For each low-entropy submatrix Before selection Principal singular value reconstruction: (7), in, Represents the original low-entropy submatrix Approximate or reconstructed value; Let j be the singular value of the i-th submatrix; This represents the j-th left singular vector of the i-th submatrix; Let represent the j-th right singular vector of the i-th submatrix.

[0035] Step 5: Use the life prediction model TST-former to extract battery aging features in stages, and output the obtained multi-scale health factors as the first-level prediction results (equipment degradation feature data): Step 5-a1: Perform a Fast Fourier Transform on the aggregated and reconstructed data (load data). By utilizing the periodicity and symmetry of the signal to reduce the amount of computation, the time-domain signal is converted to the frequency domain. According to Passevar's theorem, calculate the energy spectral density in the frequency domain. After excluding the DC component, find the frequency index corresponding to the maximum energy. Convert the frequency index into a period to obtain the main period M of the load data.

[0036] Step 5-a2: Aggregate and reconstruct the sequence (load time series). The data is transformed into a two-dimensional trajectory matrix H; a window length W is selected, and the load data is intercepted and segmented into two or more segmented signals according to the window length W, and the segmented signals are embedded into the matrix to obtain the trajectory matrix H: (8), Wherein, the number of data segments K = l - W + 1 and K > W; This represents the k-th segment signal vector, which is a subsequence of length W extracted from the original sequence by sliding window. This represents the l-th data point in the original sequence.

[0037] Step 5-a3: Perform singular value decomposition on the trajectory matrix H, decomposing it into the product of three matrices; calculate the contribution rate r of the singular spectral components, and extract noise components to reduce noise interference. The trajectory matrix X is further optimized after this processing. The contribution rate r and the filtered trajectory matrix... The formula is: (9), (10) in, It is a singular value. For the filtered first A grouped matrix; d is the upper limit of the effective rank, an intermediate parameter.

[0038] Step 5-a4, group each matrix Transform into a feature subsequence of length W ,in The w-th element in the representation; each feature subsequence represents a certain feature of the load data, such as long-term trend, seasonal trend, etc. Step 5-a5: Embed the time series, decomposing the time series into trend A and seasonal B components along the channel dimension; embed the seasonal components into E. B The trend is embedded into E by feeding the attention module and the feedforward network FFN through a subtraction mechanism. A Feed to interactive module A which monitors seasonal information. s and F s Finally, the seasonality and trend are combined to obtain the final prediction (focus module, feedforward network FFN, interaction module A). s and F s It is a component of the TST-former model. Step 5-a6, decompose the sequence data X A X BMapping from the original space to a new space (the TST-former model undergoes FFT-SSD decomposition; the sequence reconstructed by entropy-controlled principal component aggregation is input into the TST-former model for FFT-SSD decomposition to obtain X). A X B X A ∈R, X B ∈R, where R represents the set of real numbers; using direct linear layers and dropout layers to create a global covariate X. mark Trends and seasonal embeddings: (11), (12) Concat represents a connection operation. Indicates trend item embedding, Indicates seasonal item embedding; This represents a linear layer (fully connected layer) used to map the concatenated high-dimensional vector to the embedding space; Concat() uses a separate linear layer to contain... It is used on datasets containing information and different components.

[0039] Step 5-a7: Since the seasonal component without a moving average kernel can better highlight the intrinsic characteristics of time series data, the seasonal embedding E is added. s Feeding is provided to the attention and feedforward network FFN; Step 5-a8: Considering that untrainable moving average kernels can lead to unreliable trend patterns, seasonal information is fused through the Interactive Module (IM) to assist in learning trend branches (the Interactive Module is a component of the TST-former model); focus on weighted A S It effectively reflects the dependencies between multiple variables and can be transformed into a coefficient matrix to update the trend embedding; the signals discarded by the seasonal branch can be regarded as meaningful information to guide the representation of the trend embedding. (13) (14); in, This indicates that the trend embedding obtained in the previous step is the input representation of the trend branch; It is an activation function that compresses the output to the [0,1] interval, representing the weight or gating signal; This represents a 1×1 convolution, used for dimensionality reduction or channel transformation, mapping the output of multiple convolutions to a suitable dimension; Indicates multi-level convolution operations; This represents the fused trend representation, which combines seasonal information, trend embedding, and enhanced weights. This indicates a splicing operation.

[0040] Step 5-a9: At the end of the seasonal and trend output blocks, a gate mechanism σ is designed for the two streams to autonomously regulate the speed of information transmission; the gate mechanism for the seasonal and trend streams is represented as follows: (15) (16) in, Indicates the seasonal item output; Indicates the trend term output; , , , These represent different 1×1 convolutional layers used to extract or transform features; The intermediate feature representation of the seasonal branch (which may be the result of seasonal embedding after a certain layer of processing).

[0041] Step 5-a10: Use the output of the previous TwinsBlock design overlay module as the input of the next TwinsBlock, overlay N TwinsBlocks to learn seasonal and trend representations, and then add the N TwinsBlocks together through linear projection to obtain the final prediction result, i.e., equipment degradation feature data.

[0042] Step 5: Combine all predicted values ​​from the TST-former to reconstruct a dynamic evolutionary feature set as input to the DAVW-ELM dynamic adaptive weighted extreme learning machine for remaining useful life prediction model. All predicted values ​​from the combined TST-former are reconstructed into a dynamically evolving feature set, which serves as the input to the Dynamically Adaptive Weighted Extreme Learning Machine (DAVW-ELM) for remaining lifespan prediction. DAVW-ELM employs a dynamic adjustment strategy, adjusting the input layer weights based on the output of each hidden layer neuron. This mainly involves calculating the activation level of each hidden neuron i and each training sample; calculating the variance of the activation level of each hidden neuron (a larger variance indicates a stronger ability to distinguish between different samples); and adjusting the input layer weights based on the activation level variance. (5.12): Different weight update strategies have different impacts on model performance. DAWW-ELM adopts an adaptive weight update strategy, dynamically selecting the appropriate update strategy based on the model's training state. It monitors the trend of training error changes; if the error decreases slowly, it indicates that the current strategy may be inefficient. Multiple weight update strategies can be selected, and the appropriate strategy can be dynamically chosen based on the trend of training error changes. (5.13): DAVW-ELM introduces a variable weighting mechanism to treat different samples differently based on their importance. Important samples are given greater weights to ensure that the model can learn the features of these samples more accurately.

[0043] Step 6: Improve the projection iterative optimization algorithm PIMO: (6.1): First, set the population size, dimension, number of iterations, and upper and lower bounds of the search space for the PIMO algorithm; (6.2): ​​The PIMO algorithm first initiates the optimization process by improving the initial population through a Logistic chaotic mapping. Given the population size S and the variable dimension G, each projective agent is considered as a search agent for the algorithm, and the position of each agent is: (17) in, This represents the k-th position of the i-th individual in the j-th dimension. and These represent the upper and lower bounds of the search space in the j-th dimension, respectively. This represents the initial Logistic value corresponding to the i-th individual, j-th dimension, and k-th position.

[0044] (6.3): Calculate the fitness value of all search agents; (6.4): Determine if the fitness value satisfies the termination objective (if the absolute error between the current optimal search agent's fitness value and the theoretical optimal fitness value for this optimization problem is less than 10%). If the termination objective is met, the current optimal solution is output. If the objective is met, the current optimal search agent is returned as the final solution. If the objective is not met, the subsequent optimization steps are executed until the optimal fitness value is reached, and the optimal parameter configuration is output. (6.5): In the Residual Guided Projection (RGP) stage, the algorithm repeatedly projects in the solution space. The algorithm randomly selects two individuals with the lowest residual fitness as guiding agents X. v1 and X v2 And calculate the gradient based on them. The formula for calculating the gradient W is: (18) in, This is the best position at present; , It is a random number between [0,1]; C is a parameter used to adjust the weights of different parts when calculating the gradient, which can control the balance between exploration and utilization in the search process.

[0045] (6.6): The residual-guided projection (RGP) stage introduces a Jacobian matrix J, where the element in the a-th row and b-th column of the Jacobian matrix J is... Represent the components of the objective function Relative to variables The partial derivatives are then used. Subsequently, the gradient and the Jacobian matrix J are incorporated into the projection step to complete the position correction. (6.7): Double random projection process: The algorithm first randomly selects two indices V3 and V4 from the particle swarm as reference points for projection update. When the gradient update is chosen by random decision, the algorithm calculates the update direction vector k1 based on the difference between the current position and the current optimal position: (19) If we choose to use Jacobian matrix projection for updating, the update direction vector k2 is calculated as follows: (20) (twenty one), in, This is the new location after the DRP update; This indicates the step size for each iteration, controlling the "movement range" for each iteration; This represents the original position (current position) of the i-th individual before the update.

[0046] (6.8): Weighted random projection update (WRPU) guides particle position updates through random factors r7 and r8 and adaptive factor α, and dynamically selects random weighted projection and adaptive correction projection paths; (6.9): The Lévy Flight Guided Projection (LFGP) process achieves efficient global search and local optimization by simulating the random motion of organisms in nature. The trigger probability P of the Lévy Flight Projection is determined by the following formula: (twenty two), Where t represents the current iteration number and U represents the maximum iteration number.

[0047] (6.10): The Levy flight position update formula introduces a spiral flight search strategy, which combines the optimal solution. (The best particle obtained in the first three steps) and the current solution And dynamically adjusted using random weight d: (twenty three), in It is a random number in the range [0,1]; dynamic weights a and b are parameters for the spiral flight search; e represents the natural constant; z represents a random number or random variable, usually used to introduce randomness. This represents the updated position of an individual after the projection operation at iteration t+1; The random coefficient is usually a random number in the interval [0,1]. This dynamic adjustment mechanism enables LFGP to achieve a smooth transition from global exploration to local convergence.

[0048] (6.11): During the optimization process, based on the dynamic changes in fitness values ​​and the progress of the optimization objective, the Residual Guided Projection (RGP), Double Random Projection (DRP), and Weighted Random Projection Update (WRPU) steps, as well as the Lévy Flight Guided Projection (LFGP), are executed iteratively until the termination condition of the algorithm is met. During the iteration, the optimal solution is gradually approximated by dynamically adjusting the random factor, gradient information, and Jacobian matrix, ultimately yielding the improved projection iterative optimization algorithm TPIMO.

[0049] Step 8: Evaluate the accuracy of the current model using three metrics: Root Mean Square Error (RMSE), Mean Absolute Percentage Error (MAPE), and Coefficient of Determination (R²). The formulas for RMSE, MAPE, and R² are as follows: (twenty four), (25), (26) Where RMSE is the root mean square error, MAPE is the mean absolute percentage error, and R² is the coefficient of determination. For predicted values, For the actual values ​​of the training samples, is the average value, and n is the number of samples.

[0050] Compared to the traditional RNN+LSTM method for predicting the remaining service life of proton exchange membrane fuel cells (PEMFCs), this patent proposes a "TST-former+DAVW-ELM dual prediction" mode, achieving a significant improvement in prediction accuracy, response speed, and energy consumption control. The PEMFC (rated capacity 2.8Ah, operating voltage 3.0-4.2V) was subjected to cyclic aging tests at 25℃. The charge / discharge regime was 1C constant current charging to 4.2V, followed by constant voltage charging to ≤0.05C, and then 1C constant current discharging to 3.0V. The first 30 samples were selected as the training set, and the last 20 as the test set. The end-of-life (EOL) was defined as the capacity decay to 80% of the rated capacity (2.24Ah). The goal is to predict the remaining service life (RUL) based on early cycle data. The TST-former extracts multi-scale health factors in stages, dividing the entire battery lifecycle into an early aging stage (0-200 cycles), a mid-term stable stage (201-600 cycles), and a late-stage accelerated aging stage (601 cycles to EOL). TST-former captures the degradation characteristics of each stage through an attention mechanism at different time scales (short-term window = 50 cycles, mid-term window = 200 cycles, long-term window = all stages). The predicted values ​​of TST-former for each stage (including time-series predictions for 9 high-frequency (HF) cycles) are combined with key dynamic indicators from the original features (such as changes in the slope of the voltage curve and fluctuations in charge / discharge efficiency) to construct a dynamic evolution feature set. A DAVW-ELM model is constructed, dynamically and adaptively weighting each feature according to its importance in different aging stages (weights are dynamically adjusted based on the SHAP values ​​of the training set data; for example, the weight of HF3 in the late stage is increased to 0.3, while the weight of HF1 in the early stage is 0.2). In the final prediction accuracy comparison, the TST-former+DAVW-ELM combined model has a mean absolute error (MAE) of 18.7, which is lower than both traditional ELM and LSTM, with a root mean square error (RMSE) of 25.3 and a coefficient of determination (R²) of 0.96. Compared to LSTM, the MAE is reduced by 42.5%, and the RMSE is improved by 69.9%, attributed to the multi-scale HF model capturing stage-specific degradation patterns. Simulation data shows that TST-former reduces the battery pack temperature prediction error RMSE from 3.15... Decreased to 0.92 The accuracy was improved by 69.9%; the MAPE accuracy of variable operating condition temperature control decreased from 10.80% to 3.50%, an improvement of 67.6%. DAVW-ELM utilizes TPIMO to optimize hyperparameters in real time, reducing the response time from 15-20 seconds to 3-5 seconds, and the temperature overshoot was reduced from... 5.2 Narrow to 1.3 The energy consumption of the cooling system was reduced by 75.0%; at the same time, the energy consumption of the cooling system was reduced by 35.2%, the battery cycle life increased from 1200 cycles to 1850 cycles, an improvement of 54.2%; the thermal runaway early warning response was improved by 76.4%, and the system updates model parameters online every 12 hours without the need for downtime retraining. The results are shown in Table 1.

[0051] Table 1 Comparison of Optimization, Improvement and Enhancement

[0052]

[0053] This embodiment also provides a health management system for lithium iron phosphate batteries in electric vehicles, including: 1. Application Scenario: A new energy vehicle manufacturer faces the problem of inconsistent battery degradation caused by high temperature, high frequency fluctuations during fast charging, and complex operating conditions in the actual operation of lithium iron phosphate batteries (capacity 100Ah, nominal voltage 3.2V, target cycle life 2000 cycles) for pure electric passenger vehicles. It needs to realize battery health status (SOH) early warning and life optimization through real-time RUL prediction.

[0054] 2. Test subjects: 100 mass-produced vehicles equipped with this battery were selected, covering typical scenarios such as urban commuting (average daily mileage of 50km), long-distance highway driving (average daily mileage of 200km), and low-temperature environments (-10~5℃).

[0055] 3. Deployment: (1) Data acquisition points: High-precision sensors are integrated in the battery management system (BMS) to collect core parameters in real time. Electrical parameters: cell voltage (sampling accuracy ±0.5mV), total current (±1A), charging / discharging power (10Hz sampling); environmental and status parameters: cell temperature (two NTC sensors are arranged in each string of cells, range -40~85℃), battery pack humidity (5%~95%RH), cumulative cycle number; aging characteristic parameters: DC internal resistance (DCR, offline detection once a week by 1kHz AC impedance method), capacity decay rate (calibrated by standard charge and discharge test every 50 cycles).

[0056] (2) Edge computing hardware: The NVIDIA Jetson Xavier NX embedded platform (21TOPS computing power) is used, equipped with a 6-core ARM CPU and a deep learning acceleration unit, which supports real-time feature extraction and model inference. The TST-former degradation feature extraction model (1.3M parameters / 2.9ms inference latency) and the DAVW-ELM remaining lifetime prediction model (0.4M parameters / 1.7ms inference latency) are run. The total inference time of the dual-model cascade is 5.2ms, which meets the 10Hz real-time closed-loop requirement of fuel cell health management.

[0057] (3) Experimental Operation: The test period covered high temperature in summer and low temperature in winter, and a total of 120 million operational data points were collected. High temperature fast charging scenario: A vehicle was fast charged 3 times in an ambient temperature of 42℃ (charging current fluctuated from 1C to 2C to 1.5C). The model predicted the RUL change in real time: Initial SOH: 92% (RUL = 1640 cycles); After 3 fast charges, the model predicted that the SOH would decay to 90.5% and the RUL would shorten to 1580 cycles, with an error of only 0.2% compared with the offline capacity test result (SOH 90.3%). Low temperature decay scenario: A vehicle drove an average of 150km per day in an environment of -5℃ (frequent rapid acceleration / deceleration). The model predicted that the RUL decay rate was 1.8 times faster than that under normal temperature conditions through temperature-current coupling characteristics, triggering a level 2 warning in advance and suggesting that the user enable the battery preheating function. The actual RUL was extended by 120 cycles: RUL prediction error (RMSE) ≤ 5 cycles (corresponding to capacity error ≤ 0.5Ah), SOH estimation error ≤ 1.2%.

[0058] This invention provides a method, system, and apparatus for predicting the remaining service life of a proton exchange membrane fuel cell. Many methods and approaches exist for implementing this technical solution; the above description is merely a preferred embodiment of the invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of this invention, and these improvements and modifications should also be considered within the scope of protection of this invention. All components not explicitly stated in this embodiment can be implemented using existing technologies.

Claims

1. A method for predicting the remaining service life of a proton exchange membrane fuel cell, characterized in that, Includes the following steps: Step 1: Deploy micro-sensors within the fuel cell stack to continuously collect monitoring data on the charging and discharging process of the fuel cell; Step 2: The consensus feature selection method (CFS) is used to select features from the charging and discharging process monitoring data. Feature variables with high relevance are selected as input variables. The consensus feature selection method (CFS) includes: first, constructing shadow features from the charging and discharging process monitoring data using the Boruta feature selection algorithm for preliminary importance assessment, and then performing secondary screening using an extreme gradient boosting algorithm based on Monte Carlo tree search, retaining only the features that are recognized by both secondary screenings. Step 3: Decompose the fuel cell charging and discharging process monitoring data into two or more subsequences using the robust local mean decomposition based on dynamic time warping (DTW). Step 4: Reconstruct the subsequence using the entropy-controlled principal component aggregation method; Step 5: Establish the lifespan prediction model TST-former, use the lifespan prediction model TST-former to extract battery aging features in stages, and output multi-scale health factors as the first layer prediction results, i.e. equipment degradation feature data; establish the remaining lifespan prediction model Dynamic Adaptive Weighted Extreme Learning Machine DAVW-ELM, combine all the prediction values ​​of TST-former to reconstruct a dynamic evolution feature set, and use it as the input of the remaining lifespan prediction model Dynamic Adaptive Weighted Extreme Learning Machine DAVW-ELM; Step 6, improve the projection iterative optimization algorithm PIMO: use Logistic chaotic mapping to initialize the population of the projection iterative optimization algorithm; introduce a spiral flight search strategy in the algorithm position update stage, and ensure that the solution is always within the feasible region through the projection operator, thus obtaining the improved projection iterative optimization algorithm TPIMO; Step 7: The weights and biases of the dynamic adaptive weighted extreme learning machine DAVW-ELM for predicting the remaining useful life are optimized using the improved projection iterative optimization algorithm TPIMO to obtain the optimal weights and biases. The optimized dynamic adaptive weighted extreme learning machine DAVW-ELM is then used to fuse and predict the test data samples of the proton exchange membrane fuel cell to obtain the final prediction result.

2. The method according to claim 1, characterized in that, Step 1 includes: A 3×3 grid-like micro-sensor array is deployed in different areas of the fuel cell stack to simultaneously collect 12 data channels, including voltage, temperature, humidity, and gas pressure. A humidity sensor, a gas flow meter, and a membrane impedance monitoring module are also deployed to continuously collect monitoring data of the fuel cell's charging and discharging process. A dynamic adaptive sampling strategy is adopted, with triggered sampling, a threshold set, high-frequency sampling, and low power consumption mode maintained during normal operation. Ultimately, new data on fuel cell charging and discharging voltage, current, current density, temperature, capacity, internal resistance, cycle count, aging time, and output power are obtained.

3. The method according to claim 2, characterized in that, Step 2, the Boruta feature selection algorithm includes the following steps: Step 2-a1: Boruta first copies the original data feature matrix sample set R with m rows and n columns to obtain a copy sample. Then, it performs a random row transformation on each column of the copy sample to obtain a shadow feature sample L with m rows and n columns. Finally, it concatenates R and L horizontally to obtain a new matrix N=[R,L] with dimensions m×2n, which serves as the input for the feature importance evaluation of the subsequent random forest. Finally, it outputs the optimal feature subset S. Step 2-a2: Run a random forest classifier on matrix N, mask individual features in sequence and record the average loss value of model accuracy Meanlmp; calculate the average importance MeanImp of all features using SHAP value theory, mark features with importance lower than Meanlmp as unimportant, and permanently delete them from matrix N. Step 2-a3: Extract the maximum value of the shadow feature importance, Maxlmp, and compare it with the importance score of the true feature. Among the true features, those with an importance greater than Maxlmp are marked as important, and other features are marked as tentative. Step 2-a4: Select the importance score of the optimal shadow feature as the initial screening threshold, gradually compare the importance of each true feature with the optimal shadow feature, retain the true features with higher scores than the optimal shadow feature, repeatedly remove features with lower scores than the optimal shadow feature, and continuously optimize the feature set. Steps 2-a5 continue until there are no more provisional features, or the preset number of iterations for the random forest is reached, at which point the algorithm stops.

4. The method according to claim 3, characterized in that, Step 2, the MCTS-XGBoost algorithm based on Monte Carlo tree search improvement includes the following steps: Step 2-b1: Initialize a search tree, with the root node representing the optimal feature subset S; set the parameters of the Monte Carlo Tree Search (MCTS), including the number of searches N, the maximum depth of the tree D, the number of samples per node K, and the balance parameter C between search and utilization. Step 2-b2: XGBoost utilizes the feature selection capability of its tree structure to perform a secondary evaluation and ranking of the optimal feature subset S. The optimal feature subset S is used as input, and the feature weights are optimized through the built-in feature selection mechanism of XGBoost's tree structure. In each iteration, XGBoost adds a new regression tree to fit the residual of the previous prediction. Let the predicted value of the ensemble model with t-1 trees be... The newly built regression tree in round t The goal is to fit the current residual, and the iterative process is as follows: （1）， in The whole represents the predicted value of the j-th sample in round 0; Step 2-b3: Using the XGBoost feature importance evaluation mechanism, the importance of each feature in the optimal feature subset S is evaluated. When splitting a node in each tree, the optimal feature and the splitting node are selected based on the feature gain (Gain). The gain formula is: （2）， in This indicates the gain or signal parameter associated with the left channel; Indicates the gain or signal parameter associated with the right channel; For regularization parameters; These are the weighting coefficients; It is a very small positive number; α and β are adjustable power parameters; and These are the feature parameters associated with the left channel and the feature parameters associated with the right channel, respectively. Step 2-b4: Based on the feature importance evaluation results, divide the optimal feature subset S into two or more subspaces, and use each subspace as a new node in the search tree; Step 2-b5: Use a Bayesian optimized surrogate model to partition the space, evaluate the node value based on the upper bound of the confidence interval UCB algorithm, calculate the upper bound of the confidence interval UCB value, and select the node with the highest upper bound of the confidence interval UCB value as the next exploration target. Step 2-b6: Within the feature subspace corresponding to the selected node, random sampling is performed again to generate a new feature set. By iteratively calculating the upper bound of the confidence interval (UCB) for each node, the search path is dynamically determined to obtain the optimal parameter configuration. The definition of the upper bound of the confidence interval (UCB) is: （3）， Where K represents the degree of exploration during the search, and C=0 is a greedy strategy; The number of sampling points for the parent node; This represents the number of sampling points for the current node. This represents the average value currently observed in the i-th feature subspace of node j; This represents the number of times that node j has been visited or sampled in the i-th feature subspace; Step 2-b7: Select the leaf node with the highest UCB value of the upper bound of the confidence interval from the search tree, and use it as the final feature set. Output the feature set corresponding to the leaf node with the highest UCB value of the upper bound of the confidence interval as the final result of XGBoost feature selection.

5. The method according to claim 4, characterized in that, Step 3 includes the following steps: Step 3-1: Preprocess the original signal and extract the solution signal. : （4）， in It is a discrete signal; It is a local average function; It is a local envelope function; Robust Local Mean Decomposition (RLMD) is applied to the original signal. The target signal that is easier to decompose after preprocessing; Step 3-2, during the LMD iteration process, if the frequency of the frequency-modulated signal in two consecutive iterations satisfies w 1k+2 >w 1k+1 And w 1k+3 >w 1k+2 When the iteration stops, the result of the nth iteration is output, and the pure frequency modulated signal Y2(a) and the envelope signal h2(a) are obtained. The product function Z is obtained by multiplying the pure frequency modulated signal Y2(a) and the envelope signal h2(a). PF (a) The discrete signal X(a) is decomposed into f product function PF components and 1 residual signal. ;w 1k+2 The frequency value of the frequency-modulated signal is calculated with an iteration step size of 1k+2. Step 3-3: Given two sets of PF sequences, which are sequences... and sequence Construct a dimension as Cost matrix Define the distance matrix ; in, express and of Norm; Represents a sequence The m-th element in Represents a sequence The nth element in; Represents the i-th element of sequence V eigenvectors; Represents the i-th element of sequence h eigenvectors; Steps 3-4: Construct a cumulative cost matrix , used in the cost matrix Calculate the minimum cumulative distance between two PF sequences h and V from the starting point (1,1) to the ending point (m,n), and define the weighted sum of the local cost measures as r(i, j); Steps 3-5, Cost Matrix The value of the bottom right element r(m,n) is the DTW distance between PF sequences; Steps 3-6: Use DTW distance to screen out PF sequence clusters with high similarity; reconstruct a group of PF components from the sequences within the cluster, and then iteratively peel them off: if the DTW distance between a component and the other components is less than a set threshold range, then record the component as a new independent component PF1 and extract it. Repeat the iterative stripping operation on the remaining components until all PF components are separated into unrelated subsequences.

6. The method according to claim 5, characterized in that, Step 4 includes the following steps: Step 4-1, the entropy-controlled principal component aggregation method introduces an augmented matrix improved kernel principal component analysis method to process the high-dimensional correlation of the trajectory matrix, and then fuses singular spectral entropy to screen periodic or deterministic subsequences; Step 4-2, construct the trajectory matrix, assuming the original time series is... The length is N; where This represents the Nth element in the original time series x; Choosing a window length L, the one-dimensional sequence is transformed into a trajectory matrix X of L◊K using a sliding window. Through singular value decomposition and entropy analysis, the data is further aggregated, noise and non-primary structures are removed, resulting in representative subsequences. （5）； Step 4-3: Assemble the sample points x into an input trajectory matrix X, where X contains N samples, and use a nonlinear mapping... Mapping sample point x to a high-dimensional space A new matrix is ​​obtained. ; Step 4-4, In higher dimensions We then perform dimensionality reduction using augmented matrix improved kernel principal component analysis (KPCA) to obtain the high-dimensional space. The covariance matrix is ​​calculated using a kernel function. Then calculate the matrix. The eigenvectors corresponding to the larger eigenvalues ​​are obtained, and then the corresponding weight vectors are obtained. These weight vectors are then sorted in descending order according to the proportion of the eigenvalues, and the principal components are found in order of importance. Steps 4-5 involve constructing an augmented matrix by introducing nonlinear variations to expand the feature space of the original data. Represented as: （6）， Where q represents the time delay length; Steps 4-6: Construct a segmented trajectory matrix, divide the dimensionality-reduced matrix into local sub-blocks, and process them separately; Steps 4-7: Perform singular value decomposition on the trajectory matrix X; Steps 4-8 involve normalizing the singular values ​​and calculating the proportion of each singular value to the sum of all singular values. This yields the proportion of each singular value to the total energy, which serves as the probability distribution. ; Steps 4-9: Calculate the singular spectral entropy, and reconstruct two or more representative subsequences based on the entropy analysis results; Steps 4-10: For each low-entropy submatrix Before selection Principal singular value reconstruction: （7）， in, Represents the original low-entropy submatrix Approximate or reconstructed value; Let j be the singular value of the i-th submatrix; This represents the j-th left singular vector of the i-th submatrix; Let represent the j-th right singular vector of the i-th submatrix.

7. The method according to claim 6, characterized in that, Step 5, establishing the lifetime prediction model TST-former, includes the following steps: Step 5-a1: Perform a Fast Fourier Transform on the aggregated and reconstructed data. By utilizing the periodicity and symmetry of the signal, the computational load is reduced, and the time-domain signal is converted to the frequency domain. According to Passevar's theorem, calculate the energy spectral density in the frequency domain. After excluding the DC component, find the frequency index corresponding to the maximum energy, convert the frequency index into a period, and obtain the main period M of the load data. Step 5-a2, aggregate and reconstruct the sequence The data is transformed into a two-dimensional trajectory matrix H; a window length W is selected, and the load data is intercepted and segmented into two or more segmented signals according to the window length W, and the segmented signals are embedded into the matrix to obtain the trajectory matrix H: （8）， Wherein, the number of data segments K = l - W + 1 and K > W; This represents the k-th segment signal vector, which is a subsequence of length W extracted from the original sequence by a sliding window. This represents the l-th data point in the original sequence; Step 5-a3: Perform singular value decomposition on the trajectory matrix H, decomposing H into the product of three matrices; calculate the contribution rate r of the singular spectral components, and the selected trajectory matrix. The formula is: （9）， （10）， in, It is a singular value. For the filtered first A grouped matrix; d is the upper limit of the effective rank; Step 5-a4, group each matrix Transform into a feature subsequence of length W ,in This represents the w-th element; Step 5-a5: Embed the time series, decomposing the time series into trend A and seasonal B components along the channel dimension; embed the seasonal components into E. B The trend is embedded into E by feeding the attention module and the feedforward network FFN through a subtraction mechanism. A Feed to interactive module A which monitors seasonal information. s and F s Finally, the seasons and trends are combined to obtain the final prediction; Step 5-a6, decompose the sequence data X A X B Mapping from the original space to the new space yields X. A X B X A ∈R,X B ∈R, where R represents the set of real numbers; using direct linear layers and dropout layers to create a global covariate X. mark Trends and seasonal embeddings: （11）， （12）， Concat represents a connection operation. Indicates trend item embedding, Indicates seasonal item embedding; Indicates a linear layer; Step 5-a7, embed the season into E s Feeding is provided to the attention and feedforward network FFN; Step 5-a8 involves fusing seasonal information through the Interactive Module to aid in the learning of trend branches; signals discarded by seasonal branches are treated as meaningful information to guide the representation of trend embeddings. （13）， （14）； in, This indicates the trend embedding obtained in the previous step; It is an activation function; Represents a 1×1 convolution; Indicates multi-level convolution operations; This indicates the trend after fusion; Step 5-a9: At the end of the seasonal and trend output blocks, a gate mechanism σ is designed for the two streams to autonomously regulate the speed of information transmission; the gate mechanism for the seasonal and trend streams is represented as follows: （15）， （16）； in, Indicates the seasonal item output; Indicates the trend term output; , , , This represents different 1×1 convolutional layers; The intermediate features representing seasonal branches; Step 5-a10: Use the output of the previous TwinsBlock design overlay module as the input of the next TwinsBlock, overlay N TwinsBlocks to learn seasonal and trend representations, and then add the N TwinsBlocks together through linear projection to obtain the final prediction result, i.e., equipment degradation feature data.

8. The method according to claim 7, characterized in that, Step 5, establishing the dynamic adaptive weighted extreme learning machine (DAVW-ELM) for predicting remaining useful life, includes the following steps: Step 5-b1: Combine all prediction results of TST-former to reconstruct a dynamic evolution feature set, which serves as the input to the dynamic adaptive weighted extreme learning machine DAVW-ELM for the remaining useful life prediction model. The remaining useful life prediction model, Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM), employs a dynamic adjustment strategy. It adjusts the corresponding input layer weights based on the output of each hidden layer neuron, including: calculating the activation level for each hidden neuron i and each training sample; calculating the variance of the activation level for each hidden neuron; and adjusting the input layer weights based on the activation level variance. Step 5-b2: The Dynamic Adaptive Weighted Extreme Learning Machine (DAWW-ELM) for predicting remaining useful life adopts an adaptive weight update strategy. It dynamically selects an appropriate update strategy based on the training state of the DAWW-ELM and monitors the changing trend of the training error. Based on the changing trend of the training error, it dynamically selects the weight update strategy. Step 5-b3: The remaining useful life prediction model, Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM), introduces a variable weighting mechanism to differentiate between different samples based on their importance. Step 6 includes the following steps: Step 6-1: First, set the population size, dimension, number of iterations, and upper and lower limits of the search space for the Projection Iterative Optimization (PIMO) algorithm. Step 6-2: The Projective Iterative Optimization (PIMO) algorithm first initiates the optimization process by improving the initial population using a Logistic chaotic mapping. Given the population size S and the variable dimension G, each projective agent is considered a search agent for the algorithm, and the position of each agent is: （17）， in This represents the k-th position of the i-th individual in the j-th dimension; and These represent the upper and lower bounds of the search space in the j-th dimension, respectively. This represents the initial value of the Logistic regression corresponding to the i-th individual, the j-th dimension, and the k-th position. Step 6-3: Calculate the fitness value of all search agents; Step 6-4: Determine if the fitness value meets the termination objective. If it does, return the current optimal search agent as the final solution. If it does not, continue to execute subsequent optimization steps until the optimal fitness value is reached, and output the optimal parameter configuration. Step 6-5, Residual Guided Projection (RGP): The algorithm repeatedly projects in the solution space, and randomly selects two individuals with the lowest residual fitness as guiding agents X. v1 and X v2 And according to X v1 and X v2 Calculate the gradient. The calculation formula is: （18）， in, This is the best position at present; and It is a random number between [0,1]; C is a parameter used to adjust the weights of different parts when calculating the gradient; Step 6-6, the Residual Guided Projection (RGP) stage introduces the Jacobian matrix J, where the element in the a-th row and b-th column of the Jacobian matrix J... Represent the components of the objective function Relative to variables The partial derivatives are then used; subsequently, the gradient and the Jacobian matrix J are incorporated into the projection step to complete the position correction. Steps 6-7, Double Random Projection Process (DRP): The algorithm first randomly selects two indices V3 and V4 from the particle swarm as reference points for projection updates. When randomly selecting gradient updates, the algorithm calculates the update direction vector k1 based on the difference between the current position and the current optimal position. （19）， If we choose to use Jacobian matrix projection for updating, the update direction vector k2 is calculated as follows: （20）， （21）， in, It is the new position after the DRP update. This indicates the step size for each iteration; This represents the original position of the i-th individual before the update; Steps 6-8, Weighted Random Projection Update WRPU: Particle position updates are guided by random factors r7 and r8 and adaptive factor α, and random weighted projection and adaptive correction projection paths are dynamically selected; Steps 6-9, Lévy Flight Guided Projection (LFGP) process: This process achieves efficient global search and local optimization by simulating the random motion of organisms in nature. The trigger probability P of the Lévy flight projection is determined by the following formula: （22）； Where t represents the current iteration number and U represents the maximum iteration number; Steps 6-10 introduce a spiral flight search strategy into the Levy flight position update formula, combining the optimal solution. and the current solution And dynamically adjusted using random weight d: （23）， in, It is a random number in the range [0,1]; the dynamic weight z= a and b are parameters for the spiral flight search; e represents the natural constant. This represents the updated position of an individual after the projection operation at iteration t+1; Represents random coefficients; Steps 6-11: During the optimization process, based on the dynamic changes in fitness values ​​and the progress of the optimization objective, the residual guided projection (RGP), dual random projection (DRP), weighted random projection update (WRPU), and Levy flight guided projection (LFGP) are executed iteratively until the termination condition of the algorithm is met. During the iteration, the optimal solution is gradually approached by dynamically adjusting the random factor, gradient information, and Jacobian matrix, and finally the improved projection iterative optimization algorithm TPIMO is obtained. Step 7 includes the following steps: Step 7-1: Divide the device degradation feature data obtained in Step 5 into a learning sample set and a validation set, and use the dynamic evolution feature set as the test set; Step 7-2: Select a portion of data from the learning sample set to initialize the DAVW-ELM dynamic adaptive weighted extreme learning model for remaining useful life prediction. Use the remaining data from the learning sample set to drive the DAVW-ELM dynamic adaptive weighted extreme learning model for remaining useful life prediction to update and learn. After learning is complete, use the validation set to train the DAVW-ELM model for remaining useful life prediction, and obtain the trained DAVW-ELM dynamic adaptive weighted extreme learning model for remaining useful life prediction, and obtain the preliminary prediction results. Step 7-3: Initialize the parameters of the improved projection iterative optimization algorithm TPIMO, including population size, individual position, dimension, maximum number of iterations, and search space boundary. Step 7-4: The prediction bias is used as the objective function by the improved projection iterative optimization algorithm TPIMO, and the individual fitness value is calculated according to the objective function; the weights and biases of the dynamic adaptive weighted extreme learning machine DAVW-ELM for the remaining useful life prediction model are optimized. Step 7-5: Determine if the maximum number of iterations has been reached. If it has, output the optimal weights and biases; otherwise, return to step 7-4. Step 7-6: Input the test set data into the Dynamic Adaptive Weighted Extreme Learning Machine (DAVW-ELM) model, which is optimized by the improved projection iterative optimization algorithm TPIMO, to make predictions and obtain the final remaining useful life prediction results. The device degradation characteristic data include single cell voltage, voltage fluctuation amplitude and frequency, recovery voltage magnitude and recovery time, current density, limiting current, cell internal temperature, temperature gradient, reaction gas pressure, electrochemical impedance spectroscopy, membrane performance parameters, catalyst activity, and hydrogen purity.

9. A system for implementing the method according to any one of claims 1 to 8, characterized in that, include: The data acquisition module is used to continuously collect monitoring data of the charging and discharging process of the fuel cell; The feature selection module is used to perform step 2; The data decomposition module is used to execute step 3; The data reconstruction module is used to execute step 4; The module builds the module used to execute step 5; An improved PIMO (Projection Iterative Optimization) module is used to execute step 6; The prediction module is used to perform step 7.

10. An apparatus, characterized in that, It includes a processor and a memory, the memory storing program code that, when executed by the processor, causes the processor to perform the steps of the method as described in any one of claims 1 to 8.

Images