Pumped storage unit degradation trend prediction method and system based on polynomial XGBoost and double-layer prediction model

Abstract

The invention discloses a pumped storage unit degradation trend prediction method and system based on a polynomial XGBoost and a double-layer prediction model, and the method comprises the steps: building a health state evaluation model through employing the polynomial XGBoost for solving a problem of the interpretability of a health model, enhancing the fitting capability of the model for a nonlinear relation through the introduction of polynomial transformation, and achieving the prediction of the degradation trend of a pumped storage unit. Then, the contribution of each feature to model output is explained by adopting an SHAP method, so that the interpretability of the model is improved; in order to solve the possible problem that a single model is insufficient in PDI prediction trend capturing capacity, an improved IWUTP algorithm is adopted to optimize VMD parameters, and PDI signals are effectively decomposed. A double-layer prediction model is constructed, the GRU is used for capturing the time sequence characteristics of the PDI, then the GRU output is subjected to refined prediction through the IWUTP-KELM, and finally accurate prediction of the degradation trend of the pumped storage unit is achieved. According to the method, a good prediction effect is achieved on historical data, and the stability of the prediction model and the accuracy of the prediction effect are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of health management technology for pumped storage units, specifically to a method and system for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model. Background Technology

[0002] In existing technologies, the intermittency and uncertainty of renewable energy sources such as wind power and photovoltaics significantly increase the complexity of power system supply-demand balance and dispatch economics under high grid-connected conditions. Therefore, how to improve the controllability and absorption capacity of new energy sources, and ensure the stability and security of the power grid, has become a key issue that modern power systems urgently need to address.

[0003] Pumped storage, as the most mature energy storage technology with large-scale development potential, can effectively mitigate the volatility and intermittency of new energy sources such as wind and solar power through energy storage and release mechanisms, playing a crucial role in supporting large-scale grid integration of new energy sources and enhancing power system stability. However, with the increasing operating years of pumped storage units, structural fatigue and performance degradation gradually accumulate, posing potential challenges to their operational safety and system stability. Vibration trend prediction, as an important means of characterizing the evolution of unit operating status, has been widely used to assist in degradation analysis and early warning. However, relying solely on the changing trends of vibration signals is insufficient to fully reflect the dynamic evolution of the health level of pumped storage units, and its expressive power remains limited. In recent years, some studies have focused on state assessment based on health indicators, attempting to preliminarily classify equipment operating risks through monitoring signal characteristics. However, these methods mostly focus on the division of static state intervals, making it difficult to fully reflect the continuous degradation process of equipment performance over time. In contrast, time-series-based performance degradation trend modeling can deeply reveal the equipment degradation trajectory, improve the accuracy of predictions and the foresight of early warnings, and thus has become a research hotspot in the current field of health management.

[0004] However, existing methods for predicting the degradation trends of pumped storage units largely focus on constructing degradation indicators and modeling trends, paying less attention to the specific contributions of characteristic variables to the prediction results, leading to a lack of sufficient interpretability in the models. Furthermore, most methods for constructing Performance Degradation Indices (PDIs) rely on complete lifecycle data of the equipment from normal operation to failure. However, pumped storage equipment has a long lifespan and complex structure, making complete failure records difficult to obtain, resulting in a scarcity of real degradation data. This limits the application of lifecycle-based driving methods, and current degradation trend modeling often uses a single prediction model, which struggles to fully capture the complex degradation dynamics of the equipment, posing a challenge to the accuracy of the predictions. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides a method and system for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model. Regarding the interpretability of the health model, a health status assessment model is first constructed using polynomial XGBoost. The introduction of polynomial transformation enhances the model's ability to fit nonlinear relationships. Subsequently, the SHAP method is used to explain the contribution of each feature to the model output, improving the model's interpretability. To address the potential deficiency in capturing the PDI prediction trend of a single model, an improved IWUTP algorithm is used to optimize the VMD parameters, effectively decomposing the PDI signal. By constructing a two-layer prediction model, the GRU is first used to capture the temporal characteristics of the PDI, and then IWUTP-KELM is used to refine the prediction of the GRU output, ultimately achieving accurate prediction of the deterioration trend of the pumped storage unit.

[0006] The technical solution adopted in this invention is as follows:

[0007] A method for predicting the degradation trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model includes the following steps:

[0008] Step 1: Obtain the vibration signals and related operating characteristics of the pumped storage unit, construct a health status prediction model using the multinomial XGBoost algorithm, and introduce the SHAP method to evaluate the contribution of each input feature to the output of the health status prediction model, thereby realizing the interpretability analysis of the model.

[0009] Step 2: Input the vibration data of the deteriorating working conditions to be predicted into the trained health status prediction model to obtain the corresponding predicted vibration data. Construct a Predicted Displacement Indicator (PDI) based on the deviation between the predicted values ​​and the actual observed values ​​to quantify the degree of equipment deterioration;

[0010] Step 3: The constructed PDI sequence is decomposed using the IWUTP-optimized VMD method to obtain K intrinsic mode functions in order to obtain multi-scale temporal information of degradation features;

[0011] Step 4: Use GRU to predict each intrinsic mode function separately, and superimpose the prediction results of each mode at each time step to construct the corresponding prediction training set sequence and prediction test set sequence. Simultaneously, retain the actual training set sequence of PDI as a reference for subsequent modeling;

[0012] Step 5: Take the predicted training set sequence, predicted test set sequence, and actual training set sequence of PDI output from GRU as input to construct a KELM model based on IWUTP optimization, further model and predict the degradation trend, and finally output a more accurate degradation trend prediction result.

[0013] Step 1 includes the following steps:

[0014] S1.1: Acquire the original vibration signal of the pumped storage unit and key operating characteristics such as power generation head, active power, and guide vane opening;

[0015] S1.2: A health status prediction model is constructed based on the polynomial XGBoost algorithm. The calculation formulas included in the XGBoost algorithm are shown below:

[0016] XGBoost is a high-efficiency machine learning algorithm based on the gradient boosting framework. It improves prediction performance by progressively building an additive model and integrating multiple weak classification or regression decision trees. It optimizes model parameters by minimizing an objective function that includes a loss function and a regularization term, and supports parallel computation and various regularization strategies, effectively preventing overfitting and improving the model's generalization ability. Let the training samples be M = (x1, y1), (x2, y2), ..., (x... n ,y n ), i = 1, 2, 3, ..., n where x i Let y be the feature vector of the i-th sample. i For actual values, assume each decision tree is modeled. in, f(x) is the predicted value. i Let F represent the model's predicted output for the input feature x. Then, the overall objective function F is as follows:

[0017]

[0018] In formula (1): The loss function is represented by equation (2); Ω(f(x) t )) represents the regularization term, and its calculation formula is shown in equation (3); f(x) t ) indicates that the model adapts to the input feature x. t The predicted output, where t represents the time step and n represents the number of samples;

[0019] Using the squared error as the loss function, its expression is as follows:

[0020]

[0021] In formula (2): This represents the loss function, used to measure the error between the true value and the predicted value.

[0022] A regularization term is introduced to control model complexity; its expression is shown below:

[0023]

[0024] In equation (3): Ω is the model complexity penalty term; γ represents the regularization coefficient used to control the number of numbers; T is the number of regression trees; λ represents the L2 regularization coefficient of the leaf node weights; π represents the leaf node weight vector of the regression tree, which contains the feature weights composed of the predicted values ​​output by all leaf nodes in the model.

[0025] The Boosting algorithm employs an additive model approach, progressively stacking several weak regressors to construct a strong regression model. In the t-th iteration, the model's prediction is the sum of the predictions from the previous t-1 iterations and the prediction from the current t-th tree. Essentially, Boosting optimizes the prediction performance of the current weak classifier in the t-th iteration by fitting the residuals of the current model from the previous t-1 iterations, thereby minimizing the overall objective function. In constructing the health status prediction model, the hydroelectric head, active power, and guide vane opening are selected as the model input features, denoted as J = [j1, j2, j3], where j1, j2, and j3 represent the hydroelectric head, active power, and guide vane opening, respectively.

[0026] A polynomial extension method is introduced based on the original XGBoost features. Specifically, based on the original XGBoost features, a fourth-order polynomial extension is used to generate power terms and multidimensional interaction terms (such as square terms, cube terms, and products between features) to describe nonlinear degradation relationships and the interactive effects of multiple features; the specific formula is as follows:

[0027]

[0028] Where: J represents the original feature vector; a s Indicates the corresponding feature j s The polynomial power of ; φ(J) is the feature set after polynomial expansion; s represents the feature index, with values ​​ranging from 1, 2, 3; It is a set of non-negative integers.

[0029] S1.3: The SHAP method is introduced to evaluate the contribution of each input feature to the model output, thereby enabling interpretability analysis of the model. A summary of SHAP features is as follows: Figure 2 As shown. The formula for SHAP is as follows:

[0030] In a single decision tree, the SHAP value is calculated by analyzing the splitting of a feature across all decision paths. For a given feature, its SHAP value represents the expected marginal contribution of that feature to the predicted output across all possible paths, expressed as follows:

[0031]

[0032] Where, x i The feature is represented by the Shapley value; the Shapley value represents the feature x.i Contribution to the model output; Represents all cases containing x i The summation is performed on path p, where Pr(p) is the sampling probability or weight of path p; △f(p,x) i ) represents the feature x in path p. i This results in an increase in the model output.

[0033] Finally, the SHAP values ​​are summed to quantify the contribution of each feature to the overall prediction result of the model.

[0034]

[0035] In formula (6): SHAP j (g) represents feature x i global contribution, Shapley i In the nth sample, feature x i The Shapley value, where n represents the total number of samples.

[0036] Step 2 includes the following steps:

[0037] S2.1: Input the power generation head, active power, guide vane opening and vibration data under adverse conditions into the trained health state prediction model, including equations (1) to (4), and construct F(H(t), P(t), G(t)), where: H(t) is the power generation head, P(t) is the active power, and G(t) is the guide vane opening, to obtain the corresponding predicted vibration data;

[0038] S2.2: The PDI is constructed based on the deviation between the predicted value and the actual observed value. The calculation formula is as follows:

[0039]

[0040] Wherein: F t V represents the vibration value predicted by the health model at time t. t This refers to the actual measured value at the corresponding time. Step 3 includes the following steps:

[0041] S3.1: The WUTP algorithm models candidate solutions in the optimization problem as pressure gradient vectors. By simulating the water transport process from soil to leaves in plants, it uses water flow velocity to drive dynamic adjustment and iterative updates of candidate solutions in multiple stages. This process can be divided into three consecutive biomimetic stages: soil to root surface, root surface to xylem, and xylem to leaves, ultimately completing an efficient search from the initial solution to the global optimum. The specific initialization steps of the algorithm are as follows:

