A flywheel health state evaluation and residual life prediction method and system

Abstract

This invention relates to a method and system for assessing the health status and predicting the remaining life of a flywheel, relating to the field of flywheel technology. The method includes: acquiring a flywheel vibration signal; extracting time-domain statistical features from the flywheel vibration signal to obtain a multidimensional feature vector; performing principal component analysis to reduce the dimensionality of the multidimensional feature vector to obtain a health status index; reconstructing the phase space of the health status index to obtain a phase space trajectory matrix; fitting the phase space trajectory matrix with a polynomial trend to obtain a degradation trajectory curve; extrapolating the degradation trajectory curve based on a preset failure threshold to obtain a predicted failure time point; and calculating the time difference between the predicted failure time point and the current time to obtain a predicted remaining life value. This invention solves the technical problem that traditional techniques struggle to accurately characterize the nonlinear evolution of flywheel performance degradation, resulting in low accuracy in remaining life prediction.

Description

Technical Field

[0001] This invention relates to the field of flywheel technology, and in particular to a method and system for assessing flywheel health status and predicting remaining lifespan. Background Technology

[0002] During the operation of flywheel energy storage systems, health status assessment and remaining life prediction are crucial for ensuring equipment safety and reliability. Current technologies typically employ vibration signal analysis to monitor the flywheel's operating status. Health indicators are constructed by extracting time-domain or frequency-domain characteristic parameters (such as root mean square value, kurtosis, and maxima), and then combined with empirical thresholds or simple regression models to provide a rough estimate of the remaining life. These methods are widely used in condition monitoring of industrial rotating machinery, particularly in wind turbines and high-speed motors, where they have already established a certain level of application.

[0003] However, the inventors found that the existing methods often lack effective dimensionality reduction and information fusion mechanisms when dealing with high-dimensional vibration characteristics, resulting in insufficient sensitivity of health indicators to early and subtle degradation characteristics. At the same time, in the lifetime prediction stage, they often rely on linear extrapolation or fixed failure thresholds, which makes it difficult to accurately characterize the nonlinear evolution law of flywheel performance degradation, resulting in low accuracy of remaining lifetime prediction and failing to meet the requirements of prediction timeliness and accuracy in high-reliability application scenarios. Summary of the Invention

[0004] The purpose of this invention is to at least partially solve one of the technical problems existing in the prior art.

[0005] To achieve the above objectives, the present invention provides a method for assessing the health status and predicting the remaining lifespan of a flywheel, comprising the following steps: The flywheel vibration signal is acquired, and time-domain statistical features are extracted from the flywheel vibration signal to obtain a multi-dimensional feature vector; Principal component analysis is performed on the multidimensional feature vectors to reduce their dimensionality, thereby obtaining health status indicators. The health status indicators are reconstructed in phase space to obtain a phase space trajectory matrix, and the phase space trajectory matrix is ​​fitted with a polynomial trend to obtain a degradation trajectory curve. The degradation trajectory curve is extrapolated based on a preset failure threshold to obtain the predicted failure time point. The time difference between the predicted failure time point and the current time is calculated to obtain the predicted remaining service life.

[0006] This invention also provides a flywheel health status assessment and remaining life prediction system, comprising: The extraction module is used to acquire the flywheel vibration signal, perform time-domain statistical feature extraction on the flywheel vibration signal, and obtain a multi-dimensional feature vector; The dimensionality reduction module is used to perform principal component analysis to reduce the dimensionality of the multidimensional feature vectors and obtain health status indicators. The reconstruction module is used to reconstruct the phase space of the health status indicators to obtain the phase space trajectory matrix, and to perform polynomial trend fitting on the phase space trajectory matrix to obtain the degradation trajectory curve. The calculation module is used to extrapolate the degradation trajectory curve based on a preset failure threshold to obtain the predicted failure time point, and calculate the time difference between the predicted failure time point and the current time to obtain the predicted remaining service life value.

[0007] This invention provides a method for assessing the health status and predicting the remaining life of a flywheel, comprising: acquiring a flywheel vibration signal; extracting time-domain statistical features from the flywheel vibration signal to obtain a multidimensional feature vector; performing principal component analysis to reduce the dimensionality of the multidimensional feature vector to obtain a health status index; reconstructing the phase space of the health status index to obtain a phase space trajectory matrix; performing polynomial trend fitting on the phase space trajectory matrix to obtain a degradation trajectory curve; extrapolating the degradation trajectory curve based on a preset failure threshold to obtain a predicted failure time point; and calculating the time difference between the predicted failure time point and the current time to obtain the remaining life. This approach uses lifetime prediction values ​​to address the technical problem of low remaining lifetime prediction accuracy caused by the inability of traditional techniques to accurately characterize the nonlinear evolution of flywheel performance degradation. The proposed method acquires flywheel vibration signals and extracts time-domain statistical features to form a multidimensional feature vector, effectively capturing the complex dynamic responses caused by fault mechanisms such as bearing wear, rotor imbalance, or magnetic circuit asymmetry during flywheel operation. Principal component analysis is then used to reduce the dimensionality of the high-dimensional features, eliminating redundancy and noise interference, and integrating information from multiple sensitive indicators to generate a health status index with monotonicity and trend characteristics, avoiding the vulnerability of single features to fluctuations in operating conditions. Based on this, phase space reconstruction is performed on the health status index, mapping the one-dimensional time series to a high-dimensional geometric trajectory, fully revealing the nonlinear dynamic structure of the flywheel degradation process. Finally, by performing polynomial trend fitting on the phase space trajectory matrix, a degradation trajectory curve that smoothly characterizes the performance evolution trend is constructed, overcoming the problem of drastic fluctuations in the original index and the difficulty in direct modeling. Furthermore, by combining a preset failure threshold with physical constraints, extrapolation calculations are performed on the degradation trajectory curve to accurately locate the predicted failure time point. A rigorous time base correction mechanism is then used to calculate the time difference between this point and the current time, ultimately outputting a highly reliable predicted remaining service life. The entire process tightly couples signal feature extraction, nonlinear state characterization, and trend extrapolation correction, significantly improving the ability to characterize the nonlinear evolution of flywheel performance degradation. Therefore, in this embodiment, the above scheme effectively solves the problem of low remaining service life prediction accuracy caused by neglecting the nonlinear characteristics of the degradation process in traditional technologies, achieving high-precision and robust prediction of the remaining service life of high-precision rotating machinery such as magnetic levitation flywheels. Attached Figure Description

[0008] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0009] Figure 1This is a schematic diagram of a flywheel health status assessment and remaining life prediction method in one embodiment of the present invention; Figure 2 This is the present invention. Figure 1 A schematic diagram of the implementation process of S1 in the middle; Figure 3 This is the present invention. Figure 1 A schematic diagram of the implementation process of S2 in the middle; Figure 4 This is a structural block diagram of a flywheel health status assessment and remaining life prediction system according to an embodiment of the present invention; The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0010] The embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention. The step numbers in the following embodiments are set only for ease of explanation, and there is no limitation on the order between the steps. The execution order of each step in the embodiments can be adaptively adjusted according to the understanding of those skilled in the art.

[0011] The following describes in detail, with reference to the accompanying drawings, a method for assessing the health status and predicting the remaining life of a flywheel according to an embodiment of the present invention. First, the method for assessing the health status and predicting the remaining life of a flywheel according to an embodiment of the present invention will be described in detail with reference to the accompanying drawings.

[0012] Figure 1 This invention provides a method for assessing the health status and predicting the remaining life of a flywheel, comprising the following steps: Step S1: Obtain the flywheel vibration signal, extract time-domain statistical features from the flywheel vibration signal to obtain a multi-dimensional feature vector, and perform principal component analysis to reduce the dimensionality of the multi-dimensional feature vector to obtain a health status index.