[0042] y i =[lb j+rand×(ub j -lb j (8)

[0043] Where: y i Represents the i-th individual; lb j with ub j represents the lower and upper bounds of the j-th dimension, respectively, and rand represents a random variable that follows a uniform distribution on the interval [0,1].

[0044] After initializing the population, the plant leaves and moving water particles were initialized separately. The initialization steps for the plant leaves are as follows:

[0045] Leaves = [l1, l2, ..., l i ,…,l m (9);

[0046] Where: Leaves represents the collection of plant leaves; m represents the problem dimension, l i Let i be the i-th variable. The water particles in the plant are initialized as follows:

[0047] Particles = [p1, p2, p3, ..., p n (10);

[0048] Where: Particles represents the set of moving water particles; n is the population size; the water particles must have sufficient hydrostatic pressure P, p1, p2, p3, ..., p n These represent individual water particles, driving their flow to a position close to the blade. At this point, the water particles are under a pressure gradient. Under the influence of [the flow], the roots travel a distance x from the soil to the leaves. Specifically, the plant roots travel a distance x from the soil... i =0 absorbs water particles, and these particles travel a distance x within the plant leaf. f The flow distance can be determined by Obtained through calculation.

[0049] In summary, the population can be represented as follows:

[0050]

[0051] Where: Totalop represents the set of the population; P 1,m P represents the hydrostatic pressure of the first water particle in the m-th dimension; 2,m P represents the hydrostatic pressure of the second water particle in the m-th dimension; n,m x represents the hydrostatic pressure of the nth water particle in the mth dimension; 1,m x represents the distance the first water particle travels in the m-th dimension; 2,mx represents the distance the second water particle travels in the m-th dimension; n,m y represents the distance the nth water particle travels in the mth dimension; 1,m This represents the first individual in m dimensions; y 2,m This represents the second individual in m dimensions; y n,m This represents the nth individual in m dimensions.

[0052] S3.2: Calculate the fitness value based on the variable values ​​and store it in array f. i In this context, the fitness values ​​are as follows:

[0053]

[0054] Where: f1([y 1,1 y 1,2 … y 1,m [) represents the fitness set of the first water particle across all dimensions.

[0055] [y 1,1 y 1,2 … y 1,m ] represents the set of individuals of the first water particle in each dimension.

[0056] f2([y 2,1 y 2,2 … y 2,m [) represents the fitness set of the second water particle across all dimensions.

[0057] [y 2,1 y 2,2 … y 2,m ] represents the set of individuals of the second water particle in each dimension.

[0058] f n ([y n,1 y n,2 … y n,m [) represents the fitness set of the nth water particle across all dimensions.

[0059] [y n,1 y n,2 … y n,m ] represents the set of individual n-th water particles in each dimension.

[0060] i = 1, 2, 3, ..., n represents the index of the i-th water particle in the population; j = 1, 2, 3, ..., m represents the dimension index of each particle, and the fitness value f of the i-th particle. i The definition is as follows:

[0061]

[0062] Where: f(y) i,1 ,y i,2 ,…,y i,m Let represent the fitness set of the i-th water particle in each dimension.

[0063] S3.3: In the process of water transport from soil to leaves in plants, the movement of water particles under the action of hydrostatic pressure gradient is one of the main driving mechanisms. This process usually occurs in structures such as soil pores, xylem vessels, and intercellular spaces. The transport rate J of water particles in the soil, membrane structure, and root system is often described by the following mathematical model:

[0064]

[0065] Where: t is the number of iterations. and Let represent the water flow velocity of the i-th water particle in the j-th dimension at the next iteration step and the current iteration step, respectively; σ represents the pressure gradient of the i-th particle in the j-th dimension; σ is the reflection coefficient, which takes values ​​in the range [0,1]. L represents the optimal pressure gradient of the i-th particle in the j-th dimension, obtained under the combined effects of randomness and systematicity, and is calculated by equation (16); p The hydraulic conductivity of the membrane is set to 1×10⁻⁶. -9 ms -1 Pa -1 r1 and r2 are random numbers uniformly distributed in the interval [0, 1]. λ is the inertial coefficient used to regulate the amplitude of water particle flow, calculated by equation (18); y g The gravitational term is represented and its calculation method is shown below:

[0066]

[0067] Where: ρ represents the density of water in the soil, root system, and plant membrane, and is taken as 1000 kg / m³. 3 g is the acceleration due to gravity, with a value of 9.81 m / s². 2 ; This represents the osmotic pressure gradient of the i-th particle in the j-th dimension, expressed by the formula... Calculations show that, where Π represents osmotic pressure. That is, the partial derivative of the i-th particle in the j-th dimension during the t-th iteration; Represents the hydrostatic pressure gradient. Let be the partial derivative of the i-th particle in the j-th dimension during the t-th iteration. Let be the flow distance of the i-th particle in the j-th dimension during the t-th iteration.

[0068] The optimal pressure gradient of the i-th particle in the j-th dimension obtained under the combined effects of randomness and systematicity. The calculation is as follows:

[0069]

[0070] in: Let α be the optimal pressure gradient of the i-th particle in the j-th dimension, and let α be a random vector defined on the interval [0, 1]. The calculation process is as follows:

[0071]

[0072] Where: rand represents a uniformly distributed random value taking values ​​in the interval [0,1]. The value indicates rounding. The inertia coefficient λ, used to regulate the amplitude of water particle flow, is calculated as follows:

[0073]

[0074] S3.4: The hydrostatic pressure P drives the horizontal flow of water particles in the soil, and its velocity J can be described by the following formula:

[0075]

[0076] in: For the partial derivative of hydrostatic pressure, The partial derivative represents the distance water particles travel, where K is the hydraulic conductivity, with a value of 1 × 10⁻⁶. -9 When the water content changes, the water flow velocity under the influence of the water content gradient θ can be expressed as:

[0077]

[0078] Where: Q is the volumetric flux density of water in the soil, dl is the average root length per unit volume of soil, and D is the permeability coefficient, which is 1×10⁻⁶. -9 m 2 / s.

[0079] In summary, considering that the cross-sectional area of ​​the root system in the soil is πa 2 And dθ=πa 2 dP, where a is the radius of the root, and π is taken as 3.14. Based on actual conditions, the formula for calculating the water flow velocity J can be derived as follows:

[0080]

[0081] Based on the above formula The calculation expression is as follows:

[0082]

[0083] Where: p is a constant with a value of 0.5, c1, c2, and c3 are all random variables that follow a uniform distribution on the interval [0,1], and parameter a is set to 1; δ represents the water flow velocity of the i-th individual in the j-th dimension at iteration number t; rand This represents the random decay factor.

[0084] δ rand The definition is as follows:

[0085]

[0086] Where: v represents the dynamic decay factor, and γ is a constant with a value of 1×10⁻⁶. -7 t and T represent the current iteration number and the maximum iteration number, respectively.

[0087] When water particles flow vertically, considering the work done against gravity, the total water potential φ is defined as follows:

[0088]

[0089] Where: ψ is the water potential, x is the vertical distance of the water particle in the gravitational field, and ρgx represents the gravitational potential energy component;

[0090] The vertical flow of water particles can be expressed by the following formula:

[0091]

[0092] According to dθ=πa 2 The relationship of dP can be derived from the following expression:

[0093]

[0094] Based on the above formula, under the action of vertical water flow, The calculation expression is as follows:

[0095]

[0096] Where: c4, c5, and c6 are all random variables that follow a uniform distribution on the interval [0,1].

[0097] S3.5: In the soil-to-root surface stage, the algorithm introduces pressure gradients and hydraulic conduction mechanisms to simulate the basic migration process of water particles driven by hydrostatic pressure, and incorporates individual historical optimal information for directional shifting to accelerate the convergence speed of candidate solutions in the initial stage. Considering the influence of plant roots on water absorption, the migration behavior of water in the soil is more complex. Against this background, the mathematical model describing the movement process of water particles under root absorption is shown below:

[0098]

[0099] Where: S(x,t) represents the root system's water absorption term. This function takes into account the root system's distribution characteristics at depth x and time t, and is used to characterize the process of plant roots absorbing water with time and soil depth.

[0100] Considering the vertical migration of water under the influence of gravity, the motion of water particles can be described by the following formula:

[0101]

[0102] Wherein, the water potential ψ can be calculated from the difference between the hydrostatic pressure P and the osmotic pressure Π, i.e., ψ = P - Π. Based on this, the water flow velocity J can be expressed as:

[0103]

[0104] Therefore, the hydrostatic pressure gradient can be calculated using the following formula:

[0105]

[0106] Where: r5, r6, and c7 are independent random variables that are uniformly distributed on the interval [0,1].

[0107] S3.6: The kinetic model of water particle migration from the root surface to the xylem in the root surface to xylem stage is expressed as follows:

[0108]

[0109] Where: L p χ represents radial water conductivity, and χ is the dimensionless reflection coefficient. p,s With ψ p,x ψ represents the hydrostatic pressure at the root surface and in the membrane, respectively. o,s With ψ o,x These represent the osmotic pressure at the root surface and within the membrane, respectively.

[0110] Based on this, the hydrostatic pressure gradient can be expressed by the following equation:

[0111] Where: χ takes the value of 0.5, and r7, r8, and c8 are independent random variables that are uniformly distributed on the interval [0,1].

[0112] S3.7: In the xylem-to-blade stage, the algorithm introduces viscous drag and fluid enhancement mechanisms to control the steady decrease of particle velocity and promote the gradual convergence of the solution towards the optimum. At the blade end, the final release of water is simulated, corresponding to the particle convergence process. Through velocity coupling terms and mutation control mechanisms, fine-grained adjustments are made to optimal solutions, thereby improving the algorithm's global optimization performance. The water pressure gradient from xylem to blade can be expressed by the following formula:

[0113]

[0114] Where: η represents the dynamic solution viscosity, with a value of 0.0018, a is the radius of the plant container, with a value of 1, and c9 represents a random variable that follows a uniform distribution in the interval [0,1].

[0115] S3.8: Crossover is a strategy to enhance the population diversity of optimization algorithms. In the IWUTP algorithm, crossover is introduced to improve the algorithm's global search capability and reduce the risk of getting trapped in local optima. Specifically, this operation first selects two individuals from the historical best individuals as parents, and performs crossover according to the crossover rate CR, the calculation formula of which is shown below:

[0116]

[0117] Where, r d Let z represent a random number that follows a uniform distribution on the interval [0,1]. d Let be the d-th dimension component of the offspring z, with a crossover rate CR of 0.7, and p i1 With p i2 These are two parent individuals selected from the historical best individuals. (p) i2 ) d This represents the second parent individual selected from the historical best individuals in dimension d, (p i1 ) d Let represent the first parent individual selected from the historical best individuals in d dimensions.

[0118] During the mutation phase, mutation operations are performed sequentially on each dimension of the individual, resulting in the mutated individual z. d The calculation formula is as follows:

[0119] z d =lb d +u d ·(ub d -lb d (h) d <MR) (38);

[0120] Where: u d with h dThis represents a random number that follows a uniform distribution on the interval [0,1]. MR is the rate of variation, set to 0.1. lb d with ub d These represent the lower and upper bounds of the d-th dimension variable, respectively. After the mutation operation is completed, the fitness values ​​of the mutated individuals are compared with those of the original individuals, and individuals with better fitness are retained to update the population.

[0121] S3.9: Elite selection is a strategy to improve the stability and convergence accuracy of optimization algorithms. Its core idea is to retain the best individual in the current population during each generation of evolution, thereby avoiding the loss of high-quality solutions due to search perturbations in subsequent iterations. In the IWUTP algorithm, the elite strategy is used to record the individual with the best fitness in the current water particle swarm, and after each iteration, this elite individual replaces the individual with the worst fitness in the population. The mathematical expression for elite individual extraction is as follows:

[0122]

[0123] Where: i * The index p represents the elite individual. i Let f(p) be the historical optimal position of the i-th particle, argmin represent the independent variable that minimizes the objective function, and N be the population size; i ) represents the fitness value of the i-th particle.

[0124] After determining the optimal particle index, the corresponding optimal particle position and its fitness are calculated using the following formula:

[0125]

[0126] Where: f elite f(x) represents the optimal particle fitness value. elite ) represents the fitness value of an elite individual, x elite p is the individual with the best fitness in the current population. i * indicates the position of an elite individual. Meanwhile, the index of the worst-fit individual in the population is determined by the following mathematical expression:

[0127]

[0128] Where: j * f(x) represents the index of the worst-fit individual in the population. i ) represents the fitness of the i-th individual in the population.

[0129] Let x i For the i-th individual in the current population, if This represents the fitness value of the worst-fitting individual in the population. Then, an elite individual replacement is performed, using the following formula:

[0130]

[0131] in: This indicates the position of the individual with the worst fitness in the population.

[0132] Furthermore, the current best individual is updated to the historical best solution using the following formula.

[0133]

[0134] Based on equations (8) to (42), the IWUTP algorithm is finally obtained;

[0135] S3.10: Variational Mode Decomposition (VMD) utilizes variational model theory to adaptively decompose a signal into several eigenmode functions with specific spectral characteristics through optimization methods. VMD solves the following variational constraint problem through optimization, and its formula is shown below:

[0136]

[0137] In the above formula: K represents the number of modes to be decomposed, and k is the number of mode decompositions; u k (t) is the k-th mode; ω k δ(t) is the center frequency of the k-th mode; δ(t) is the Dirac function; It is a multiplication operator; u k Let f(t) represent the time-domain representation of the k-th mode, and let f(t) represent the original signal.

[0138] S3.11: To optimize the solution of unconstrained variational problems, we introduce the Lagrange multiplier λ and the penalty factor α. The specific formulas are shown below:

[0139]

[0140] In equation (45): L({u k},{ω k},λ) denote the extended Lagrangian function, {u k} represents K modal functions, {ω k Let} denote the center frequencies of the K modes, and λ(t) be the Lagrange multiplier that varies with time.

[0141] Among them, the modal component u k (t) and center frequency ω k The update formula is as follows:

[0142]

[0143] In the above formula, Let represent the frequency domain representation of the k-th modal component in the (n+1)-th iteration. This represents the center frequency of the k-th modal component in the (n+1)-th iteration. This is represented as the frequency domain representation of the i-th mode. Let f(t) and u represent respectively. k The frequency domain expressions for u(t) and u(t); ω represents the frequency; n represents the number of iterations.

[0144] S3.12: The number of modes K and the penalty factor α of variational mode decomposition (VMD) are used as the two dimensions of the IWUTP algorithm, and the permutation entropy is used as the fitness function. Specifically, during the optimization process, each IWUTP individual corresponds to a set of (K, α) parameters. The population of the IWUTP algorithm is initialized by formulas (8) to (11). The algorithm will perform VMD decomposition based on the current parameters, calculate the permutation entropy of each mode, and use the minimum or average value of these entropy values ​​as the fitness evaluation index, thereby realizing the adaptive adjustment of parameters.

[0145] The formula for permutation entropy is expressed as follows:

[0146]

[0147] In equation (48): nu represents the embedding dimension, l i Let represent the probability of the i-th permutation pattern, obtained statistically, and nu! represent the total number of all possible permutations.

[0148] In step 4, GRU is used to predict each intrinsic mode function, and the prediction results of each mode are superimposed point by point at each time step. Let the input sequence be X = (x1, x2, ..., x...). t ), x1, x2, ..., x t Input data at time steps 1, 2, ..., t.

[0149] The hidden state is H = (h1, h2, ..., h t ), h1, h2, ..., h t Let represent the hidden state vectors at time steps 1, 2, ..., t, respectively.

[0150] The calculation and update process of GRU is defined by the following formula:

[0151] r t =σ(W r ·[h t-1 ,x t ]+b r (49);

[0152] z t =σ(W z·[h t-1 ,x t ]+b z (50);

[0153]

[0154] Among them, z t and r t These represent updating the door and resetting the door, respectively; h t and represents the hidden state and the candidate hidden state, respectively; W and b are the weight matrix and bias vector of the corresponding terms, respectively; the symbol · represents dot multiplication; t and t-1 represent the current time step and the previous time step, respectively; tanh is the hyperbolic tangent activation function; the symbol * represents element-wise multiplication. t-1 W represents the hidden state at the previous moment. r W represents the weight matrix input to the reset gate. z Let b represent the weight matrix input to the update gate. r b represents the bias vector input to the reset gate. z This represents the bias vector input to the update gate, [h] t-1 ,x t [h] indicates concatenating two vectors. t-1 *r t ,x t ] indicates that h t-1 With r t Multiply element by element, then multiply with vector x t The concatenation is performed, where σ(·) represents the Sigmoid activation function. In step 4, the intrinsic mode functions obtained in step 3 are first input into the GRU model for training. A sliding window is used to divide the dataset, specifically selecting the last 100 time steps as the test set and the rest as the training set, thus obtaining the input training set X. train Input test set X test and the actual training set sequence Y train During training, the Adam optimization algorithm is used to iteratively update the trainable parameters in the GRU model, including weights and biases, with the goal of minimizing the loss function, thereby improving the model's predictive performance. All hyperparameter optimizations are based on experimental tuning to ensure optimal training results.

[0155] After training, the prediction training set sequence of the first-layer model is obtained. The test set from step 3 is then input into the trained GRU model to obtain the prediction test set sequence of the first-layer model. This allows for the construction of corresponding prediction training set sequences and prediction test set sequences.

[0156] Simultaneously retain the actual training set sequence Y of PDI.train This will serve as a reference for subsequent modeling.

[0157] In step 5, for an ELM model with L hidden nodes, let the training sample set contain N samples, denoted as... Where: x i Let y be the input feature of the i-th sample. i Indicates the corresponding target output; (x i ,y i ) represents the data pair of the i-th training sample. Represents a g-dimensional real space. It represents a c-dimensional real space.

[0158] The model aims to efficiently solve for the optimal output weight β by understanding the relationship between the hidden layer output and the output weights. Therefore, the weight matrix between the hidden layer and the output layer can be expressed as:

[0159]

[0160] In equation (53): H represents the output matrix after the sample is mapped through the hidden layer, H T Let H be the transpose of H, T be the output target matrix, I be the unit diagonal matrix, and C be the regularization parameter.

[0161] After introducing the kernel function, the resulting KELM model can map the input samples from the original n-dimensional space to a higher-dimensional feature space. The kernel matrix expression is defined as follows:

[0162]

[0163] In equation (54): Ω(i,j) represents the kernel matrix element, h(x i h(x) j ) are samples x i x j Mapping value in the hidden layer; K(x) i ,x j ) represents sample x i x j The kernel function.

[0164] Common kernel functions used in KELM models include linear kernels, polynomial kernels, and Gaussian kernels. To effectively handle nonlinear data characteristics and enhance the model's ability to learn complex patterns, this invention selects the Gaussian kernel function as the kernel function for the KELM model, and its expression is defined as follows:

[0165]

[0166] In equation (55): σ is the kernel parameter of the Gaussian kernel function.

[0167] Based on the computational mechanism of the KELM model, the output of the model can be expressed as:

[0168]

[0169] In equation (56): f(x) represents the model's predicted output for the input sample x, h(x) represents the feature matrix of the input sample x after mapping, and K(x,x1) represents the kernel function matrix of samples x and x1. N ) represents the sample x and x N The kernel function matrix.

[0170] The KELM model is optimized using the IWUTP algorithm obtained from embedded equations (8) to (42) and used as the second-layer prediction model. Specifically, the root mean square error is used as the fitness function of the IWUTP algorithm to evaluate the regularization coefficient and kernel parameters in the KELM model, and the two hyperparameter values ​​corresponding to the minimum fitness function are selected as the final optimization results. The prediction training set, prediction test set, and actual PDI training set sequence output by GRU in step 4 are used as input to finally obtain the actual PDI prediction results.

[0171] It also includes step 6: using root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAPE), and coefficient of determination (R²). 2 The four evaluation indicators are calculated using the following formulas:

[0172]

[0173] In the above formula: n represents the number of samples, y i and Let these represent the actual value and the predicted value of the i-th sample, respectively. It is the mean of the actual values, that is

[0174] A deterioration trend prediction system for pumped storage units based on multinomial XGBoost and a two-layer prediction model is proposed. The system includes: a data acquisition module, a data generation module, a model training module, and a model prediction module.

[0175] The data acquisition module acquires the original vibration data of the pumped storage unit, as well as key operating characteristics such as power generation head, active power, and guide vane opening. The historical data consists of the original vibration data of the pumped storage unit, as well as power generation head, active power, and guide vane opening, acquired at 30-minute intervals. This data will serve as the key input for the subsequent prediction model.

[0176] The data generation module constructs a health status prediction model based on the multinomial XGBoost algorithm and introduces the SHAP method to evaluate the contribution of each input feature to the model output, thereby realizing the interpretability analysis of the model.

[0177] Vibration data under deteriorated operating conditions to be predicted are input into a trained health status prediction model to obtain corresponding predicted vibration data. A Predicted Displacement Indicator (PDI) is constructed based on the deviation between the predicted and actual observed values ​​to quantify the degree of equipment degradation. The constructed PDI sequence is decomposed using the IWUTP-optimized VMD method to obtain K intrinsic mode functions (IMFs) for acquiring multi-scale temporal information of degradation characteristics.

[0178] The model training module uses GRU to predict each intrinsic mode function (EMF) separately, and then superimposes the prediction results of each mode at each time step to construct the corresponding prediction training set sequence and prediction test set sequence. Simultaneously, the actual training set sequence of PDI is retained. The prediction training set, prediction test set output by GRU, and the actual training set of PDI are used as input to construct a KELM model based on IWUTP optimization, further modeling and predicting the degradation trend, and finally outputting a more accurate degradation trend prediction result. This trains a vibration trend prediction model for pumped storage units based on GRU-IWUTP-KELM.

[0179] The model prediction module takes the PDI data of the pumped storage unit to be predicted and, after going through the same data preprocessing process, imports it as input data into the trained pumped storage unit vibration trend prediction model based on GRU-IWUTP-KELM to predict the deterioration trend of the pumped storage unit in the future time period.

[0180] This invention discloses a method for predicting the degradation trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model. The technical effects are as follows:

[0181] 1) This invention utilizes the multinomial XGBoost algorithm to construct a health status prediction model and introduces the SHAP method to evaluate the contribution of each input feature to the model output. The vibration data of the deteriorating working conditions to be predicted is input into the trained health status prediction model to obtain the corresponding predicted vibration data. A PDI (Predicted Displacement Index) is constructed based on the deviation between the predicted and actual observed values ​​to quantify the degree of equipment deterioration.

[0182] 2) This invention employs the IWUTP-VMD method to decompose the constructed PDI sequence into several sub-sequence components, reducing nonlinearity and nonstationarity. Simultaneously, a GRU is used to predict each intrinsic mode function, and the prediction results for each mode are superimposed point-by-point at each time step to construct corresponding prediction training set and prediction test set sequences. The actual PDI training set sequence is retained. The prediction training set, prediction test set, and actual PDI training set output by the GRU are used as input to construct an IWUTP-optimized KELM model, further modeling and predicting the degradation trend, ultimately outputting a more accurate degradation trend prediction result.

[0183] 3) By comparing and analyzing experimental results, this invention verifies the superiority of the method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model. It has achieved good prediction results on historical data, and the proposed method significantly improves the stability of the prediction model and the accuracy of the prediction results. Attached Figure Description

[0184] The present invention will be further described below with reference to the accompanying drawings and embodiments:

[0185] Figure 1 This is a flowchart of the present invention.

[0186] Figure 2 This is a summary diagram of the importance of SHAP features in this invention.

[0187] Figure 3 This is a diagram showing the decomposition results of the IWUTP-VMD method in this invention.

[0188] Figure 4 This is a comparison chart of the predicted and actual values ​​using the model proposed in this invention.

[0189] Figure 5 This is a comparison chart of error box curves using the model proposed in this invention. Detailed Implementation

[0190] A method for predicting the degradation trend of pumped storage units based on polynomial XGBoost and a two-stage prediction model is proposed. First, a health status assessment model is constructed using polynomial XGBoost. By incorporating features from polynomial transformation, the model's nonlinear fitting ability is significantly enhanced. The SHAP method is then used to explain the contribution of each feature to the model output, improving the model's interpretability. Based on the prediction results of the health status model, vibration data is extracted to construct the degradation trend indicator (PDI). Subsequently, an improved IWUTP algorithm is introduced to optimize the VMD parameters and efficiently decompose the PDI signal. Finally, a two-stage prediction model is constructed. First, a gated cyclic unit (GRU) is used to capture the temporal characteristics of the PDI. Then, IWUTP-KELM is used to refine the prediction of the GRU output, thereby obtaining an accurate prediction of the degradation trend. The prediction process is as follows: Figure 1 As shown, the specific steps include the following:

[0191] Step a: Obtain the raw vibration data of the pumped storage unit, as well as key operating characteristics such as power generation head, active power, and guide vane opening, as shown in Table 1.

[0192] Table 1. Data Collection Objects for Pumped Storage Unit Operations

[0193] Collection objects unit Hydropower head m Active power MW Guide vane opening ％ vibration amplitude μm

[0194] A health status prediction model is constructed based on the multinomial XGBoost algorithm, and the SHAP method is introduced to evaluate the contribution of each input feature to the model output, thereby realizing the interpretability analysis of the model:

[0195] Step a1: Obtain the raw vibration data of the pumped storage unit, as well as key operating characteristics such as power generation head, active power, and guide vane opening.

[0196] Step a2: Construct a health status prediction model based on the multinomial XGBoost algorithm. The calculation formula for XGBoost is shown below:

[0197] XGBoost is a high-efficiency machine learning algorithm based on the gradient boosting framework. It improves prediction performance by progressively building an additive model and integrating multiple weak classification or regression decision trees. It optimizes model parameters by minimizing an objective function that includes a loss function and a regularization term, and supports parallel computation and various regularization strategies, effectively preventing overfitting and improving the model's generalization ability. Let the training samples be T = (x1, y1), (x2, y2), ..., (x... n ,y n ), i = 1, 2, 3, ..., n where x i Let y be the feature vector of the i-th sample. i For actual values, assume each decision tree is modeled. in For the predicted values, the overall objective function is as follows:

[0198]

[0199] Specifically, the squared error is used as the loss function, and its expression is as follows:

[0200]

[0201] In equation (2), L represents the loss function, which measures the error between the true value and the predicted value, and n is the number of samples.

[0202] A regularization term is introduced to control model complexity; its expression is shown below:

[0203]

[0204] Where T is the number of regression trees, Ω is the model complexity penalty term, γ is the regularization coefficient used to control the number of trees, λ represents the L2 regularization coefficient of the leaf node weights, and π is the leaf node weight vector of the regression trees, which contains the feature weights composed of the predicted values ​​output by all leaf nodes in the model.

[0205] The Boosting algorithm employs an additive model approach, progressively stacking several weak regressors to construct a strong regression model. In the t-th iteration, the model's prediction is the sum of the predictions from the previous t-1 iterations and the prediction from the current t-th tree. Essentially, Boosting optimizes the predictive performance of the current weak classifier in the t-th iteration by fitting the residuals of the current model from the previous t-1 iterations, thereby minimizing the overall objective function. In constructing the healthy model, hydroelectric head, active power, and guide vane opening are selected as input features, denoted as J = [j1, j2, j3]. j1, j2, and j3 represent hydroelectric head, active power, and guide vane opening, respectively.

[0206] A polynomial extension method is introduced based on the original XGBoost features. Specifically, based on the original XGBoost features, a fourth-order polynomial extension is used to generate power terms and multidimensional interaction terms (such as square terms, cube terms, and products between features) to describe nonlinear degradation relationships and the interactive effects of multiple features; the specific formula is as follows:

[0207]

[0208] Where J represents the original feature vector, a s Indicates the corresponding feature j s Let φ(J) be the polynomial power, φ(J) be the feature set after polynomial expansion, and s be the feature index, with values ​​ranging from 1 to 3. It is a set of non-negative integers.

[0209] Step a3: The SHAP method is introduced to evaluate the contribution of each input feature to the model output, thereby achieving the interpretability analysis of the model. SHAP is described as follows:

[0210] To enhance the interpretability of health models, this invention introduces the SHAP method based on Shapley value theory in game theory. SHAP calculates the marginal contribution of each feature to the model output, thereby making the health model more interpretable. The SHAP method quantifies the contribution of input features to the health state prediction model output, focusing on identifying the impact of key operating parameters such as hydroelectric head, active power, and guide vane opening on the prediction results, thus improving the model's interpretability. Based on the model's structural characteristics, the Tree Explainer method is used to calculate and interpret the SHAP values ​​of XGBoost and Polynomial XGBoost models. By constructing a summary graph of SHAP feature importance, the average contribution and fluctuation range of key parameters such as hydroelectric head, active power, and guide vane opening to the model output are visually displayed, revealing the core role of these features in health state modeling.

[0211] In a single decision tree, the SHAP value is calculated by analyzing the splitting of a feature across all decision paths. For a given feature, its SHAP value represents the expected marginal contribution of that feature to the predicted output across all possible paths, expressed as follows:

[0212]

[0213] Where, the Shapley value is the feature x i Contribution to the model output Represents all cases containing x i The summation is performed on path p, where Pr(p) is the sampling probability or weight of path p, and Δf(p,x) is the summation. i ) represents the feature x in path p. i This results in an increase in the model output.

[0214] Finally, the SHAP values ​​are summed to quantify the contribution of each feature to the overall prediction result of the model.

[0215]

[0216] Step b: Input the vibration data of the deteriorating working conditions to be predicted into the trained health status prediction model to obtain the corresponding predicted vibration data. Construct a Predicted Displacement Indicator (PDI) based on the deviation between the predicted and actual observed values ​​to quantify the degree of equipment deterioration. The calculation formula is as follows:

[0217]

[0218] Where F t V represents the vibration value predicted by the health model at time t. t These are the actual measured values ​​at the corresponding times.

[0219] Step c: Improve the WUTP algorithm by introducing elite selection and crossover mutation to obtain the IWUTP algorithm, the formula of which is shown below:

[0220] Step c1: The WUTP algorithm models candidate solutions in the optimization problem as pressure gradient vectors. By simulating the water transport process from soil to leaves in plants, it uses water flow velocity to drive dynamic adjustment and iterative updates of candidate solutions in multiple stages. This process can be divided into three consecutive biomimetic stages: soil to root surface, root surface to xylem, and xylem to leaves, ultimately completing an efficient search from the initial solution to the global optimum. The specific initialization steps of the algorithm are as follows:

[0221] y i =[lb j +rand×(ub j -lb j )]

[0222] Where lb j with ub j represents the lower and upper bounds of the j-th dimension, respectively, and rand represents a random variable that follows a uniform distribution on the interval [0,1].