[0013] Specifically, after acquiring the flywheel vibration signal, the signal is first subjected to time-domain statistical feature extraction, including calculating conventional time-domain parameters such as mean, variance, peak value, kurtosis, skewness, and root mean square, thus constructing a feature vector containing multiple dimensions. Subsequently, principal component analysis (PCA) is performed on this multi-dimensional feature vector for dimensionality reduction. By constructing the covariance matrix and solving its eigenvalues ​​and eigenvectors, the top principal components with a cumulative contribution rate exceeding a preset threshold (e.g., 95%) are selected, projecting the original high-dimensional data into a low-dimensional space, ultimately obtaining a health status index that can characterize the overall operating state of the flywheel. For example, in a high-speed magnetic levitation flywheel energy storage system, the subtle vibration changes caused by early bearing wear or rotor imbalance are often submerged in high-dimensional redundant features. The aforementioned principal component analysis can effectively compress noise interference and highlight the degradation-sensitive direction. In this embodiment, the above scheme significantly improves the sensitivity and stability of the health status index to early performance degradation by fusing multi-dimensional time-domain features and using principal component analysis to achieve information condensation.

[0014] Step S2: Reconstruct the phase space of the health status index to obtain the phase space trajectory matrix, and perform polynomial trend fitting on the phase space trajectory matrix to obtain the degradation trajectory curve.

[0015] Specifically, when reconstructing the phase space of the health status indicators, the time delay τ and the embedding dimension m are first selected based on Takens' embedding theorem. τ can be determined by the first zero-crossing of the autocorrelation function or by mutual information, while m is determined by the spurious nearest neighbor method to determine the minimum effective embedding dimension. Then, the one-dimensional health status indicator sequence x(t) is used to construct a phase space trajectory matrix according to the formula X(t)=[x(t)x(t+τ)...x(t+(m-1)τ)]^T. Each row of this matrix represents an m-dimensional phase point, reflecting the overall evolution trajectory of the system in the reconstruction space. Next, polynomial trend fitting is performed on the principal component directions of this phase space trajectory matrix (usually the first principal component or the time series mean trajectory). For example, the least squares method is used to fit a quadratic or cubic polynomial to approximate the degradation trend of the flywheel performance over time, ultimately forming a continuous and smooth degradation trajectory curve. In the monitoring of high-speed magnetic levitation flywheel operation, this curve can effectively capture the nonlinear transition from the stable period to the accelerated degradation stage. In this embodiment, the above scheme overcomes the limitation that a single scalar index is difficult to characterize the dynamic degradation process by combining phase space reconstruction with polynomial fitting, making the degradation trajectory more consistent with the actual physical evolution law.

[0016] Step S3: Extrapolate the degradation trajectory curve based on the preset failure threshold to obtain the predicted failure time point, and calculate the time difference between the predicted failure time point and the current time to obtain the predicted remaining service life value.

[0017] Specifically, after obtaining the degradation trajectory curve, a preset failure threshold is introduced based on the trend of the health status index reflected by the curve over time. This threshold is determined in advance based on the historical operating data of the flywheel system or engineering safety specifications, and is used to define the critical point at which the equipment can no longer operate safely. The degradation trajectory curve is extrapolated along the time axis until it intersects with the failure threshold. The time coordinate corresponding to the intersection point is the predicted failure time point. This process can be achieved through interpolation or numerical approximation. Especially when the degradation trajectory changes non-linearly, an iterative search strategy is required to locate the time position of the first crossing of the threshold. Subsequently, the predicted failure time point is compared with the current actual operating time of the system, and the difference between the two constitutes the predicted remaining service life value. For example, in a high-speed magnetic levitation flywheel energy storage system, when the health status index continues to decline due to rotor imbalance or bearing degradation, and the extrapolation reaches the failure threshold, the failure window can be predicted in advance. In this embodiment, the above scheme effectively improves the interpretability and engineering practicality of the remaining service life prediction by combining the degradation trajectory with a clear failure boundary for time extrapolation.

[0018] In a specific embodiment, such as Figure 2 As shown, time-domain statistical feature extraction is performed on the flywheel vibration signal to obtain a multi-dimensional feature vector, including: S11, after taking the absolute value of the flywheel vibration signal point by point, calculate the sequence average to obtain the average rectified value, and calculate the peak-to-peak value by subtracting the minimum value from the maximum value of the flywheel vibration signal. S12, the flywheel vibration signal is zero-point offset corrected based on the average rectified value to obtain a centered sequence, and the square root average of the centered sequence is obtained by summing the squares based on the peak-to-peak value to obtain the root mean square value. S13, based on the average rectified value, the peak-to-peak value and the root mean square value, a multidimensional feature vector is obtained by numerical concatenation.

[0019] Specifically, when extracting time-domain statistical features from the flywheel vibration signal, the absolute value of each sampling point of the signal is first taken, and then the arithmetic mean of the resulting absolute value sequence is calculated to obtain the average rectified value. Simultaneously, the maximum and minimum amplitude values ​​are identified from the original flywheel vibration signal, and the difference between them is used to obtain the peak-to-peak value. Subsequently, the original flywheel vibration signal is zero-point offset corrected using the aforementioned average rectified value, i.e., the average rectified value is subtracted from the entire signal to form a centered sequence. Based on this, using the peak-to-peak value as a normalization reference scale, the centered sequence is squared point by point, accumulated, divided by the data length, and the square root is taken to finally obtain the root mean square value. The average rectified value, peak-to-peak value, and root mean square value are concatenated in a fixed order to form a multidimensional feature vector. During the operation of the high-speed magnetic levitation flywheel, even slight imbalances or local bearing defects will cause coordinated changes in these statistical quantities, thus effectively capturing them. In this embodiment, the above scheme integrates rectification characteristics, dynamic range, and energy measurement to construct a multi-dimensional feature vector that is sensitive to early degradation and computationally stable, providing a reliable input for subsequent health status assessment.

[0020] In a specific embodiment, time-domain statistical feature extraction is performed on the flywheel vibration signal to obtain a multi-dimensional feature vector, including: The flywheel vibration signal is divided into sampling points to obtain multiple sampling point sequences; The mean of each sampling point sequence is calculated to obtain the mean of each sequence, and the variance of each sampling point sequence is calculated to obtain the variance of each sequence. The mean reflects the average vibration amplitude of the flywheel vibration signal in each sampling point sequence, and the variance reflects the fluctuation of the vibration amplitude of the flywheel vibration signal in each sampling point sequence. Based on the mean and variance of each sequence, the peak value of the flywheel vibration signal is calculated to obtain the peak value of each sampling point sequence, wherein the peak value reflects the maximum vibration amplitude of the flywheel vibration signal at each sampling point sequence; The mean, variance, and peak values ​​of each sequence are combined to obtain a multidimensional feature vector.

[0021] Specifically, when extracting time-domain statistical features from the flywheel vibration signal, the continuously acquired vibration signal is first divided into several segments of fixed length, each segment containing the same number of sampling points, thus forming multiple sampling point sequences. Then, the arithmetic mean of each sampling point sequence is calculated to obtain the mean of each sequence, which represents the average vibration amplitude of the flywheel vibration signal within that time period. Simultaneously, the variance of each sampling point sequence is calculated to quantify the dispersion of the vibration amplitude around the mean, i.e., the fluctuation of the vibration amplitude. Based on this, the maximum absolute value of the same sampling point sequence is further extracted as the peak value of that segment, reflecting the local maximum vibration intensity. Finally, the mean, variance, and peak values ​​of all sampling point sequences are sequentially concatenated in chronological or channel order to form a high-dimensional multidimensional feature vector. For example, during the acceleration of the magnetic levitation flywheel to its rated speed, if the variance of a certain sequence suddenly increases while the peak value rises synchronously, it may indicate a transient imbalance in the rotor. In this embodiment, the above scheme constructs a multi-dimensional feature vector with time-series resolution by jointly characterizing the central trend, discrete characteristics and extreme response of the vibration signal, which significantly enhances the ability to identify early degradation modes.