[0223] After initializing the population, the plant leaves and moving water particles were initialized separately. The initialization steps for the plant leaves are as follows:

[0224] Leaves = [l1, l2, ..., l i ,…,l m ]

[0225] Where m represents the problem dimension, l i Let i be the i-th variable. The water particles in the plant are initialized as follows:

[0226] Particles = [p1, p2, p3, ..., p n ]

[0227] Where n is the population size, and the water particles need sufficient hydrostatic pressure P to drive them to a position close to the blade. At this point, the water particles are under a pressure gradient... Under the influence of [the flow], the roots travel a distance x from the soil to the leaves. Specifically, the plant roots travel a distance x from the soil... i =0 absorbs water particles, and these particles travel a distance x within the plant leaf. f The flow distance can be determined by The results were obtained through calculation. In summary, the population can be represented as follows:

[0228]

[0229] Step c2: Calculate the fitness value based on the variable values ​​and store it in array f. i In this context, the fitness values ​​are as follows:

[0230]

[0231] Where i = 1, 2, 3, ..., n represents the index of the i-th water particle in the population, j = 1, 2, 3, ..., m represents the dimension index of each particle, and the fitness value of the i-th particle is defined as follows:

[0232] f i =f(y i,1 ,y i,2 ,…,y i,m )

[0233] Step c3: In the process of water transport from soil to leaves, the movement of water particles under the action of hydrostatic pressure gradient is one of the main driving mechanisms. This process usually occurs in structures such as soil pores, xylem vessels, and cell wall interstitials. The transport rate J of water particles in the soil, membrane structures, and root system is often described using the following mathematical model:

[0234]

[0235] Where t is the number of iterations. and Let represent the water flow velocity of the i-th water particle in the j-th dimension at the next iteration step and the current iteration step, respectively. Let σ represent the pressure gradient of the i-th particle in the j-th dimension. σ is the reflection coefficient, which takes values ​​in the range [0,1]. Let L represent the optimal pressure gradient of the i-th particle in the j-th dimension, obtained under the combined effects of randomness and systematicity. p The hydraulic conductivity of the membrane is set to 1×10⁻⁶. -9 ms -1 Pa -1 r1 and r2 are random numbers uniformly distributed in the interval [0,1]. λ is the inertia coefficient used to regulate the amplitude of water particle flow, calculated by equation (16). g The gravitational term is represented and its calculation method is shown below:

[0236]

[0237] Where ρ represents the density of water in the soil, roots, and plant membrane, and is taken as 1000 kg / m³. 3 g is the acceleration due to gravity, with a value of 9.81 m / s².2 , This represents the osmotic pressure gradient of the i-th particle in the j-th dimension, expressed by the formula... Calculations show that Π represents osmotic pressure. This represents the hydrostatic pressure gradient.

[0238]

[0239] in Let α be the optimal pressure gradient of the i-th particle in the j-th dimension, and let α be a random vector defined on the interval [0,1]. The calculation process is as follows:

[0240]

[0241] Where rand represents a uniformly distributed random value that takes values ​​in the interval [0,1]. This indicates rounding.

[0242]

[0243] Step c4: The hydrostatic pressure P drives the horizontal flow of water particles in the soil, and its velocity J can be described by the following formula:

[0244]

[0245] in: For hydrostatic pressure gradient, This represents the distance water particles travel, where K is the hydraulic conductivity, with a value of 1 × 10⁻⁶. -9 When the water content changes, the water flow velocity under the influence of the water content gradient θ can be expressed as:

[0246]

[0247] Where Q is the volumetric flux density of water in the soil, dl is the average root length per unit volume of soil, and D is the permeability coefficient, which is 1 × 10⁻⁶. -9 m 2 / s.

[0248] In summary, considering that the cross-sectional area of ​​the root system in the soil is πa 2 And dθ=πa 2 dP, where a is the radius of the root, and π is taken as 3.14. Based on the actual situation, the formula for calculating the water flow velocity J can be derived as follows:

[0249]

[0250] Based on the above formula The calculation expression is as follows:

[0251]

[0252] Where: p is a constant with a value of 0.5, c1, c2, and c3 are all random variables uniformly distributed on the interval [0,1], parameter a is set to 1, and δ rand The definition is as follows:

[0253]

[0254] Where v represents the dynamic decay factor, and γ is a constant with a value of 1 × 10⁻⁶. -7 t and T represent the current iteration number and the maximum iteration number, respectively.

[0255] When water particles flow vertically, considering the work done against gravity, the total water potential φ is defined as follows:

[0256] φ=ψ+ρgh=P-Π+ρgx

[0257] Where ψ represents the water potential, x is the vertical distance of the water particle in the gravitational field, and ρgx represents the gravitational potential energy component. The vertical flow of water particles can be expressed by the following formula:

[0258]

[0259] According to dθ=πa 2 The relationship of dP can be derived from the following expression:

[0260]

[0261] Based on the above formula, under the action of vertical water flow, The calculation expression is as follows:

[0262]

[0263] Where c4, c5, and c6 are all random variables that follow a uniform distribution on the interval [0,1].

[0264] Step c5: In the soil-to-root surface stage, the algorithm introduces pressure gradients and hydraulic conduction mechanisms to simulate the basic migration process of water particles driven by hydrostatic pressure, and combines individual historical optimal information for directional shifting to accelerate the convergence speed of candidate solutions in the initial stage. Considering the influence of plant roots on water absorption, the migration behavior of water in the soil is more complex. Against this background, the mathematical model describing the movement process of water particles under root absorption is shown below:

[0265]

[0266] Where S(x,t) represents the root system's water absorption term. This function takes into account the root system's distribution characteristics at depth x and time t, and is used to characterize the process of plant roots absorbing water with time and soil depth.

[0267] Considering the vertical migration of water under the influence of gravity, the motion of water particles can be described by the following formula:

[0268]

[0269] The water potential ψ can be calculated from the difference between the hydrostatic pressure P and the osmotic pressure Π, i.e., ψ = P - Π. Based on this, the water flow velocity J can be expressed as:

[0270]

[0271] Therefore, the hydrostatic pressure gradient can be calculated using the following formula:

[0272]

[0273] Where r5, r6, and c7 are independent random variables that are uniformly distributed on the interval [0,1].

[0274] Step c6: In the root surface to xylem stage, the algorithm introduces random perturbations and jump mechanisms to enhance search diversity and reduce the risk of getting trapped in local optima. The dynamic model of water particle migration from the root surface to the xylem is expressed as follows:

[0275] Q = L p (ψ p,s -ψ p,x +χ(ψ o,s -ψ o,x ))

[0276] Where L p χ represents radial water conductivity, and χ is the dimensionless reflection coefficient. p,s With ψ p,x ψ represents the hydrostatic pressure at the root surface and in the membrane, respectively. o,s With ψ o,x Let represent the osmotic pressure at the root surface and within the membrane, respectively. Based on this, the hydrostatic pressure gradient can be expressed by the following equation:

[0277]

[0278] Where χ takes the value of 0.5, and r7, r8, and c8 are independent random variables that are uniformly distributed on the interval [0,1].

[0279] Step c7: In the xylem-to-blade stage, the algorithm introduces viscous drag and fluid enhancement mechanisms to control the steady decrease of particle velocity and promote the gradual convergence of the solution towards the optimal value. At the blade end, the final release of water is simulated, corresponding to the particle convergence process. Through velocity coupling terms and mutation control mechanisms, fine-grained adjustments are made to optimal solutions, thereby improving the algorithm's global optimization performance. The water pressure gradient from xylem to blade can be expressed by the following formula:

[0280]

[0281] Where η represents the dynamic solution viscosity, with a value of 0.0018, a is the radius of the plant container, with a value of 1, and c9 represents a random variable that follows a uniform distribution in the interval [0,1].

[0282] Step c8: Crossover is a strategy to enhance the population diversity of optimization algorithms. In the IWUTP algorithm, crossover is introduced to improve the algorithm's global search capability and reduce the risk of getting trapped in local optima. Specifically, this operation first selects two individuals from the historical best individuals as parents, and performs crossover according to the crossover rate CR, the calculation formula of which is shown below:

[0283]

[0284] Where, r d Let z represent a random number that follows a uniform distribution on the interval [0,1]. d Let be the d-th dimension component of the offspring z, with a crossover rate CR of 0.7, and p i1 With p i2 These are two parent individuals selected from the best individuals in history.

[0285] During the mutation phase, mutation operations are performed sequentially on each dimension of the individual, and the calculation formula is as follows:

[0286] z d =lb d +u d ·(ub d -lb d (h) d <MR)

[0287] Where: u d with h d This represents a random number that follows a uniform distribution on the interval [0,1]. MR is the rate of variation, set to 0.1. lb d with ub d These represent the lower and upper bounds of the d-th dimension variable, respectively. After the mutation operation is completed, the fitness values ​​of the mutated individuals are compared with those of the original individuals, and individuals with better fitness are retained to update the population.

[0288] Step c9: Elite selection is a strategy to improve the stability and convergence accuracy of optimization algorithms. Its core idea is to retain the best individual in the current population during each generation of evolution, thereby avoiding the loss of high-quality solutions due to search perturbations in subsequent iterations. In the IWUTP algorithm, the elite strategy is used to record the individual with the best fitness in the current water particle swarm, and after each iteration, this elite individual replaces the individual with the worst fitness in the population. The mathematical expression for elite individual extraction is as follows:

[0289]

[0290] Where p i Let be the historical best position of the i-th particle, argmin represent the independent variable that minimizes the objective function, and N is the population size.

[0291] After determining the optimal particle index, the corresponding optimal particle position and its fitness are calculated using the following formula:

[0292]

[0293] Where f represents the fitness function. elite ) represents the fitness value of an elite individual, x elite This represents the individual with the best fitness in the current population. Simultaneously, the index of the individual with the worst fitness in the population is determined using the following mathematical expression:

[0294]

[0295]

[0296]

[0297] Furthermore, the current best individual is updated to the historical best solution using the following formula:

[0298]

[0299] Step c10: VMD utilizes variational model theory to adaptively decompose the signal into several eigenmode functions with specific spectral characteristics through optimization methods. VMD solves the following variational constraint problem through optimization, and its formula is shown below:

[0300]

[0301] In the formula: k is the number of mode decompositions, u k (t) is the k-th mode, ω k δ(t) is the center frequency of the k-th mode, and δ(t) is the Dirac function. It is a multiplication operator.

[0302] Step c11: To optimize the solution of the unconstrained variational problem, the Lagrange multiplier λ and the penalty factor α are introduced. The specific formulas are shown below:

[0303]

[0304] Among them, the modal component u k (t) and center frequency ω k The update formula is as follows:

[0305]

[0306] In the formula, Let f(t) and u represent respectively. k The frequency domain expressions for (t) and u(t), where ω represents the frequency and n represents the number of iterations.

[0307] Step c12: The number of modes and the penalty factor of VMD are used as two dimensions of IWUTP, and the permutation entropy is used as the fitness function to achieve adaptive adjustment of parameters.