[0022] In a specific embodiment, such as Figure 3 As shown, the step of performing principal component analysis to reduce the dimensionality of the multidimensional feature vector to obtain health status indicators includes: S21, perform data standardization processing on the multidimensional feature vector to obtain standardized feature vector, and calculate the covariance between each pair of features in the standardized feature vector to obtain the covariance matrix; S22, perform eigenvalue decomposition on the covariance matrix to obtain eigenvalues ​​and corresponding eigenvectors, and calculate the proportion of each eigenvalue to the sum of all eigenvalues ​​to obtain the variance contribution rate of each principal component. According to the order of variance contribution rate from large to small, select the eigenvectors corresponding to the first few eigenvalues ​​whose cumulative variance contribution rate reaches the preset contribution rate threshold to form a projection matrix. S23, the standardized feature vector is linearly transformed using the projection matrix to obtain a low-dimensional principal component vector. The low-dimensional principal component vector is then weighted and summed based on the variance contribution rate of each principal component to obtain a health status index.

[0023] Specifically, when performing principal component analysis (PCA) to reduce the dimensionality of the multidimensional feature vectors to obtain health status indicators, the original multidimensional feature vectors first need to undergo data standardization. This process is performed independently for each feature dimension. Specifically, the mean of all sample values ​​for that dimension is subtracted, and then the result is divided by its standard deviation. This eliminates the influence of differences in units or amplitude scales among different features, ultimately forming a standardized feature vector. Subsequently, the covariance between any two features is calculated based on this standardized feature vector, constructing a symmetric covariance matrix. This matrix fully characterizes the linear correlation structure between the features, providing a mathematical basis for subsequent dimensionality reduction. Next, eigenvalue decomposition is performed on the covariance matrix to obtain a set of real eigenvalues ​​and their corresponding unit eigenvectors, where the magnitude of each eigenvalue reflects the amount of data variance carried by the corresponding principal component. Furthermore, all eigenvalues ​​are summed, and the proportion of each eigenvalue to the total sum is calculated, i.e., the variance contribution rate of each principal component. These eigenvalues ​​are then sorted from largest to smallest according to this proportion, and accumulated until the cumulative variance contribution rate reaches a preset threshold (e.g., 95%). At this point, the eigenvectors corresponding to the selected k largest eigenvalues ​​are combined into a projection matrix. The column vectors of this projection matrix represent the retained principal component directions. Then, the aforementioned standardized eigenvectors are multiplied by this projection matrix to perform a linear transformation, outputting a low-dimensional principal component vector with significantly reduced dimensionality. It is worth noting that although the principal components are orthogonal to each other and sorted by variance after dimensionality reduction, directly using the first principal component sometimes fails to fully reflect the equipment degradation trend, especially when secondary components contain sensitive fault information. Therefore, in this scheme, instead of only using the first principal component, the variance contribution rate of each principal component is used as a weighting coefficient to perform a weighted summation of each component in the low-dimensional principal component vector, ultimately fusing them to generate a single scalar—the health status index. Taking a high-speed magnetic levitation flywheel system as an example, its multidimensional feature vector may contain dozens of statistics. After principal component analysis, if the cumulative contribution rate of the top three principal components reaches 96.2%, then their variance contribution rates of 0.72, 0.18, and 0.06 are used as weights to weight the scores of the three principal components. The resulting analysis can more robustly track performance degradation caused by bearing wear or rotor imbalance, avoiding monitoring blind spots caused by ignoring minor but sensitive components. In this embodiment, the above scheme introduces a variance contribution rate weighting mechanism, which retains the main information while taking into account potential sensitive degradation characteristics. This enables the health status indicators to not only have noise reduction and compression capabilities, but also effectively characterize the performance evolution trend of the flywheel system throughout its entire life cycle, significantly improving the continuity of condition monitoring and the reliability of prediction.

[0024] In a specific embodiment, the step of reconstructing the phase space of the health status indicators to obtain the phase space trajectory matrix includes: The health status index is segmented to obtain multiple degradation stage data segments, and an autocorrelation function is calculated for each degradation stage data segment to obtain an autocorrelation function value sequence, wherein the autocorrelation function value sequence includes the autocorrelation function values ​​corresponding to each delay time. The initial delay time for each degradation stage is determined based on the delay time corresponding to the first drop of the autocorrelation function value to a preset threshold. For each data segment in the degradation stage, starting from the initial delay time, the delay time is gradually increased, and false nearest neighbor analysis is performed at each delay time. The proportion of false nearest neighbors is calculated. The proportion of false nearest neighbors is obtained by calculating the ratio of the number of false nearest neighbors in the phase space to the total number of neighbors under a given embedding dimension. Based on the delay time and embedding dimension corresponding to the first drop of the proportion of false nearest neighbors to a preset proportion threshold, the optimal delay time and embedding dimension for each degradation stage are determined. Based on the optimal delay time and embedding dimension of each degradation stage, the phase space of the data segments of each degradation stage is reconstructed to obtain the phase space trajectory sub-matrix of each degradation stage. The phase space trajectory sub-matrixes of each degradation stage are then concatenated in chronological order to obtain the phase space trajectory matrix.

[0025] Specifically, when reconstructing the phase space of the health status indicators to obtain the phase space trajectory matrix, the continuous health status indicator time series is first divided into multiple degradation stage data segments according to the performance degradation characteristics during actual equipment operation. Each data segment corresponds to the evolution behavior of the magnetic levitation flywheel system within a specific lifespan. This segmentation is not a simple equal-length division, but rather based on indicator trend changes, operating condition switching, or historical fault records to ensure that the internal dynamic characteristics of each segment are relatively consistent.

[0026] Subsequently, for each data segment in the degradation stage, its autocorrelation function value sequence is calculated. Specifically, starting from zero and gradually increasing the delay time, the normalized covariance between the data segment and its own delayed sequence is calculated for each delay time, thus obtaining a series of autocorrelation function values. As the delay time increases, the autocorrelation function value typically shows a decreasing trend. When this value first drops to a preset threshold (e.g., 0.1), the corresponding delay time is determined as the initial delay time for that degradation stage. This initial delay time serves only as the starting point for subsequent parameter optimization and is not directly used for the final reconstruction.

[0027] Next, starting from the initial delay time, the delay time is gradually increased, and a false nearest neighbor analysis is performed at each candidate delay time. This analysis requires setting an initial embedding dimension (e.g., 3). The data segment of the current degradation stage is initially embedded into the phase space according to the given embedding dimension and the current delay time, forming a trajectory point set. For each trajectory point, its nearest neighbor at the current dimension is found; then the embedding dimension is increased by one dimension. If the distance between two points in the new dimension increases significantly (usually more than 10 times the original distance is used as the criterion), it is determined to be a false nearest neighbor. The ratio of all such false nearest neighbors to the total number of neighbors is calculated, which is the false nearest neighbor ratio. If this ratio cannot be reduced to a preset ratio threshold (e.g., 5%) at the current embedding dimension, the embedding dimension needs to be increased and the above process repeated. When the false nearest neighbor ratio first drops below this threshold, the corresponding delay time and embedding dimension are determined as the optimal delay time and embedding dimension for this degradation stage.

[0028] It is worth noting that the dynamic complexity varies across different degradation stages. For example, in the early micro-wear stage of a flywheel bearing, health indicators change gradually, autocorrelation decays slowly, the initial delay time is relatively large, and a high embedding dimension is required to fully expand the phase space. However, in the near-failure stage, indicators fluctuate dramatically, autocorrelation decays rapidly, the optimal delay time shortens, and the embedding dimension may actually decrease. Therefore, determining parameters independently at each stage can more accurately match the local dynamic structure.

[0029] After obtaining the optimal delay time and embedding dimension for each degradation stage, phase space reconstruction is performed on the data segments of each degradation stage. Specifically, based on the optimal parameters, a one-dimensional health status index sequence is mapped to multi-dimensional trajectory points using a sliding window method. Each point consists of the index values ​​at the current time and several delayed time points, thus forming the phase space trajectory sub-matrix for that stage. Finally, the phase space trajectory sub-matrices of all degradation stages are concatenated end-to-end strictly according to the original time order to form a complete phase space trajectory matrix.

[0030] In this embodiment, the above scheme effectively avoids the phase space folding or redundant expansion problem caused by using globally unified parameters by determining the delay time and embedding dimension in stages and combining the initial selection of autocorrelation function and the fine adjustment mechanism of false nearest neighbor ratio. This enables the reconstructed phase space trajectory matrix to reflect the nonlinear dynamic characteristics of the magnetic levitation flywheel system in different degradation stages throughout its entire life cycle in a true and complete manner, providing a high-fidelity state characterization basis for subsequent degradation pattern recognition and remaining lifetime prediction.

[0031] In a specific embodiment, the step of performing polynomial trend fitting on the phase space trajectory matrix to obtain the degenerate trajectory curve includes: The phase space trajectory matrix is ​​segmented to obtain multiple trajectory segment matrices, and the mean of each trajectory segment matrix is ​​calculated to obtain the mean vector of each trajectory segment. Based on the mean vector of each trajectory segment, a piecewise polynomial fitting is performed on the phase space trajectory matrix to obtain the polynomial fitting curve of each trajectory segment. Based on the start and end time points of each trajectory segment, the polynomial fitting curves of each trajectory segment are spliced ​​together to obtain the preliminary degenerate trajectory curve after splicing. The initial degradation trajectory curve is smoothed by using the sliding window averaging method to perform smoothing calculations on each point on the initial degradation trajectory curve, resulting in a smoothed degradation trajectory curve.

[0032] Specifically, when performing polynomial trend fitting on the phase space trajectory matrix to obtain the degradation trajectory curve, the matrix first needs to be segmented along the time axis. Specifically, the entire phase space trajectory matrix is ​​divided into several continuous and non-overlapping trajectory segment matrices according to a preset time window length or based on the degradation rate change points. Each trajectory segment matrix corresponds to the phase space evolution state of the flywheel system during a certain operating time. For example, during the long-term operation of a magnetic levitation flywheel, degradation is slow in the early stage, accelerates in the middle stage, and abruptly changes in the late stage. Based on this, the trajectory matrix can be divided into three trajectory segment matrices, each corresponding to a different degradation rate interval.

[0033] Subsequently, the mean is calculated for each trajectory segment matrix. This operation does not calculate the global mean over the entire matrix, but rather takes the column-wise arithmetic mean of all row vectors (i.e., trajectory points in phase space) within each trajectory segment matrix, resulting in a mean vector with the same dimension as the embedding dimension. This mean vector can be considered the "centroid" of the trajectory segment in phase space, representing the central tendency of the system state within that time period. For example, if a trajectory segment matrix consists of 500 6-dimensional trajectory points, its mean vector is a 6-dimensional vector, where each dimension represents the average value of the corresponding coordinate components.

[0034] Based on this, a piecewise polynomial fitting is performed on the phase space trajectory matrix according to the time center point corresponding to the mean vector of each trajectory segment. In practice, instead of directly fitting the high-dimensional matrix as a whole, the phase space trajectory is projected onto a dominant direction (such as the direction of the first principal component) or a key dimension (such as the reconstructed component corresponding to the original dimension of the health status indicator) and scalarized into a one-dimensional time series. Then, using the start time, end time, and center time of each trajectory segment as independent variables, and the projection value of the corresponding mean vector in the scalar direction as the dependent variable, a low-order polynomial (usually quadratic or cubic) is independently fitted to each trajectory segment using the least squares method, thereby obtaining the polynomial fitting curve of each trajectory segment. Each curve is only valid within its corresponding time period, ensuring accurate capture of local degradation trends.

[0035] Next, based on the start and end times of each trajectory segment, the polynomial fitting curves of adjacent trajectory segments are spliced ​​together at the time boundaries. Temporal continuity is maintained during splicing, meaning the end time of the previous segment is strictly aligned with the start time of the next segment, avoiding jumps or overlaps. This forms a preliminary degradation trajectory curve covering the entire life cycle. Although this curve reflects the overall degradation trend, due to independent fitting of each segment, discontinuities or slight fluctuations in derivatives may occur at the connection points between segments, manifesting as local "bends" or jitters.

[0036] To eliminate such discontinuities, the initial degradation trajectory curve needs to be smoothed. Specifically, a sliding window averaging method is used: a fixed-length window (e.g., containing 11 consecutive time points) is set and slid along the time axis. The arithmetic mean of the function values ​​of all curve points within the window is taken, and this average is assigned to the center point of the window, thus generating a new smoothing point. The choice of window length must balance noise suppression and trend fidelity—too long will blur the degradation inflection points, while too short will result in insufficient smoothing. After this processing, high-frequency disturbances on the original initial degradation trajectory curve are effectively filtered out, and the curve as a whole presents a smooth, continuous monotonic or non-monotonic degradation pattern, which better conforms to the actual evolution law of the physical system.

[0037] In this embodiment, the above-described scheme achieves a robust mapping from a high-dimensional phase space trajectory to a one-dimensional degradation trajectory curve through steps such as segmenting the phase space trajectory matrix, extracting the mean vector, piecewise polynomial fitting, and sliding window smoothing. This method preserves the local dynamic characteristics of different degradation stages and eliminates discontinuities at the fitting boundaries through smoothing, resulting in a degradation trajectory curve with good smoothness, continuity, and physical interpretability, providing a stable and reliable input benchmark for subsequent remaining lifetime prediction models.

[0038] In a specific embodiment, the step of extrapolating the degradation trajectory curve based on a preset failure threshold to obtain the predicted failure time point includes: The degraded trajectory curve is subjected to polynomial fitting parameter inversion to obtain trajectory polynomial coefficients, and the time axis extension calculation is performed on the degraded trajectory curve based on the trajectory polynomial coefficients to obtain the future predicted trajectory sequence. The predicted trajectory sequence is numerically compared with a preset failure threshold to obtain the predicted intersection time point, and the second derivative of the trajectory polynomial coefficients at the predicted intersection time point is calculated to obtain the trajectory divergence curvature. Sensitivity correction calculations are performed on the predicted intersection time points based on the trajectory divergence curvature to obtain failure time correction values, and the predicted failure time points are determined based on the failure time correction values.

[0039] Specifically, when extrapolating the degradation trajectory curve to determine the predicted failure time point, the first step is to perform parametric inversion on the polynomial expression corresponding to the curve. Since the degradation trajectory curve has already been obtained through piecewise polynomial fitting and smoothing in the previous steps, its overall form can be approximated by a unified low-order polynomial (such as cubic or quartic). Parametric inversion utilizes the least squares method, taking all time-state value pairs on the smoothed degradation trajectory curve as observation data, and solving for the coefficients of each term of the polynomial, thus obtaining the complete trajectory polynomial coefficients. These coefficients uniquely determine the mathematical expression of the degradation trajectory, providing an analytical basis for subsequent timeline extension.