[0308] Step d: GRU is a recurrent neural network that incorporates a gating mechanism to address the vanishing and exploding gradient problems encountered by traditional recurrent neural networks when processing long sequences. GRU regulates information flow through update and reset gates, effectively capturing long-term dependencies in time series data. The update gate controls the degree of fusion between the current state and the previous state, thus determining the extent to which historical information is retained in the current state. The reset gate adjusts the degree to which information from the previous time step is forgotten when calculating the current state. This structure allows GRU to selectively transfer and store information with a relatively small number of parameters, significantly improving training efficiency and model convergence speed. Let the input sequence be X = (x1, x2, ..., x...). t The hidden state is H = (h1, h2, ..., h t The calculation and update process of GRU is defined by the following formula:

[0309] r t =σ(W r ·[h t-1 ,x t ]+b r )

[0310]

[0311]

[0312] Among them, z t and r tThese represent the update gate and the reset gate, respectively. t and Let represent the hidden state and the candidate hidden state, respectively. W and b are the weight matrix and bias vector of the corresponding terms, respectively. The symbol · represents dot multiplication. t and t-1 represent the state at the current time step and the state at the previous time step, respectively. tanh is the hyperbolic tangent activation function. The symbol * represents element-wise multiplication.

[0313] Step e: Extreme Learning Machine (ELM) is a highly efficient single-hidden-layer feedforward neural network that shortens training time by randomly generating the weights and biases of the hidden layers. However, ELM has limited performance when dealing with nonlinear problems. To address this, KELM introduces a kernel function into ELM, replacing the random mapping in ELM with a stable kernel mapping, thereby enhancing the model's ability to handle nonlinear data and improving prediction accuracy and generalization performance. For an ELM model with L hidden nodes, let the training sample set contain N samples, denoted as... Where x i Let y be the input feature of the i-th sample. i This represents the corresponding target output. The model objective is to effectively solve for the optimal output weight β by understanding the relationship between the hidden layer output and the output weights. Therefore, the weight matrix between the hidden layer and the output layer can be expressed as:

[0314]

[0315] Where H represents the output matrix after the sample is mapped through the hidden layer, H T Let H be the transpose of H, T be the target output matrix, I be the unit diagonal matrix, and C be the regularization parameter.

[0316] By introducing a kernel function, KELM can map input samples from the original n-dimensional space to a higher-dimensional feature space, where the kernel matrix expression is defined as follows:

[0317] Ω=HH T

[0318] Ω(i,j)=h(x i )h(x j )=K(x i ,x j )

[0319] In the formula, Ω(i,j) represents the kernel matrix element, and h(x i h(x) j ) are samples x i x j Mapping values ​​in the hidden layer.

[0320] Common kernel functions in KELM include linear kernels, polynomial kernels, and Gaussian kernels. To effectively handle nonlinear data characteristics and enhance the model's ability to learn complex patterns, this invention selects the Gaussian kernel function as the kernel function for KELM, and its expression is defined as follows:

[0321]

[0322] In the formula, σ is the kernel parameter of the Gaussian kernel function.

[0323] Based on the computational mechanism of the KELM model, the output of the model can be expressed as:

[0324]

[0325] Finally, the regularization coefficient and kernel parameters of KELM are used as the two dimensions of IWUTP, and RMSE is used as the fitness function to achieve adaptive adjustment of the parameters.

[0326] Step f: To verify the effectiveness of the proposed method, this invention selects operating parameters, lower frame Y-axis vibration data, and lower guide vane Y-axis swing data collected by the monitoring system of Unit 2 of a domestic pumped storage power station for experiments. The sampling time interval is 30 minutes. A total of 4000 data vectors F(H(t), P(t), G(t)) of lower frame Y-axis vibration data and operating parameters from May 31, 2008 to December 25, 2009 are selected to construct a health status prediction model, where H(t) is the generating head, P(t) is the active power, and G(t) is the guide vane opening. These three are used as feature inputs for model training and evaluation. A total of 2854 samples of lower frame Y-axis vibration data and lower guide vane Y-axis swing data from January 2, 2010 to December 15, 2011 are selected to construct the PDI and conduct a degradation trend prediction experiment.

[0327] To construct a health status prediction model, the collected data was divided into training and testing sets. The first 3000 samples (corresponding to the period from May 31, 2008 to November 24, 2009) were used for model training, and the last 1000 samples (corresponding to November 24, 2009 to December 25, 2009) were used for model testing. To comprehensively evaluate the performance of different modeling methods in health status representation, this invention designed three modeling schemes: Backpropagation Neural Network (BPNN), XGBoost regression model, and Multinomial Extended XGBoost model. The BPNN uses a single hidden layer structure with 100 neurons, the ReLU activation function, and the Adam optimizer. The maximum number of iterations is set to 100, and the random seed is set to 42 to ensure the reproducibility of the experiment. The XGBoost model is set to 100 trees with a maximum depth of 2, a learning rate of 0.1, and the objective function is squared error regression. The random seed is also set to 42. In addition, to enhance the model's ability to express nonlinear features, the polynomial-extended XGBoost model introduces a fourth-order polynomial extension on the basis of the original input features.

[0328] In this invention, the PDI is constructed based on the relative deviation between the model's predicted values ​​and the actual observed values. Regarding data partitioning, 2000 sampling points from January 2, 2010 to February 6, 2011 were selected as the training set, and 854 sampling points from February 6, 2011 to December 15, 2011 were selected as the test set.

[0329] In the PDI prediction experiment, a sliding window was used to segment the data, with the number of windows set to 10. During the GRU model training phase, the first 744 sampling points were selected as the training set, and the subsequent 100 sampling points were used as the test set. In the KELM model, the training set contained 510 sampling points, the test set was also the last 100 sampling points, and the remaining data was used as the validation set. The learning rate for the GRU model was uniformly set to 1×10⁻⁶. -3 To prevent overfitting of the GRU layer in the combined model from affecting KELM training, this invention employs three-fold cross-validation for training the GRU layer in the combined model, and also applies the same three-fold validation strategy to the GRU model in the control group.

[0330] The root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAE), and coefficient of determination (R²) were used. 2 Four evaluation metrics are used to assess and analyze the model's prediction results. The specific formulas for the evaluation metrics are as follows:

[0331]

[0332] In the formula: n represents the number of samples, y i and Let these represent the actual value and the predicted value of the i-th sample, respectively. It is the mean of the actual values, that is

[0333] The decomposition results of the proposed IWUTP-VMD method are as follows: Figure 3 As shown. To fully demonstrate the effectiveness of the GRU-IWUTP-KELM pumped storage unit vibration trend prediction method, the proposed model GRU-IWUTP-KELM was compared with five other models: WUTP-KELM, IWUTP-KELM, MOPSO-KELM, and GRU. The comparison results of the proposed model are shown below. Figure 4 As shown, RMSE, MAE, MAPE, and R are used. 2 The predictive performance of all experimental models was quantitatively evaluated. The evaluation index results for different models are shown in Table 2, and the error box plots are shown below. Figure 5 As shown.

[0334] from Figure 4 As can be seen, the pumped storage PDI prediction curve of the proposed model GRU-IWUTP-KELM largely matches the actual curve, demonstrating high prediction accuracy and confirming the superiority and effectiveness of the proposed model.

[0335] Table 2 Evaluation metrics for different models

[0336]

[0337]

[0338] Table 2 shows that, in the Y direction of the lower rack, the WUTP-KELM model optimized by WUTP achieves good prediction results, with an RMSE of 1.604, a MAE of 1.436, a MAPE of 21.07%, and an R² of 1.604. 2 The R² value was 0.9891. After introducing the improved IWUTP algorithm, the RMSE and MAE of the IWUTP-KELM model decreased to 1.587 and 1.419, respectively, and the MAPE decreased to 20.90%, showing that the algorithm improved the model's ability to escape local optima, but the improvement on KELM performance was limited. Further optimization of the MOPSO-KELM model reduced the RMSE and MAE to 1.409 and 1.257, respectively, and the MAPE to 18.62%, R²... 2The accuracy was improved to 0.9913. After introducing the GRU model, the prediction performance was significantly improved, with RMSE decreasing to 1.1177, a reduction of approximately 20.7% compared to MOPSO-KELM; MAE decreased to 0.9514, a reduction of approximately 24.4%. Based on this, the proposed model was further optimized, with RMSE decreasing to 0.8439, a reduction of approximately 24.5% compared to the traditional GRU; MAE decreasing to 0.6893, a reduction of approximately 27.4%; and MAPE decreasing from 24.17% to 17.19%, a reduction of approximately 28.9%. In summary, the proposed model outperforms the traditional GRU and other comparative models in both prediction accuracy and fitting ability, demonstrating significant performance advantages.

[0339] In the Y direction of the lower guide bearing, the WUTP-KELM model optimized by WUTP achieved good prediction results, with an RMSE of 2.5756, a MAE of 2.126, a MAPE of 17.7587%, and R0. 2 The value was 0.9854. Compared to WUTP-KELM, the RMSE and MAE of the IWUTP-KELM model decreased slightly to 2.5654 and 2.1094, respectively. The MOPSO-KELM model performed further well, with its RMSE decreasing to 2.2714, a reduction of approximately 11.46% compared to IWUTP-KELM; MAE decreasing to 1.8532, a reduction of approximately 12.20%; and MAPE decreasing to 15.8731%. 2 The optimal value was achieved at 0.9878. Further incorporating the GRU model significantly improved prediction performance, with RMSE decreasing to 1.4638 (approximately 35.56% reduction compared to MOPSO-KELM), MAE decreasing to 1.2119 (approximately 34.57% reduction), and MAPE decreasing to 13.0246%. Compared to the traditional GRU model, the proposed model further improved performance, with RMSE decreasing from 1.4638 to 1.2228 (approximately 16.50% reduction), MAE decreasing from 1.2119 to 0.9996 (approximately 17.46% reduction), and MAPE decreasing from 13.0246% to 11.6514%.

[0340] In summary, the proposed model outperforms the traditional GRU and other comparative models in both prediction accuracy and fitting ability, demonstrating significant performance advantages.

Claims

1. A method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model, characterized in that... Includes the following steps: Step 1: Obtain the vibration signal and related operating characteristics of the pumped storage unit, construct a health status prediction model using the multinomial XGBoost algorithm, and introduce the SHAP method to evaluate the contribution of each input feature to the output of the health status prediction model. Step 2: Input the vibration data of the deteriorating working conditions to be predicted into the trained health state prediction model to obtain the corresponding predicted vibration data; A Predicted Degradation Index (PDI) is constructed based on the deviation between predicted and actual observed values ​​to quantify the degree of equipment degradation. Step 3: The constructed PDI sequence is decomposed using the IWUTP-optimized VMD method to obtain K intrinsic mode functions in order to obtain multi-scale temporal information of degradation features; Step 4: Use GRU to predict each intrinsic mode function separately, and superimpose the prediction results of each mode at each time step to construct the corresponding prediction training set sequence and prediction test set sequence; at the same time, retain the actual training set sequence of PDI as a reference for subsequent modeling. Step 5: Take the predicted training set sequence, predicted test set sequence, and actual training set sequence of PDI output from GRU as input to construct a KELM model based on IWUTP optimization, further model and predict the degradation trend, and finally output a more accurate degradation trend prediction result.

2. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 1, characterized in that: Step 1 includes the following steps: S1.1: Acquire the original vibration signal of the pumped storage unit and key operating characteristics such as power generation head, active power, and guide vane opening; S1.2: A health status prediction model is constructed based on the polynomial XGBoost algorithm. The calculation formulas included in the XGBoost algorithm are shown below: Let the training samples be M = (x1, y1), (x2, y2), ..., (x n ,y n ), i = 1, 2, 3, ..., n where x i Let y be the feature vector of the i-th sample. i For actual values, assume each decision tree is modeled. in, f(x) is the predicted value. i Let F represent the model's predicted output for the input feature x. Then, the overall objective function F is as follows: In formula (1): The loss function is represented by equation (2); Ω(f(x) t )) represents the regularization term, and its calculation formula is shown in equation (3); f(x) t ) indicates that the model adapts to the input feature x. t The predicted output, where t represents the time step and n represents the number of samples; Using the squared error as the loss function, its expression is as follows: A regularization term is introduced to control model complexity; its expression is shown below: In equation (3): γ represents the regularization coefficient used to control the number of numbers; T is the number of regression trees; λ represents the L2 regularization coefficient of the leaf node weights; π represents the leaf node weight vector of the regression tree, which contains the feature weights composed of the predicted values ​​output by all leaf nodes in the model; When constructing the health status prediction model, the power generation head, active power and guide vane opening are selected as the model input features, denoted as J=[j1,j2,j3], where j1, j2 and j3 represent the power generation head, active power and guide vane opening, respectively; S1.3: The SHAP method is introduced to evaluate the contribution of each input feature to the model output, thereby enabling the interpretability analysis of the model; the details are as follows: In a single decision tree, the SHAP value is calculated by analyzing the splitting of a feature across all decision paths. For a given feature, its SHAP value represents the expected marginal contribution of that feature to the predicted output across all possible paths, as expressed below: Where, x i Represents a feature; Shapley value is the feature x. i Contribution to the model output; Represents all cases containing x i The summation is performed on path p, where Pr(p) is the sampling probability or weight of path p; △f(p,x) i ) represents the feature x in path p. i The resulting increase in model output; Finally, the SHAP values ​​are summed to quantify the contribution of each feature to the overall prediction result of the model. In formula (6): SHAP j (g) represents feature x i global contribution, Shapley i In the nth sample, feature x i The Shapley value, where n represents the total number of samples.

3. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 2, characterized in that: In S1.2, a polynomial extension method is introduced based on the original XGBoost features. Specifically, based on the original XGBoost features, a power term and multidimensional interaction term are generated by extending the features using a fourth-order polynomial to describe the nonlinear degradation relationship and the interaction effects of multiple features; the specific formula is as follows: Where: J represents the original feature vector; a s Indicates the corresponding feature j s The polynomial power of ; φ(J) is the feature set after polynomial expansion; s represents the feature index, with values ​​ranging from 1, 2, 3; It is a set of non-negative integers.

4. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 3, characterized in that: Step 2 includes the following steps: S2.1: Input the power generation head, active power, guide vane opening and vibration data under adverse conditions into the trained health state prediction model, construct F(H(t), P(t), G(t)), where: H(t) is the power generation head, P(t) is the active power, and G(t) is the guide vane opening, and obtain the corresponding predicted vibration data; S2.2: The PDI is constructed based on the deviation between the predicted value and the actual observed value. The calculation formula is as follows: Wherein: F t V represents the vibration value predicted by the health model at time t. t These are the actual measured values ​​at the corresponding times.

5. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 4, characterized in that: Step 3 Includes the following steps: S3.1: The specific initialization steps of the WUTP algorithm are as follows: y i =[lb j +rand×(ub j -lb j )] (8); Where: y i Represents the i-th individual; lb j with ub j These represent the lower and upper bounds of the j-th dimension, respectively, and rand represents a random variable that follows a uniform distribution on the interval [0,1]. After initializing the population, the plant leaves and moving water particles were initialized separately; the initialization steps for the plant leaves are as follows: Leaves=[l1,l2,…,l i ,…,l m ] (9); Where: Leaves represents the collection of plant leaves; m represents the problem dimension, l i Let i be the i-th variable. The water particles in the plant are initialized as follows: Particles=[p1,p2,p3,…,p n ] (10); Where: Particles represents the set of moving water particles; n is the population size; the water particles must have sufficient hydrostatic pressure P, p1, p2, p3, ..., p n These represent individual water particles, driving their flow to a position close to the blades; at this point, the water particles are under a pressure gradient. Under the influence of the flow, the roots travel a distance x from the soil to the leaves; specifically, the plant roots travel a distance x from the soil to the leaves. i =0 absorbs water particles, and these particles travel a distance x within the plant leaf. f The flow distance can be determined by Calculated; In summary, the population can be represented as follows: Where: Totalop represents the set of the population; P 1,m P represents the hydrostatic pressure of the first water particle in the m-th dimension; 2,m P represents the hydrostatic pressure of the second water particle in the m-th dimension; n,m x represents the hydrostatic pressure of the nth water particle in the mth dimension; 1,m x represents the distance the first water particle travels in the m-th dimension; 2,m x represents the distance the second water particle travels in the m-th dimension; n,m y represents the distance the nth water particle travels in the mth dimension; 1,m This represents the first individual in m dimensions; y 2,m This represents the second individual in m dimensions; y n,m This represents the nth individual in m dimensions; S3.2: Calculate the fitness value based on the variable values ​​and store it in array f. i In this context, the fitness values ​​are as follows: Where: f1([y 1,1 y 1,2 … y 1,m [) represents the fitness set of the first water particle across all dimensions. [y 1,1 y 1,2 … y 1,m [] represents the set of individuals of the first water particle in each dimension; f2([y 2,1 y 2,2 … y 2,m [) represents the fitness set of the second water particle across all dimensions. [y 2,1 y 2,2 … y 2,m [] represents the set of individual water particles in each dimension; f n ([y n,1 y n,2 … y n,m [) represents the fitness set of the nth water particle across all dimensions. [y n,1 y n,2 … y n,m [] represents the set of individual water particles in each dimension; i = 1, 2, 3, ..., n represents the index of the i-th water particle in the population; j = 1, 2, 3, ..., m represents the dimension index of each particle, and the fitness value f of the i-th particle. i The definition is as follows: f i =f(y i,1 ,y i,2 ,…,y i,m ) (13); Where: f(y) i,1 ,y i,2 ,…,y i,m Let represent the fitness set of the i-th water particle across all dimensions; S3.3: The transport rate J of water particles in the soil, membrane structure, and root system is described by the following mathematical model: Where: t is the number of iterations. and Let represent the water flow velocity of the i-th water particle in the j-th dimension at the next iteration step and the current iteration step, respectively; σ represents the pressure gradient of the i-th particle in the j-th dimension; σ is the reflection coefficient, which takes values ​​in the range [0,1]. L represents the optimal pressure gradient of the i-th particle in the j-th dimension, obtained under the combined effects of randomness and systematicity, and is calculated by equation (16); p The hydraulic conductivity of the membrane is represented by r1 and r2, which are random numbers uniformly distributed in the interval [0, 1]. λ is the inertia coefficient used to regulate the amplitude of water particle flow, calculated by equation (18). g The gravitational term is represented and its calculation method is shown below: Where: ρ represents the density of water in the soil, roots, and plant membrane; g is the acceleration due to gravity; This represents the osmotic pressure gradient of the i-th particle in the j-th dimension, expressed by the formula... Calculations show that, where Π represents osmotic pressure. That is, the partial derivative of the i-th particle in the j-th dimension during the t-th iteration; Represents the hydrostatic pressure gradient. Let be the partial derivative of the i-th particle in the j-th dimension during the t-th iteration. Let be the flow distance of the i-th particle in the j-th dimension during the t-th iteration; The optimal pressure gradient of the i-th particle in the j-th dimension obtained under the combined effects of randomness and systematicity. The calculation is as follows: in: Let α be the optimal pressure gradient of the i-th particle in the j-th dimension, and let α be a random vector defined on the interval [0, 1]. The calculation process is as follows: Where: rand represents a uniformly distributed random value taking values ​​in the interval [0,1]. The value indicates rounding; the inertia coefficient λ used to regulate the amplitude of water particle flow is calculated as follows: S3.4: The hydrostatic pressure P drives the horizontal flow of water particles in the soil, and its velocity J can be described by the following formula: in: For the partial derivative of hydrostatic pressure, The partial derivative represents the distance water particles travel, where K is the hydraulic conductivity, with a value of 1 × 10⁻⁶. -9 When the water content changes, the water flow velocity under the influence of the water content gradient θ can be expressed as: Where: Q is the volumetric flux density of water in the soil, dl is the average root length per unit volume of soil, and D is the permeability coefficient; In summary, considering that the cross-sectional area of ​​the root system in the soil is πa 2 And dθ=πa 2 dP, where a is the radius of the root; based on actual conditions, the formula for calculating the water flow velocity J can be derived as follows: Based on the above formula The calculation expression is as follows: Where: p is a constant, c1, c2, and c3 are all random variables that follow a uniform distribution on the interval [0,1], and parameter a is set to 1; δ represents the water flow velocity of the i-th individual in the j-th dimension at iteration number t; rand Indicates the random decay factor; δ rand The definition is as follows: δ rand =v(lb j -(lb j -ub j )×rand) (24); Where: v represents the dynamic decay factor, γ is a constant, and t and T represent the current iteration number and the maximum iteration number, respectively; When water particles flow vertically, considering the work done against gravity, the total water potential φ is defined as follows: φ=ψ+ρgh=P-Π+ρgx (26);; Where: ψ is the water potential, x is the vertical distance of the water particle in the gravitational field, and ρgx represents the gravitational potential energy component; The vertical flow of water particles can be expressed by the following formula: According to dθ=πa 2 The relationship of dP can be derived from the following expression: Based on the above formula, under the action of vertical water flow, The calculation expression is as follows: Where: c4, c5, and c6 are all random variables that follow a uniform distribution on the interval [0,1]. S3.5: The mathematical model of the movement of water particles under the action of root absorption is shown below: Where: S(x,t) represents the root system's water absorption term. This function takes into account the root system's distribution characteristics at depth x and time t, and is used to characterize the process of plant roots absorbing water with time and soil depth. Considering the vertical migration of water under the influence of gravity, the motion of water particles can be described by the following formula: Wherein, the water potential ψ can be calculated from the difference between the hydrostatic pressure P and the osmotic pressure Π, i.e., ψ = P - Π. Based on this, the water flow velocity J can be expressed as: Therefore, the hydrostatic pressure gradient can be calculated using the following formula: Where: r5, r6, and c7 are independent random variables that are uniformly distributed on the interval [0,1]. S3.6: The kinetic model of water particle migration from the root surface to the xylem in the root surface to xylem stage is expressed as follows: Q=L p (ψ p,s -ψ p,x +x(ψ o,s -ψ o,x )) (34); Where: L p χ represents radial water conductivity, and ψ is the dimensionless reflection coefficient; p,s With ψ p,x ψ represents the hydrostatic pressure at the root surface and in the membrane, respectively. o,s With ψ o,x These represent the osmotic pressure at the root surface and within the membrane, respectively. Based on this, the hydrostatic pressure gradient can be expressed by the following equation: Where: χ takes the value of 0.5, r7, r8, and c8 are independent random variables that are uniformly distributed on the interval [0,1]. S3.7: The water pressure gradient from the xylem to the leaf is expressed by the following formula: Where: η represents the dynamic solution viscosity, a is the radius of the plant container, and c9 represents a random variable that follows a uniform distribution in the interval [0,1]. S3.8: The crossover operation first selects two individuals from the historical best individuals as parents, and performs crossover according to the crossover rate CR, the calculation formula of which is shown below: Where, r d Let z represent a random number that follows a uniform distribution on the interval [0,1]. d Let be the d-th dimension component of the offspring z, with a crossover rate CR of 0.7, and p i1 With p i2 These are two parent individuals selected from the historical best individuals; (p i2 ) d This represents the second parent individual selected from the historical best individuals in dimension d, (p i1 ) d Let represent the first parent individual selected from the historical best individuals in d dimensions; During the mutation phase, mutation operations are performed sequentially on each dimension of the individual, resulting in the mutated individual z. d The calculation formula is as follows: z d =lb d +u d ·(ub d -lb d )(h d <MR) (38); Where: u d with h d This represents a random number that follows a uniform distribution on the interval [0,1]. MR is the rate of variation, set to 0.

1. d with ub d These are the lower and upper bounds of the d-th dimension variable, respectively; after the mutation operation is completed, the fitness values ​​of the mutated individuals and their original individuals are compared, and individuals with better fitness are retained to update the population; S3.9: In the IWUTP algorithm, an elite strategy is used to record the individual with the best fitness in the current water particle swarm, and after each iteration, the individual with the worst fitness in the population is replaced by this elite individual; the mathematical expression for the extraction of elite individuals is as follows: Where: i * The index p represents the elite individual. i Let f(p) be the historical optimal position of the i-th particle, argmin represent the independent variable that minimizes the objective function, and N be the population size; i ) represents the fitness value of the i-th particle; After determining the optimal particle index, the corresponding optimal particle position and its fitness are calculated using the following formula: Where: f elite f(x) represents the optimal particle fitness value. elite ) represents the fitness value of an elite individual, x elite The individual with the best fitness in the current population. This indicates the position of elite individuals; simultaneously, the index of the worst-fit individual in the population is determined using the following mathematical expression: Where: j * f(x) represents the index of the worst-fit individual in the population. i ) represents the fitness of the i-th individual in the population; Let x i For the i-th individual in the current population, if Let represent the fitness value of the worst-fitting individual in the population; then, elite individuals are replaced, and the replacement formula is: in: Indicates the location of the individual with the worst fitness in the population; Furthermore, the current best individual is updated to the historical best solution using the following formula. Based on equations (8) to (42), the IWUTP algorithm is finally obtained.

6. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 5, characterized in that: It also includes S3.10: Variational Mode Decomposition (VMD) solves the following variational constraint problem through optimization, and its formula is shown below: In the above formula: K represents the number of modes to be decomposed, and k is the number of mode decompositions; u k (t) is the k-th mode; ω k δ(t) is the center frequency of the k-th mode; δ(t) is the Dirac function; It is a multiplication operator; u k Let f(t) represent the time-domain representation of the k-th mode, and let f(t) represent the original signal. S3.11: To optimize the solution of unconstrained variational problems, we introduce the Lagrange multiplier λ and the penalty factor α. The specific formulas are shown below: In equation (45): L({u k },{ω k },λ) denote the extended Lagrangian function, {u k } represents K modal functions, {ω k } represents the center frequencies of the K modes, and λ(t) is the Lagrange multiplier that varies with time; Among them, the modal component u k (t) and center frequency ω k The update formula is as follows: In the above formula, Let represent the frequency domain representation of the k-th modal component in the (n+1)-th iteration. This represents the center frequency of the k-th modal component in the (n+1)-th iteration. This is represented as the frequency domain representation of the i-th mode. Let f(t) and u represent respectively. k Frequency domain expressions for u(t) and u(t); ω represents the frequency; n represents the number of iterations; S3.12: The number of modes K and the penalty factor α of variational mode decomposition (VMD) are used as the two dimensions of the IWUTP algorithm, and the permutation entropy is used as the fitness function. Specifically, during the optimization process, each IWUTP individual corresponds to a set of (K, α) parameters. The population of the IWUTP algorithm is initialized by formulas (8) to (11). The algorithm will perform VMD decomposition based on the current parameters, calculate the permutation entropy of each mode, and use the minimum or average value of these entropy values ​​as the fitness evaluation index, thereby realizing the adaptive adjustment of parameters. The formula for permutation entropy is expressed as follows: In equation (48): nu represents the embedding dimension, l i Let represent the probability of the i-th permutation pattern, obtained statistically, and nu! represent the total number of all possible permutations.

7. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 6, characterized in that: In step 4, GRU is used to predict each intrinsic mode function, and the prediction results of each mode are superimposed point by point at each time step. Let the input sequence be X = (x1, x2, ..., x...). t ), x1, x2, ..., x t Input data at time steps 1, 2, ..., t respectively; The hidden state is H = (h1, h2, ..., h t ), h1, h2, ..., h t Let represent the hidden state vectors at time steps 1, 2, ..., t, respectively; The calculation and update process of GRU is defined by the following formula: r t =σ(W r ·[h t-1 ,x t ]+b r ) (49); z t =σ(W z ·[h t-1 ,x t ]+b z ) (50); Among them, z t and r t These represent updating the door and resetting the door, respectively; h t and Represent the hidden state and candidate hidden state, respectively; W and b are the weight matrix and bias vector of the corresponding terms, respectively; the symbol · denotes dot multiplication; t and t-1 represent the current time step and the previous time step, respectively; tanh is the hyperbolic tangent activation function; the symbol * denotes element-wise multiplication; h t-1 W represents the hidden state at the previous moment. r W represents the weight matrix input to the reset gate. z Let b represent the weight matrix input to the update gate. r b represents the bias vector input to the reset gate. z This represents the bias vector input to the update gate, [h] t-1 ,x t [h] indicates concatenating two vectors. t-1 *r t ,x t ] indicates that h t-1 With r t Multiply element by element, then multiply with vector x t The concatenation is performed, and σ(·) represents the Sigmoid activation function.

8. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 7, characterized in that: In step 4, the intrinsic mode functions obtained in step 3 are first input into the GRU model for training; the dataset is divided using a sliding window, specifically selecting the last 100 time steps as the test set and the rest as the training set, thus obtaining the input training set X. train Input test set X test and the actual training set sequence Y train ; During training, the Adam optimization algorithm is used to iteratively update the trainable parameters in the GRU model, including weights and biases, with the goal of minimizing the loss function, thereby improving the model's predictive performance. After training, the prediction training set sequence of the first-layer model is obtained. The test set from step 3 is then input into the trained GRU model to obtain the prediction test set sequence of the first-layer model. This allows us to construct the corresponding prediction training set sequence and prediction test set sequence, while retaining the actual training set sequence Y of PDI. train This will serve as a reference for subsequent modeling.

9. The method for predicting the deterioration trend of pumped storage units based on polynomial XGBoost and a two-layer prediction model according to claim 8, characterized in that: In step 5, for an ELM model with L hidden nodes, let the training sample set contain N samples, denoted as... Where: x i Let y be the input feature of the i-th sample. i Indicates the corresponding target output; (x i ,y i ) represents the data pair of the i-th training sample. Represents a g-dimensional real space. Represents a c-dimensional real space; The model aims to efficiently solve for the optimal output weight β by understanding the relationship between the hidden layer output and the output weights; therefore, the weight matrix between the hidden layer and the output layer can be expressed as: In equation (53): H represents the output matrix after the sample is mapped through the hidden layer, H T Let H be the transpose of H, T be the output target matrix, I be the unit diagonal matrix, and C be the regularization parameter. After introducing the kernel function, the resulting KELM model can map the input samples from the original n-dimensional space to a higher-dimensional feature space. The kernel matrix expression is defined as follows: In equation (54): Ω(i,j) represents the kernel matrix element, h(x i h(x) j ) are samples x i x j Mapping values ​​in the hidden layer; K(x i ,x j ) represents sample x i x j Kernel function; The Gaussian kernel function is chosen as the kernel function for the KELM model, and its expression is defined as follows: In equation (55): σ is the kernel parameter of the Gaussian kernel function; Based on the computational mechanism of the KELM model, the output of the model can be expressed as: In equation (56): f(x) represents the model's predicted output for the input sample x, h(x) represents the feature matrix of the input sample x after mapping, and K(x,x1) represents the kernel function matrix of samples x and x1. N ) represents the sample x and x N The kernel function matrix; The KELM model is optimized by the IWUTP algorithm obtained by embedding (8) to Equation (42) and used as the second-layer prediction model. Specifically, the root mean square error is used as the fitness function of the IWUTP algorithm to evaluate the regularization coefficient and kernel parameters in the KELM model, and the two hyperparameter values ​​corresponding to the minimum fitness function are selected as the final optimization results. The prediction training set, prediction test set and actual training set sequence of PDI output by GRU in step 4 are input together to finally obtain the actual PDI prediction results.

10. A deterioration trend prediction system for pumped storage units based on polynomial XGBoost and a two-layer prediction model, characterized in that: The system includes: a data acquisition module, a data generation module, a model training module, and a model prediction module, wherein: The data acquisition module acquires the original vibration data of the pumped storage unit, as well as key operating characteristics such as power generation head, active power, and guide vane opening. Historical data, which includes the original vibration data of the pumped storage unit, as well as power generation head, active power, and guide vane opening, will serve as key inputs for subsequent prediction models. The data generation module constructs a health status prediction model based on the multinomial XGBoost algorithm and introduces the SHAP method to evaluate the contribution of each input feature to the model output, thereby realizing the interpretability analysis of the model. The vibration data of the deteriorated working conditions to be predicted are input into the trained health status prediction model to obtain the corresponding predicted vibration data; PDI is constructed based on the deviation between the predicted value and the actual observed value to quantify the degree of equipment deterioration; the constructed PDI sequence is decomposed using the IWUTP optimized VMD method to obtain K intrinsic mode functions to obtain multi-scale time series information of deterioration characteristics. The model training module uses GRU to predict each intrinsic mode function separately, and superimposes the prediction results of each mode at each time step to construct the corresponding prediction training set sequence and prediction test set sequence; at the same time, it retains the actual training set sequence of PDI; the prediction training set, prediction test set output by GRU and the actual training set of PDI are used as input to construct a KELM model based on IWUTP optimization, further modeling and predicting the degradation trend, and finally outputting a more accurate degradation trend prediction result; the vibration trend prediction model of pumped storage unit based on GRU-IWUTP-KELM is trained. The model prediction module takes the PDI data of the pumped storage unit to be predicted and, after going through the same data preprocessing process, imports it as input data into the trained pumped storage unit vibration trend prediction model based on GRU-IWUTP-KELM to predict the deterioration trend of the pumped storage unit in the future time period.

Images