[0040] After obtaining the trajectory polynomial coefficients, they are substituted into the polynomial function, and the independent variable (time) is gradually extended from the current maximum observation time to several future periods to calculate the corresponding state prediction values, thereby generating a future predicted trajectory sequence. This sequence reflects the evolution trend of the system's health status indicators under no-intervention conditions. For example, in the application of a magnetic levitation flywheel, if it is currently running at 8000 hours, the time can be extrapolated to 10000 hours, with a step size of 10 hours, resulting in 200 future state prediction points, which constitute the future predicted trajectory sequence.

[0041] Subsequently, each predicted value in the future trajectory sequence is compared point-by-point with a preset failure threshold. The failure threshold is typically set based on equipment safety margins, historical fault data, or industry standards; for example, a health status indicator dropping to 0.3 is considered a functional failure. During the comparison, once a predicted point is found to be below (or above, depending on the indicator definition) the failure threshold for the first time, its corresponding time point is recorded as the predicted intersection time point. However, because polynomial extrapolation may exhibit non-physical oscillations or excessively fast / slow decay when far from the fitting interval, relying solely on a single intersection point can easily introduce significant errors.

[0042] Therefore, it is necessary to further calculate the second derivative of the trajectory polynomial coefficients at the predicted intersection time point; this value is the trajectory divergence curvature. The second derivative reflects the trend of the degradation rate: if it is negative and the absolute value is large, it indicates that the degradation is accelerating; if it is close to zero, the degradation tends to be stable. In the flywheel bearing wear scenario, the near-failure stage often manifests as a sharp decline in health indicators, at which point the second derivative is significantly negative. The magnitude of the trajectory divergence curvature is directly related to the reliability of the extrapolation results—the larger the curvature, the more drastic the trajectory changes in that region, and the more sensitive the original predicted intersection time point is to model perturbations.

[0043] Based on this, a sensitivity correction calculation is introduced: the predicted intersection time point is fine-tuned using the trajectory divergence curvature as a weighting factor. Specifically, if the absolute value of the trajectory divergence curvature exceeds a certain empirical threshold (e.g., 0.001), the predicted intersection time point is scaled back by a certain proportion (e.g., 5%~10%) towards the current time to compensate for the uncertainty of the extrapolation model in the high curvature region; conversely, if the curvature is small, the original predicted intersection time point is retained or only slightly adjusted. The final result is the failure time correction value.

[0044] This failure time correction value is then determined as the predicted failure time point to guide maintenance decisions. The entire process is not a simple linear extrapolation, but rather a multi-level calibration mechanism that integrates polynomial analytical expression, threshold cross-detection, and curvature sensitivity correction.

[0045] In this embodiment, the above scheme achieves analytical extrapolation through trajectory polynomial coefficient inversion, combines failure threshold cross-detection to locate the initial failure time, and introduces trajectory divergence curvature to dynamically correct the prediction results, effectively suppressing the instability of high-order polynomial extrapolation in long-term prediction. This method significantly improves the robustness and engineering practicality of failure time estimation in the remaining life prediction of high-reliability rotating machinery such as magnetic levitation flywheels, avoiding premature or late warnings caused by ignoring degradation acceleration changes.

[0046] In a specific embodiment, the step of performing polynomial fitting parameter inversion on the degraded trajectory curve to obtain trajectory polynomial coefficients includes: The degradation trajectory curve is discretized over time to obtain a time sampling sequence, and the time sampling sequence is calculated by exponential increment to obtain the Vandermonde matrix; The column vectors of the Vandermonde matrix are orthogonalized to obtain an orthogonal coefficient matrix. Based on the orthogonal coefficient matrix, the Vandermonde matrix is ​​transformed to obtain a condition number optimization matrix. The condition number optimization matrix is ​​solved by least squares to obtain the intermediate polynomial coefficients, and the intermediate polynomial coefficients are then subjected to higher-order term damping weighting to obtain the trajectory polynomial coefficients.

[0047] Specifically, when performing polynomial fitting parameter inversion on the degradation trajectory curve to obtain the trajectory polynomial coefficients, the degradation trajectory curve first needs to be subjected to time discretization sampling. This operation is not a simple equal-interval truncation, but rather, based on the timestamps of the original health status indicators, N time points are uniformly extracted within their effective intervals to form a time sampling sequence. For example, in a magnetic levitation flywheel system, if the degradation trajectory covers 8000 to 8500 hours, samples can be taken every 10 hours, resulting in 51 sampling points. Subsequently, the time sampling sequence is calculated exponentially, that is, for each time point... Calculate in sequence ,in Given a preset polynomial order (e.g., 5), construct a... Vandermonde matrix , its first Behavior This matrix directly links the coefficients to be determined with the observed values ​​and forms the basis for subsequent parameter calculations.

[0048] However, high-order Vandermonde matrices are highly susceptible to severe ill-conditioning at large time scales, with condition numbers potentially reaching as high as [missing value]. The magnitude of the noise makes the least squares solution extremely sensitive to measurement noise. To alleviate this problem, the Vandermonde matrix needs to be orthogonalized in its column vectors. Specifically, a modified Gram-Schmidt orthogonalization algorithm is used: starting from the first column, each column is projected onto the orthogonalized subspace and the projected components are subtracted, ultimately obtaining a set of orthonormal bases, which constitute the orthogonal coefficient matrix. This process not only preserves the column space of the original matrix but also significantly improves numerical stability. Next, a matrix transformation is performed on the Vandermonde matrix based on the orthogonal coefficient matrix, i.e., the calculation... The result is the condition number optimization matrix. Since... Since it is an orthogonal matrix, this transformation does not change the solution of the least squares problem, but it transforms the original ill-conditioned system into a well-conditioned triangular system, which greatly reduces the solution error.

[0049] Based on this, the condition number optimization matrix is ​​solved using least squares. Let the sampled values ​​of the health status indicators corresponding to the degradation trajectory be a vector. Then by solving The intermediate polynomial coefficients can be obtained. While the coefficients are numerically stable, higher-order terms (such as 4th and 5th order) can still cause non-physical oscillations due to overfitting, especially leading to severe fluctuations during extrapolation at the end of the trajectory. Therefore, it is necessary to apply damped weighting to the intermediate polynomial coefficients for higher-order terms. Specifically, this is achieved by introducing a weight vector. The lower-order terms have a weight of 1, and the higher-order terms have a weight of 1. (For example, take 0.6, 0.4, 0.2), and then multiply them element by element to obtain the final trajectory polynomial coefficients. This weighting strategy suppresses the dominant role of higher-order terms, enabling the fitted curve to maintain its trend while possessing smooth extrapolation capabilities.

[0050] In this embodiment, the above scheme generates a time sampling sequence by discretizing the degradation trajectory curve over time. After constructing the Vandermonde matrix, column vector orthogonalization and matrix transformation are performed to obtain the condition number optimization matrix. Then, through least squares solution and higher-order term damping weighting, the final output is a numerically stable and physically reasonable trajectory polynomial coefficient. This method effectively overcomes the ill-conditioned and overfitting problems of traditional polynomial fitting in higher-order cases, providing a reliable model foundation for subsequent extrapolation prediction. It significantly improves the robustness and extrapolation consistency of degradation trajectory modeling in the prediction of the remaining life of magnetic levitation flywheels.

[0051] In a specific embodiment, the step of performing time-axis extension calculations on the degraded trajectory curve based on the trajectory polynomial coefficients to obtain a future predicted trajectory sequence includes: Based on the trajectory polynomial coefficients, a polynomial evaluation calculation is performed on the preset future time series to obtain an initial extrapolation scatter plot. Then, a first-order difference calculation is performed on the adjacent data points of the initial extrapolation scatter plot to obtain an extrapolation slope curve. The extrapolation slope curve includes the future time series generated by extending backward from the end time point of the degenerate trajectory curve by a fixed step size. The extrapolated slope curve is subjected to extreme value limiting and negative value zeroing calculations to obtain a monotonic slope envelope. Based on the monotonic slope envelope, the terminal values ​​of the degenerate trajectory curve are accumulated point by point to obtain the future predicted trajectory sequence.

[0052] Specifically, when performing time-axis extension calculations on the degenerate trajectory curve based on the trajectory polynomial coefficients to obtain the future predicted trajectory sequence, a preset future time series must first be constructed. This future time series is not arbitrarily set, but rather starts from the end time point of the degenerate trajectory curve and extends forward several cycles at fixed step sizes (e.g., every 5 or 10 hours) to form an equally spaced time vector. For example, if the observation ends after the magnetic levitation flywheel system has been running for 8500 hours, and prediction is needed up to 10000 hours, then starting from 8500, with a step size of 10 hours, 150 future time points are generated: 8510, 8520, ..., 10000, constituting the future time series.

[0053] Subsequently, the coefficients of the trajectory polynomial are substituted into the standard polynomial expression, and polynomial evaluation is performed for each time point in the future time series to obtain a set of corresponding predicted health status indicators. These predicted values ​​form a discrete set of points on the time-state plane, i.e., an initial extrapolation scatter plot. Because higher-order polynomials may produce non-physical oscillations or anomalous sharp drops in the extrapolation region, although the scatter plot is mathematically continuous, its local trends may not conform to actual degradation patterns. Therefore, further analysis of its dynamic characteristics is necessary.

[0054] The specific operation involves performing a first-order difference calculation on adjacent data points in the initial extrapolated scatter plot. Let the first difference be... The state value of each predicted point is The corresponding time is Then the first difference between it and the previous point is divided into This results in an extrapolated slope curve reflecting the change in the degradation rate. This curve visually presents the instantaneous degradation rate at each future moment, but in the region dominated by higher-order terms, there may be alternating positive and negative values ​​or a maximum peak value, which contradicts the actual monotonic degradation behavior of the magnetic levitation flywheel system.

[0055] To correct such unreasonable fluctuations, physical constraints need to be imposed on the extrapolated slope curve. First, the maximum value of the first derivative of the degradation trajectory curve over the entire historical observation period is extracted (obtained through numerical differentiation or analytical differentiation), and this is used as the upper limit threshold for extreme value limiting. For example, if the historical maximum degradation rate is a decrease of 0.002 units per hour, any portion of the extrapolated slope exceeding this value is clipped to this upper limit. Second, considering that the health status index of the magnetic levitation flywheel will not exhibit a "restorative rise" during normal degradation, all negative slopes (i.e., index rebound) have no physical meaning, so they are forcibly set to zero. After this dual processing of extreme value limiting and setting negative values ​​to zero, the result is a monotonic slope envelope—a non-negative, bounded, and corrected slope sequence that conforms to the upper limit of the historical degradation rate.

[0056] Finally, the terminal values ​​of the degradation trajectory curve are calculated point-by-point based on the monotonic slope envelope. Specifically, the actual state values ​​of the degradation trajectory curve at the terminal time point are used. Using the initial value, the slopes of each element in the monotonic slope envelope are multiplied by their corresponding time step, and these multiplications are accumulated and added to the predicted value of the previous time step, thus recursively generating a new state prediction sequence. For example, if the step size is 10 hours and the first correction slope is 0.0015, then the next prediction point is... And so on. This process ensures that the future predicted trajectory sequence not only inherits the historical degradation trend, but also strictly satisfies the engineering prior knowledge of monotonically non-increasing and rate-bounded.

[0057] In this embodiment, the above scheme obtains an initial extrapolation scatter plot through polynomial calculation, and then constructs a monotonic slope envelope by combining first-order difference, extreme value limiting and negative value zeroing, and uses this to drive point-by-point accumulation to generate a future predicted trajectory sequence. This effectively suppresses the non-physical understanding brought about by pure mathematical extrapolation, makes the prediction results more consistent with the actual degradation behavior of the magnetic levitation flywheel system, and significantly improves the credibility and engineering applicability of the remaining life prediction.

[0058] In a specific embodiment, calculating the time difference between the predicted failure time point and the current time to obtain the predicted remaining useful life includes: The difference between the predicted failure time and the current time is calculated to obtain the preliminary remaining lifetime, and the optimal delay time is accumulated to obtain the reconstruction lag time. Based on the reconstruction lag time, the preliminary remaining lifetime time is corrected using a time base to obtain a time anchor offset. Based on the time anchor offset, the preliminary remaining lifetime time is compensated to obtain a predicted remaining lifetime value.

[0059] Specifically, when calculating the time difference between the predicted failure time and the current time to obtain the predicted remaining service life, a difference operation is first performed between the predicted failure time and the current time. This operation is a direct subtraction, that is, subtracting the current system operating time (e.g., 8500 hours) from the predicted failure time (e.g., 10250 hours). The result is the preliminary remaining service life (here, 1750 hours). Although this value is intuitive, it does not consider the time delay effect introduced during phase space reconstruction. If it is directly used for maintenance decisions, the prediction may be premature or delayed due to reference offset.

[0060] The root of the problem lies in the fact that a time-delay embedding method is required when constructing the phase space trajectory matrix, and its core parameters include the embedding dimension. and optimal delay time For a 6-dimensional phase space reconstruction ( The actual state vector used is derived from the original time series at intervals of... The phase space consists of sampling points. This means that the "current state" represented by each trajectory point in phase space actually corresponds to the earliest moment of the original time series, rather than the moment of the latest moment of the component. Therefore, the entire degradation trajectory curve exhibits a systematic lag on the time axis.

[0061] To quantify this lag, the optimal delay time needs to be accumulated. Specifically, the reconstruction lag duration is defined as the cumulative time of the delay time at each stage during the phase space reconstruction process, and its mathematical expression is: The physical meaning of k is the order number of the delayed embedding in the phase space reconstruction, used to traverse all delay terms from the first order to the (m-1)th order to calculate the total reconstruction lag time. For example, when the embedding dimension... The optimal delay time determined by the mutual information method When the time interval is 1 hour, the reconstruction lag time is 1 hour. Hours. This value represents the inherent offset between the latest observation point of the original time series and the actual physical moment corresponding to the end of the phase space trajectory.

[0062] Subsequently, the preliminary remaining lifetime is corrected based on the reconstruction lag duration. Since the end of the degradation trajectory curve actually lags behind the current real time by a reconstruction lag duration, and the predicted failure time is extrapolated from this lag trajectory, its time base is also shifted overall. Therefore, the reconstruction lag duration needs to be considered as a time anchor offset to correct the prediction result. The specific compensation calculation method is: add this time anchor offset to the preliminary remaining lifetime duration, i.e. In the example above, 1750 hours plus 750 hours results in a predicted remaining useful life of 2500 hours.

[0063] It is important to note that this compensation is not a simple addition, but rather a restoration of the time reference frame. If this step is ignored, the predicted failure time point will be incorrectly anchored on the lag time axis of the phase space reconstruction, leading to an underestimation of the remaining lifetime. This is especially problematic in high-dimensional reconstruction or large-delay scenarios (such as when extracting early, subtle bearing fault features requires significant time adjustment). This deviation cannot be ignored.

[0064] In this embodiment, the above scheme introduces the reconstruction lag time as a time anchor offset to restore the preliminary remaining lifespan to a physical time reference, effectively eliminating the systematic time delay error introduced by phase space reconstruction. This method ensures that the predicted remaining lifespan is strictly aligned with the actual operating time of the equipment, significantly improving the engineering operability and maintenance scheduling accuracy of the prediction results in the health management of high-precision rotating machinery such as magnetic levitation flywheels.

[0065] In a specific embodiment, calculating the time difference between the predicted failure time point and the current time to obtain the predicted remaining useful life includes: Based on the optimal delay time and the embedding dimension, the current moment is shifted to the left along the time axis to obtain a phase space time window. Based on the phase space time window, the difference in the horizontal coordinate between the initial degradation trajectory curve and the smoothed degradation trajectory curve is calculated to obtain the alignment time node. The timeline segments of the alignment time node and the predicted failure time point are truncated to obtain the lifetime time segment, and the length value of the lifetime time segment is read to obtain the predicted value of the remaining lifetime.

[0066] Specifically, when calculating the time difference between the predicted failure time and the current time to obtain the predicted remaining service life, the current time must first be shifted to the left along the time axis. This shift is not arbitrary but strictly based on the optimal delay time used for phase space reconstruction. With Embedding Dimension Determined. Specifically, a complete state vector in phase space is composed of... It consists of observation points arranged in order of time delay, with the latest time being the current time. The earliest time is Therefore, the starting position of the entire phase space time window is this earliest moment, and its length is... For example, in a magnetic levitation flywheel system, if the optimal delay time is determined to be 50 hours using the mutual information method and the embedding dimension is 6, then the width of the phase space time window is 250 hours, and its starting scale is the current time (e.g., 8500 hours) minus 250 hours, i.e., 8250 hours.

[0067] Subsequently, a benchmark anchoring calculation is performed on the difference in the abscissa between the initial degraded trajectory curve and the smoothed degraded trajectory curve based on the phase space time window. Here, "difference in the abscissa" refers to the slight misalignment of the time axes of the two curves due to data preprocessing (such as filtering and interpolation) under the same physical process. Since smoothing may introduce phase lag (especially when using non-causal filters), the initial degraded trajectory curve (connecting the original sampling points) and the smoothed curve may have a shift of several hours at their ends. To eliminate this effect, the two curves need to be aligned within the phase space time window: using the start time of the phase space time window as a reference, the time offset of the end point of the smoothed curve relative to the end point of the initial curve is calculated, and this offset is used for correction. The final aligned time node is the physical time that strictly corresponds to the end of the phase space state vector after phase correction. The specific scale position of this aligned time node on the time coordinate system can be represented as follows: ,in To smooth out the introduced phase lag, it is typically obtained through cross-correlation peak detection or group delay estimation. For example, if smoothing results in a terminal lag of 12 hours, the alignment time node is 8488 hours.

[0068] After obtaining the alignment time point, it is placed on the same time axis as the predicted failure time point (e.g., 10250 hours), and a time axis segment is extracted. This operation essentially delineates a continuous interval in the time coordinate system from the alignment time point to the predicted failure time point, forming a lifetime time segment. The physical meaning of this segment is clear: it represents the complete residual evolution process from the current effective state characterization moment (corrected phase and reconfiguration lag) to system functional failure. Finally, the length of this lifetime time segment is numerically read, i.e., the difference between the two time scales is directly calculated. . It is a key intermediate variable in the prediction of remaining useful life, used to calculate the remaining useful life by subtracting it from the current time. In the numerical example above, 10250 minus 8488 yields 1762 hours as the final prediction result.

[0069] It is important to emphasize that this method avoids simply equating the "current moment" with the data acquisition moment. Instead, it restores the true physical time reference of the state representation by anchoring it to the reference through a phase space time window. If the alignment step is ignored and the original current moment is directly subtracted, tens to hundreds of hours of systematic error may be introduced due to smoothing phase distortion or reconstruction window offset.

[0070] In this embodiment, the above scheme constructs a phase-space time window based on the optimal delay time and embedding dimension, and uses the difference in the abscissa between the initial and smooth degradation trajectory curves for benchmark anchoring to accurately determine the alignment time node. Then, it extracts a lifetime time segment from this point to the predicted failure time point and reads its length, thereby obtaining a predicted remaining lifetime value synchronized with the actual dynamic state of the system. This method effectively integrates the phase-space geometric constraints in nonlinear time series analysis with the phase correction mechanism in signal processing, significantly improving the accuracy and engineering practicality of the predicted time reference in the lifetime management of high-reliability equipment such as magnetic levitation flywheels.

[0071] It should be noted that the flywheel energy storage system uses magnetic levitation bearing technology, which reduces mechanical wear, and the flywheel rotor material (such as carbon fiber composite material) has high fatigue resistance, resulting in a long overall service life.

[0072] The typical design life is around 20 years, but some advanced systems, through optimized materials and control systems, can achieve a lifespan of 25 years or even longer. The data cited above are for better understanding and are not limited to the examples given.

[0073] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention. It should be noted that the information interaction, execution process, etc. between the above devices / units are based on the same concept as the method embodiments of this application. Their specific functions and technical effects can be found in the embodiment section of the control device, and will not be repeated here.

[0074] Please see Figure 4 , Figure 4 This is a schematic diagram of the framework of an embodiment of the flywheel health status assessment and remaining life prediction device of this application. Figure 4 As shown, the flywheel health status assessment and remaining life prediction device includes a loading module 1, used to simulate the service conditions of the bolt specimen under test using a hydraulic torsion-tension composite cylinder to obtain the dynamic torsional deformation curve of the bolt; a first calculation module 2, used to calculate the torsional stiffness of the bolt specimen under test based on the dynamic torsional deformation curve of the bolt to obtain a torsional stiffness curve; an extraction module 3, used to extract the envelope extreme value based on the torsional stiffness curve to obtain the stiffness degradation, and to calculate the secant slope of the dynamic torsional deformation curve of the bolt based on the stiffness degradation to obtain the preload attenuation; and a second calculation module 4, used to calculate the fatigue damage accumulation of the bolt specimen under test based on the preload attenuation to obtain the number of loosening failure cycles.

[0075] The above module is used to perform the steps of the flywheel health status assessment and remaining life prediction method.

Claims

1. A flywheel health condition assessment and remaining life prediction method, characterized in that, Includes the following steps: The flywheel vibration signal is acquired, and time-domain statistical features are extracted from the flywheel vibration signal to obtain a multi-dimensional feature vector; Principal component analysis is performed on the multidimensional feature vectors to reduce their dimensionality, thereby obtaining health status indicators. The health status indicators are reconstructed in phase space to obtain a phase space trajectory matrix, and the phase space trajectory matrix is ​​fitted with a polynomial trend to obtain a degradation trajectory curve. The degradation trajectory curve is extrapolated based on a preset failure threshold to obtain the predicted failure time point. The time difference between the predicted failure time point and the current time is calculated to obtain the predicted remaining service life.

2. The flywheel health state assessment and remaining useful life prediction method of claim 1, wherein, The time-domain statistical feature extraction of the flywheel vibration signal yields a multi-dimensional feature vector, including: The average rectified value is obtained by taking the absolute value of each point of the flywheel vibration signal and averaging the sequence. The peak-to-peak value is obtained by subtracting the minimum value from the maximum value of the flywheel vibration signal. The flywheel vibration signal is zero-point offset corrected based on the average rectified value to obtain a centered sequence, and the square root average of the centered sequence is obtained by summing the squares based on the peak-to-peak value to obtain the root mean square value. A multidimensional feature vector is obtained by concatenating the average rectified value, the peak-to-peak value, and the root mean square value.

3. The flywheel health state assessment and remaining life prediction method of claim 1, wherein, Time-domain statistical feature extraction is performed on the flywheel vibration signal to obtain a multi-dimensional feature vector, including: The flywheel vibration signal is divided into sampling points to obtain multiple sampling point sequences; The mean of each sample point sequence is calculated to obtain the mean of each sequence, and the variance of each sample point sequence is calculated to obtain the variance of each sequence. Based on the mean and variance of each sequence, the peak value of the flywheel vibration signal is calculated to obtain the peak value of each sampling point sequence; The mean, variance, and peak values ​​of each sequence are combined to obtain a multidimensional feature vector.

4. The flywheel health state assessment and remaining life prediction method of claim 1, wherein, The step of performing principal component analysis to reduce the dimensionality of the multidimensional feature vectors to obtain health status indicators includes: The multidimensional feature vector is subjected to data standardization to obtain a standardized feature vector, and the covariance between each pair of features in the standardized feature vector is calculated to obtain a covariance matrix. The covariance matrix is ​​decomposed into eigenvalues ​​to obtain eigenvalues ​​and corresponding eigenvectors. The proportion of each eigenvalue to the sum of all eigenvalues ​​is calculated to obtain the variance contribution rate of each principal component. According to the order of variance contribution rates from largest to smallest, the eigenvectors corresponding to the first few eigenvalues ​​whose cumulative variance contribution rates reach the preset contribution rate threshold are selected to form a projection matrix. The standardized feature vector is linearly transformed using the projection matrix to obtain a low-dimensional principal component vector. The low-dimensional principal component vector is then weighted and summed based on the variance contribution rate of each principal component to obtain a health status index.

5. The flywheel health state assessment and remaining life prediction method of claim 1, wherein, The phase space reconstruction of the health status indicators to obtain the phase space trajectory matrix includes: The health status index is segmented to obtain multiple degradation stage data segments, and an autocorrelation function is calculated for each degradation stage data segment to obtain an autocorrelation function value sequence, wherein the autocorrelation function value sequence includes the autocorrelation function values ​​corresponding to each delay time. The initial delay time for each degradation stage is determined based on the delay time corresponding to the first drop of the autocorrelation function value to a preset threshold. For each data segment in the degradation stage, starting from the initial delay time, the delay time is gradually increased, and false nearest neighbor analysis is performed at each delay time. The proportion of false nearest neighbors is calculated. The proportion of false nearest neighbors is obtained by calculating the ratio of the number of false nearest neighbors in the phase space to the total number of neighbors under a given embedding dimension. Based on the delay time and embedding dimension corresponding to the first drop of the proportion of false nearest neighbors to a preset proportion threshold, the optimal delay time and embedding dimension for each degradation stage are determined. Based on the optimal delay time and embedding dimension of each degradation stage, the phase space of the data segments of each degradation stage is reconstructed to obtain the phase space trajectory sub-matrix of each degradation stage. The phase space trajectory sub-matrixes of each degradation stage are then concatenated in chronological order to obtain the phase space trajectory matrix.

6. The flywheel health state assessment and remaining life prediction method of claim 1, wherein, The step of performing polynomial trend fitting on the phase space trajectory matrix to obtain the degenerate trajectory curve includes: The phase space trajectory matrix is ​​segmented to obtain multiple trajectory segment matrices, and the mean of each trajectory segment matrix is ​​calculated to obtain the mean vector of each trajectory segment. Based on the mean vector of each trajectory segment, a piecewise polynomial fitting is performed on the phase space trajectory matrix to obtain the polynomial fitting curve of each trajectory segment. Based on the start and end time points of each trajectory segment, the polynomial fitting curves of each trajectory segment are spliced ​​together to obtain the preliminary degenerate trajectory curve after splicing. The initial degradation trajectory curve is smoothed by using the sliding window averaging method to perform smoothing calculations on each point on the initial degradation trajectory curve, resulting in a smoothed degradation trajectory curve.

7. The flywheel health state assessment and remaining useful life prediction method according to any one of claims 1-6, characterized in that, The extrapolation calculation of the degradation trajectory curve based on a preset failure threshold to obtain the predicted failure time point includes: The degraded trajectory curve is subjected to polynomial fitting parameter inversion to obtain trajectory polynomial coefficients, and the time axis extension calculation is performed on the degraded trajectory curve based on the trajectory polynomial coefficients to obtain the future predicted trajectory sequence. The predicted trajectory sequence is numerically compared with a preset failure threshold to obtain the predicted intersection time point, and the second derivative of the trajectory polynomial coefficients at the predicted intersection time point is calculated to obtain the trajectory divergence curvature. Sensitivity correction calculations are performed on the predicted intersection time points based on the trajectory divergence curvature to obtain failure time correction values, and the predicted failure time points are determined based on the failure time correction values.

8. The flywheel health state assessment and remaining useful life prediction method of claim 7, wherein, The step of performing polynomial fitting parameter inversion on the degraded trajectory curve to obtain trajectory polynomial coefficients includes: The degradation trajectory curve is discretized over time to obtain a time sampling sequence, and the time sampling sequence is calculated by exponential increment to obtain the Vandermonde matrix; The column vectors of the Vandermonde matrix are orthogonalized to obtain an orthogonal coefficient matrix. Based on the orthogonal coefficient matrix, the Vandermonde matrix is ​​transformed to obtain a condition number optimization matrix. The condition number optimization matrix is ​​solved by least squares to obtain the intermediate polynomial coefficients, and the intermediate polynomial coefficients are then subjected to higher-order term damping weighting to obtain the trajectory polynomial coefficients.

9. The method for assessing flywheel health status and predicting remaining lifespan according to claim 7, characterized in that, The step of performing time-axis extension calculations on the degraded trajectory curve based on the trajectory polynomial coefficients to obtain a future predicted trajectory sequence includes: Based on the trajectory polynomial coefficients, a polynomial evaluation calculation is performed on the preset future time series to obtain an initial extrapolation scatter plot. Then, a first-order difference calculation is performed on the adjacent data points of the initial extrapolation scatter plot to obtain an extrapolation slope curve. The extrapolation slope curve includes the future time series generated by extending backward from the end time point of the degenerate trajectory curve by a fixed step size. The extrapolated slope curve is subjected to extreme value limiting and negative value zeroing calculations to obtain a monotonic slope envelope. Based on the monotonic slope envelope, the terminal values ​​of the degenerate trajectory curve are accumulated point by point to obtain the future predicted trajectory sequence.

10. A flywheel health condition assessment and remaining life prediction system, characterized in that, The method for performing flywheel health status assessment and remaining life prediction according to any one of claims 1 to 9 includes: The extraction module is used to acquire the flywheel vibration signal, perform time-domain statistical feature extraction on the flywheel vibration signal, and obtain a multi-dimensional feature vector; The dimensionality reduction module is used to perform principal component analysis to reduce the dimensionality of the multidimensional feature vectors and obtain health status indicators. The reconstruction module is used to reconstruct the phase space of the health status indicators to obtain the phase space trajectory matrix, and to perform polynomial trend fitting on the phase space trajectory matrix to obtain the degradation trajectory curve. The calculation module is used to extrapolate the degradation trajectory curve based on a preset failure threshold to obtain the predicted failure time point, and calculate the time difference between the predicted failure time point and the current time to obtain the predicted remaining service life value.

Images