Method for predicting residual service life of industrial equipment based on time sequence decomposition

Abstract

The invention discloses a method for predicting the residual service life of industrial equipment based on time series decomposition, and the method comprises the steps: separating a seasonal component and a residual component from single-variable health indexes from complex equipment operation time series data through a time series decomposition technology, and carrying out the prediction of the residual service life of the industrial equipment; according to the method, the pure trend signal which only reflects the inherent degradation process of the equipment performance is extracted, and prediction is carried out based on the pure trend signal, so that the problems of reduced life prediction precision and poor robustness caused by the fact that the health indexes contain a large amount of non-degradation related fluctuations in a strong-periodicity or variable-working-condition operation environment in the prior art are solved.

Description

Technical Field

[0001] This invention relates to the field of remaining service life prediction technology for industrial equipment, and in particular to a method for predicting the remaining service life of industrial equipment based on time series decomposition. Background Technology

[0002] With the advancement of intelligent manufacturing, predictive maintenance of key industrial equipment and accurate prediction of its remaining service life are crucial for avoiding unplanned downtime, reducing maintenance costs, and ensuring production safety.

[0003] Current life prediction and maintenance methods typically model raw sensor data directly, failing to effectively separate equipment performance degradation trends from periodic operating condition fluctuations and random noise, resulting in low prediction accuracy and poor robustness.

[0004] The above content is only used to help understand the technical solution of the present invention and does not represent an admission that the above content is prior art. Summary of the Invention

[0005] The main objective of this invention is to provide a method for predicting the remaining service life of industrial equipment based on time series decomposition, in order to solve the problem that the current method fails to effectively separate the trend of equipment performance degradation from periodic operating condition fluctuations and random noise, resulting in low prediction accuracy and poor robustness.

[0006] To achieve the above objectives, the present invention provides a method for predicting the remaining useful life of industrial equipment based on time series decomposition, the method comprising: A multi-source dataset is acquired from multiple sensors of industrial equipment, and the data in the multi-source dataset is preprocessed to obtain a multivariate time series dataset. From the multivariate time series dataset, extract or construct a univariate health indicator that can comprehensively reflect the overall health status of industrial equipment; The univariate health index is decomposed into three components using a time series decomposition algorithm: a trend component, a seasonal component, and a residual component. The seasonal component and the residual component are then separated from the univariate health index. The trend component is evaluated for quality to obtain a quality score corresponding to the trend component. Based on the quality score, physical constraints are adaptively applied to the trend components to optimize the trend components and obtain the optimized trend component sequence. Input the optimized trend component sequence of a preset length before the current moment into the prediction model to predict the trend value in the future period. A failure threshold is determined, the predicted trend value is compared with the failure threshold, and the remaining life of the industrial equipment is calculated based on the comparison result.

[0007] Furthermore, the step of constructing a univariate health index that can comprehensively reflect the overall health status of industrial equipment includes: The multivariate data in the multivariate time series dataset is standardized to obtain standardized multivariate data. Calculate the covariance matrix for the standardized multivariate data; The covariance matrix is ​​decomposed according to the eigenvalue decomposition formula, and eigenvalues ​​and eigenvectors are extracted. Take the largest eigenvalue and the largest eigenvector of the largest eigenvalue, and calculate the univariate health index based on the largest eigenvalue and the largest eigenvector.

[0008] Furthermore, the time series decomposition algorithm is the STL (Seasonal and Trend decomposition using Loess) algorithm, and the step of decomposing the univariate health indicator into three components using the time series decomposition algorithm includes: Determine the disassembly cycle of the industrial equipment; The trend component and the seasonal component are calculated through an inner loop iteration; Calculate the residual components and update the robustness weights based on the trend components and seasonality components; Based on the robustness weights, determine whether to continue the inner loop iteration. If the determination result is not to continue the inner loop iteration, then output the decomposed trend component, the seasonal component, and the residual component.

[0009] Furthermore, the step of calculating the trend component and the seasonal component through inner loop iteration includes: Obtain the initial trend component and calculate the detrended sequence; The detrended sequence is decomposed into multiple sub-sequences according to the decomposition period; The multiple subsequences are smoothed by local weighted regression smoothing (Loess) to obtain multiple smoothed subsequences; The multiple subsequences are subjected to two moving averages, and the moving average results after the two moving averages are subjected to a second local weighted regression smoothing to obtain the seasonal component; The seasonal component is subtracted from the detrended sequence to obtain the deseasonal sequence, and the trend component is obtained by smoothing the deseasonal sequence using local weighted regression.

[0010] Further, the step of performing a quality assessment on the trend component to obtain a quality score corresponding to the trend component includes: Calculate the monotonicity index, signal-to-noise ratio, and smoothness of the trend component; A quality score is calculated based on the monotonicity index, the signal-to-noise ratio, and the smoothness.

[0011] Further, the physical constraint optimization process includes at least one of monotonicity constraints, degradation rate boundary constraints, and smoothness constraints. The step of adaptively applying physical constraint optimization processing to the trend component based on the quality score to obtain the optimized trend component sequence includes: The weight parameters of the monotonicity constraint, the degradation rate boundary constraint, and the smoothness constraint are calculated based on the quality score. The lower the quality score, the stronger the constraint is. The physical constraint optimization process is achieved by solving the following optimization problem: ,in, Weights constrained by monotonicity. For the rate boundary constraint weights, To smooth the constraint weights, For monotonicity penalty terms, For the rate boundary penalty term, To smooth out the penalty term, To rate the quality, For the original trend component, The optimized trend component, For a certain moment.

[0012] Furthermore, based on the quality score, the monotonicity constraint weight, the rate boundary constraint weight, and the smoothness constraint weight are adaptively adjusted according to the following rules: in, For maximum monotonicity weight, The maximum rate boundary weights, For maximum smoothing weights, , , These are the attenuation coefficients, Rate the quality.

[0013] Furthermore, the step of inputting the optimized trend component sequence of a preset length prior to the current time into the prediction model to predict the predicted trend value over a future period includes: Get the current time value and the length of the history window; Multiple data points are extracted from the optimized trend component sequence based on the time value and the historical window length to form an input sequence; The input sequence is normalized and then input into the prediction model.

[0014] Furthermore, the step of determining the failure threshold includes: Acquire historical data of industrial equipment throughout its entire lifecycle, from start of operation to failure, and extract the failure trend component at the moment of failure from the historical data; The failure threshold is calculated using an arithmetic average algorithm based on the failure trend components.

[0015] Further, the step of comparing the predicted trend value with the failure threshold and calculating the remaining life of the industrial equipment based on the comparison result includes: Find the first crossing point, wherein the first crossing point is the time point at which the first of the plurality of predicted trend values ​​reaches or exceeds the failure threshold; The remaining lifespan of industrial equipment is calculated using a preset formula, wherein the preset formula is: ,in, This is the predicted trend value of the step preceding the k-th step (i.e., the k-1-th step).

[0016] This application uses time series decomposition technology to separate seasonal and residual components from univariate health indicators in complex equipment operation sequence data, extracting a pure trend signal that only reflects the inherent degradation process of equipment performance, and making predictions based on this pure trend signal. This solves the problem in the prior art that the accuracy and robustness of life prediction are reduced due to the large number of non-degradation-related fluctuations in health indicators under strong periodic or variable operating conditions. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating an embodiment of the method for predicting the remaining useful life of industrial equipment based on time series decomposition according to this application. The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0019] This invention further provides a method for predicting the remaining useful life of industrial equipment based on time series decomposition. (Refer to...) Figure 1 , Figure 1 This is a flowchart illustrating an embodiment of the industrial equipment remaining service life prediction method based on time series decomposition according to the present invention.

[0020] In this embodiment, the execution entity of the industrial equipment remaining useful life prediction method based on time series decomposition is an industrial equipment remaining useful life prediction method system based on time series decomposition. This system includes an industrial equipment remaining useful life prediction method device or apparatus based on time series decomposition, which can be a PC, PDA, or other terminal device. This application uses time series decomposition technology to separate seasonal and residual components from univariate health indicators in complex equipment operating sequence data, extracting a pure trend signal that only reflects the inherent degradation process of equipment performance. Prediction is then performed based on this pure trend signal, solving the problem in existing technologies where, under strongly periodic or variable operating conditions, the health indicators contain a large number of non-degradation-related fluctuations, leading to decreased accuracy and poor robustness in life prediction.

[0021] The steps of this method for predicting the remaining useful life of industrial equipment based on time series decomposition include: Step S10: Obtain multi-source datasets collected by multiple sensors of industrial equipment, and perform data preprocessing on the data in the multi-source datasets to obtain multivariate time series datasets; In this embodiment, the multi-source dataset refers to the raw data collection from sensors of different types, locations, and physical quantities installed on the target industrial equipment. Examples include vibration sensors (accelerometers, velocity sensors, displacement sensors) for collecting time-domain / frequency-domain vibration signals; temperature sensors (thermocouples, RTDs) for monitoring the temperature of key parts of the equipment; pressure sensors for monitoring hydraulic and pneumatic system pressure; current / voltage sensors for monitoring the electrical parameters of motors and drives; acoustic sensors for collecting equipment operating noise; speed sensors for monitoring the rotational speed of rotating machinery; and position / displacement sensors for monitoring the relative position changes of mechanical components. Data preprocessing includes data cleaning, time alignment, filtering and denoising, and normalization of the collected multi-source dataset. Data cleaning handles missing values, including interpolation filling or marker deletion. Time alignment unifies sensor data with different sampling frequencies to the same timestamp through interpolation. Filtering and denoising removes noise in specific frequency bands using digital filtering (such as bandpass filtering). Normalization scales data with different dimensions and ranges to a uniform scale (e.g., the [0,1] interval).

[0022] A multivariate time series dataset refers to a collection of data containing observations of multiple variables (sensors) arranged in chronological order after data preprocessing. For example, an industrial device has M sensors, and the observations at time t are: The preprocessed multivariate time series dataset is then: in, This represents the preprocessed observation values ​​of the Mth sensor at N time points.

[0023] Step S20: Extract or construct a univariate health indicator from the multivariate time series dataset that can comprehensively reflect the overall health status of the equipment. In this embodiment, the Univariate Health Indicator (HI) is a single numerical sequence that comprehensively reflects the overall health status of a device over time. It compresses complex, multi-dimensional sensor information into a scalar that varies over time, characterizing the degradation process of a device from health to failure.

[0024] Extracting a single-variable health indicator that comprehensively reflects the overall health status of the equipment involves selecting a physical quantity from multiple sensors that is most directly related to equipment degradation. For example, for bearings, the root mean square value of vibration acceleration would be chosen. .in, The health index at time t represents the health status of the equipment at time t. The larger the value, the more severe the degradation. For the sliding window time; For vibration sensors in time Observed values; Here is the index within the window, where i is an integer. .

[0025] Principal component analysis can be used to construct a univariate health index that can comprehensively reflect the overall health status of equipment. First, the multivariate data in the multivariate time series dataset is standardized to obtain standardized multivariate data. Then, the covariance matrix is ​​calculated on the standardized multivariate data. The covariance formula is then decomposed to extract eigenvalues ​​and eigenvectors. The largest eigenvalue and the largest eigenvector of the largest eigenvalue are selected. Finally, the univariate health index is constructed.

[0026] The standardized formula is ,in, For the standardized m-th sensor in the first position The values ​​at each time point can be standardized using a formula to eliminate differences in the dimensions and numerical ranges of different sensors. Let be the mean of all observations of the m-th sensor during the initial health phase of the device; Let be the standard deviation of all observations of the m-th sensor during the initial health phase of the equipment. After the equipment is installed and debugged and has been running normally for a period of time, the sensor observations are collected during this period. These observations serve as a reference for each sensor in its normal state, thus providing a sufficient sample size to reliably estimate the statistical characteristics.

[0027] The formula for calculating the covariance matrix is: N is the total number of time points. , which is the standardized observation vector at a certain point in time.

[0028] The eigenvalue decomposition formula is , For the i-th eigenvalue, For the corresponding eigenvectors, The index of the eigenvalue or eigenvector.

[0029] The formula for constructing a univariate health index is: , The eigenvector corresponding to the largest eigenvalue. The j-th component of the largest eigenvector corresponding to the largest eigenvalue. This is the standardized value of the j-th sensor at the k-th time point.

[0030] Step S30: Decompose the univariate health index into three components through time series decomposition: trend component, seasonal component, and residual component, and separate the seasonal component and residual component from the univariate health index. In this embodiment, the trend component reflects the long-term, monotonous, and irreversible degradation process of equipment performance, the seasonal component reflects repetitive fluctuations caused by periodic work cycles, environmental changes, etc., and the residual component reflects random noise, measurement errors, and unpredictable transient disturbances. This embodiment utilizes only the trend component to separate the seasonal and residual components from the univariate health index. , As a univariate health indicator, As a trend component, As a seasonal component, For the residual components.

[0031] Step S30 includes: Step S31: Determine the disassembly cycle of the industrial equipment; In this embodiment, the decomposition period p is set according to the equipment operating characteristics (determined according to the equipment's working cycle), for example, p=24 (hour data) or p=24×60 (minute data), p=168 (hour data), p==720 (hour data).

[0032] Step S32: Calculate the trend component and the seasonal component through an inner loop iteration; In this embodiment, the initial trend component is first obtained, and the detrended sequence is calculated. Let the original sequence be... The period is p, and the initial sequence is... The seasonal component is Calculate the detrended sequence. ,in, This is the original univariate health indicator sequence. This is the trend component estimate for the (j-1)th iteration, where N is the total sequence length. In the first iteration... .

[0033] The detrending sequence is decomposed into multiple subsequences according to the decomposition period. For example, a detrending sequence of length N is folded into p subsequences according to the decomposition period p, and each subsequence corresponds to a position in the period: .

[0034] Multiple subsequences are smoothed using local weighted regression to obtain smoothed subsequences, with a window length of [value missing]. polynomial order , The smoothed subsequences are then rearranged in their original chronological order. .in, The period length, For seasonal smoothing window, For the subsequence at the k-th period position, The value of the k-th subsequence after smoothing. and Same length, The result is a smoothed version of the recombined periodic subsequences. For robust weight vectors, belong .

[0035] Perform two moving averages on multiple subsequences. The first moving average is: ,in, This is the result of the first moving average. Symmetrical expansion or truncation is used. First moving average: ,in, This is the result of the second average shift.

[0036] A second local weighted regression smoothing is performed on the moving average results after two moving averages to obtain the seasonal component. The second moving average is as follows: ,in, This is the output result after the second moving average. For window length, usually The pure seasonal component is separated from the output after the second moving average using the following formula: To perform centralized processing, the average value at each period position is first calculated using the following formula: Subtract the average value corresponding to the period position from each value. And verify the property, the verification formula is: .in, For decentralized seasonal components, The average value of the k-th period position. For centralized seasonal components, This represents the number of data points at the k-th cycle position.

[0037] Subtracting the seasonal component from the detrended sequence yields the deseasonalized sequence. The purpose of this deseasonalized sequence is to remove the seasonal component from the original sequence, preparing for trend estimation. This is represented as... , This is a deseasonalized sequence.

[0038] The trend component is obtained by smoothing the deseasonalized series using locally weighted regression. The purpose is to estimate the trend component from the deseasonalized series. The formula is as follows: ,in, For trend smoothing windows, usually =1.5p (adjusted to an odd number), It is a linear local polynomial. For robustness weights.

[0039] For each point The solution is obtained by using local weighted regression to achieve a smooth fit, and the formula is as follows: ,in, It is a cubic weight function. , It is a positive real number. For bandwidth, For robustness weights , These are local linear coefficients. Final trend component. .

[0040] Step S33: Calculate the residual component and update the robustness weight based on the trend component and seasonal component; In this embodiment, the formula for calculating the residual component is as follows: ,in, Let be the residual sequence of the k-th outer loop. The original univariate health indicator sequence, For the estimation of the trend component after the k-th outer loop, This is the estimate of the seasonal component after the kth external cycle.

[0041] The purpose of calculating the median absolute deviation of the residuals is to obtain a robust estimate of the dispersion to outliers. The calculation formula is as follows: Specifically, first calculate the median of the residuals. , Then calculate the absolute value deviation. , Finally, the median of the absolute deviation is calculated. , .in, The MAD value of the kth outer loop iteration. The standard residual is calculated using the following formula: ,in, For standard residuals, This is an adjustment constant, typically 6, used to control the strictness of outlier detection. When When <1, the data point is considered normal. When the value is ≥1, the data point is considered an outlier.

[0042] Updating robustness weights aims to assign weights to each data point and reduce the impact of outliers. , The weight of the (k+1)th outer loop iteration. It is a double square function. This represents the standard residual. The double square function is: For normal data points (u<1), the weight is close to 1, and they fully participate in the next round of calculation; for slightly abnormal data points (u≈0.5−0.9), the weight is reduced, and the influence is weakened; for significantly abnormal data points (u≥1), the weight is 0, and they are completely excluded.

[0043] Step S34: Determine whether to continue the inner loop iteration based on the robustness weight. If the determination result is not to continue the inner loop iteration, output the decomposed trend component, seasonal component and residual component. In this embodiment, a fixed number of loop iterations can be set to determine the number of inner loop iterations and initiate the inner loop iteration termination condition. Alternatively, an adaptive termination condition can be used to determine whether to continue the inner loop iteration. In the adaptive termination condition, the rate of change of the weights is first calculated using the following formula: The termination condition is: ,in, For tolerance (e.g., 0.01), This represents the maximum number of iterations, for example, 3.

[0044] When the loop terminates, take the result of the last iteration: Step S40: Perform a quality assessment on the trend component to obtain the quality score corresponding to the trend component. In this embodiment, the monotonicity index, signal-to-noise ratio, and smoothness of the trend component are first calculated, and then the quality score is calculated using the monotonicity index, signal-to-noise ratio, and smoothness.

[0045] Step S40 above includes: Step S41: Calculate the monotonicity exponent, signal-to-noise ratio, and smoothness of the trend component; In this embodiment, the formula for calculating the monotonicity index is: , ,in, This is an indicator function that takes the value 1 when the condition is true and 0 when the condition is false. It is a monotonicity index. The length of the trend component sequence. Let i be the trend component value at the i-th time point. This represents the difference between adjacent trend components. The tolerance threshold is usually taken as... .

[0046] In the signal-to-noise ratio calculation process, the variance of the trend component is first calculated using the following formula: , Next, calculate the variance of the seasonal component using the following formula: ; Calculate the variance of the residual components using the following formula: The signal-to-noise ratio of the trend component is calculated using the following formula: Finally, the signal-to-noise ratio can be normalized. The normalization formula is as follows: .in, Let Variance be the variance of the potential component. The variance of the functional components, Let Variance be the variance of the residual components. The signal-to-noise ratio of the trend component. This represents the maximum empirical value for the signal-to-noise ratio. The minimum empirical value for signal-to-noise ratio.

[0047] The smoothness calculation process requires calculating the second-order difference and the roughness index, and then calculating the smoothness index. The specific calculation formula is as follows: Where α is the adjustment parameter, usually , The second difference of the trend component. Roughness index This is the smoothness index.

[0048] Step S42: Calculate the quality score based on the monotonicity index, signal-to-noise ratio, and smoothness.

[0049] In this embodiment, if necessary, each indicator can be normalized first, and then a weighted sum can be performed, as follows: The weighting index satisfies: .For example, It is 0.4. It is 0.4. It is 0.2.

[0050] Step S50: Based on the quality score, the trend components are adaptively subjected to physical constraints for optimization to obtain the optimized trend component sequence; In this embodiment, the physical constraint optimization process includes at least one of monotonicity constraints, degradation rate boundary constraints, and smoothness constraints. The monotonicity constraint ensures that the optimized trend component sequence is non-decreasing, the degradation rate boundary constraint limits the degradation rate between adjacent time points to a preset range, and the smoothness constraint controls the smoothness of the optimized trend component sequence. Weight parameters for the monotonicity constraints, degradation rate boundary constraints, and smoothness constraints are calculated based on the quality score, where a lower quality score results in a stronger constraint.

[0051] Assuming the original trend components are The optimized trend component is The optimized trend component sequence is obtained by optimizing each trend component. The optimization objective function is: in, Weights constrained by monotonicity. For the rate boundary constraint weights, To smooth the constraint weights, For monotonicity penalty terms, For the rate boundary penalty term, To smooth out the penalty term, To rate the quality, Greater than or equal to 0, less than or equal to 1.

[0052] Based on the quality score Q, the constraint weights are adaptively adjusted according to the following rules: in, The maximum monotonicity weight is the weight when Q=0 (worst quality), and Q can be 10; Q can be 5, which is the maximum rate boundary weight. For maximum smoothing weights, Q can be 2; , , These are the attenuation coefficients, both of which can be 2.

[0053] Monotonicity penalty item for: Parameter adaptation: in, To minimize the degradation increment, For the minimum allowable increment, For the maximum allowed increment, For exponential parameters, The value is between 0.001 and 0.005. It is 0.02-0.05. It is 1.5-3.0. When < When Q decreases, a penalty is imposed; the amount of the penalty is proportional to the square of the decrease; the smaller Q is, the larger ε is, and the stricter the constraint (allowing for a larger minimum increment).

[0054] Rate boundary penalty term for: Rate boundary adaptation: in, The sampling time interval, such as 1 hour or 1 day; , As the baseline degradation rate boundary, To extend the range to the boundary, This is the exponential parameter, which can range from 1.5 to 3.0. For example, =0.001, When =0.01 and Q=0.3: =−0.002, =0.013 (relaxed boundary); Q=0.8: =0.0005, =0.0105 (tightening the boundary).

[0055] Smoothing penalty for: Smoothing parameter adaptation: in, For smoothness, For the minimum smoothing parameter, It can be 0.1-0.3. For the maximum smoothing parameter, It can be 1.5-2.0. For exponential parameters, It can be 1.5-3.0. This represents the original trend standard deviation.

[0056] The above optimization problem is solved using the iterative reweighted least squares method. After 5 iterations, the solution converges, yielding the optimized trend component sequence.

[0057] Step S60: Based on the quality score, the trend components are adaptively subjected to physical constraints for optimization to obtain the optimized trend component sequence; In this embodiment, the current time value t is first obtained. c And the historical window length L (preset parameter), extract from the optimized trend component sequence (t) c -L+1) to t c A total of L data points constitute the input sequence. The input sequence is normalized and then input into the prediction model. The specific steps are as follows: Input sequence construction: ,in The input sequence is .

[0058] Normalization process: in, The mean of the training set. The standard deviation of the training set. For the data points in the training set, This is the input sequence after normalization.

[0059] Will After being input into the prediction model, normalized predicted values ​​are obtained. : The normalized predicted values ​​are then denormalized to obtain the denormalized predicted values. The inverse normalized predicted value Forecast trend values ​​for a future period of time: ,Right now Step S70: Determine the failure threshold, compare the predicted trend value with the failure threshold, and calculate the remaining life of the industrial equipment based on the comparison result.

[0060] In this embodiment, the failure threshold can be directly specified as an absolute value based on the equipment's design safety margin, industry standards, or the experience of domain experts, or the failure threshold of industrial equipment can be calculated.

[0061] The failure threshold is calculated as follows: First, historical data of the industrial equipment throughout its entire lifecycle, from start of operation to failure, is acquired. The failure trend component at the moment of failure is then extracted from the historical data. Based on this failure trend component, the failure threshold is calculated using an arithmetic average algorithm. The calculation of the failure trend component can be referenced from steps S10-S60 and will not be elaborated upon here.

[0062] ,in, This is the failure threshold.

[0063] Find the first crossing point, where the first crossing point is the time point at which the first of the plurality of predicted trend values ​​reaches or exceeds the failure threshold. This is done by iterating through all the predicted trend values ​​(predicted sequences) in step S60. Search for conditions The smallest positive integer k. This value of k represents the value from the current time t. c Index of the number of prediction steps from the start of the prediction until the predicted failure occurs.

[0064] The remaining life of industrial equipment is calculated using a preset formula. Considering that the predicted value is an estimate at discrete time points, while the actual degradation crossover may occur between two discrete time points, linear interpolation is used to improve the accuracy of the remaining life calculation. The preset formula is as follows: in, The remaining lifespan of the industrial equipment is given by k, where k is the first time the failure threshold is reached or exceeded in the prediction sequence. The step index (j=k). This is the predicted trend value at step k (i.e., the first crossing point). This is the predicted trend value of the step preceding the k-th step (i.e., the (k-1)-th step). When k=1, this is the trend value at the current moment. , The temporal resolution of the prediction model is the actual physical time interval corresponding to each prediction step j. This parameter is consistent with the data acquisition cycle and the historical window step size for constructing the input sequence in step S50 (for example, if data is collected once per hour, and each time point represents 1 hour, then...). =1 hour).

[0065] Boundary case handling: Case 1, if for all j = 1, 2, ..., H, we have < Then it is determined that the equipment is within the predicted time range. The internal condition will not fail. In this case, RUL>H×Δt can be output; Case 2, if the current health status has deteriorated significantly, that is... ≥ If the device is in or has exceeded the failure state, it is determined that the device is in or has exceeded the failure state, and the output is RUL=0.

[0066] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or system that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or system. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or system that includes that element.

[0067] The sequence numbers of the above embodiments of the present invention are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0068] Through the above description of the embodiments, those skilled in the art can clearly understand that the methods of the above embodiments can be implemented by means of software plus necessary general-purpose hardware platforms. Of course, they can also be implemented by hardware, but in many cases the former is a better implementation method. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk) as described above, and includes several instructions to cause a terminal device (which may be a mobile phone, computer, server, air conditioner, or network device, etc.) to execute the methods described in the various embodiments of the present invention.

[0069] The above are merely preferred embodiments of the present invention and do not limit the scope of the patent. Any equivalent structural or procedural transformations made based on the description and drawings of the present invention, or direct or indirect applications in other related technical fields, are similarly included within the scope of patent protection of the present invention.

Claims

1. A method for predicting the remaining useful life of industrial equipment based on time series decomposition, characterized in that, The method includes: A multi-source dataset is acquired from multiple sensors of industrial equipment, and the data in the multi-source dataset is preprocessed to obtain a multivariate time series dataset. From the multivariate time series dataset, extract or construct a univariate health indicator that can comprehensively reflect the overall health status of industrial equipment; The univariate health index is decomposed into three components using a time series decomposition algorithm: a trend component, a seasonal component, and a residual component. The seasonal component and the residual component are then separated from the univariate health index. The trend component is evaluated for quality to obtain the quality score corresponding to the trend component; Based on the quality score, physical constraints are adaptively applied to the trend components to optimize the trend components and obtain the optimized trend component sequence. Input the optimized trend component sequence of a preset length before the current moment into the prediction model to predict the trend value in the future period. A failure threshold is determined, the predicted trend value is compared with the failure threshold, and the remaining life of the industrial equipment is calculated based on the comparison result.

2. The method as described in claim 1, characterized in that, The steps for constructing a univariate health index that can comprehensively reflect the overall health status of industrial equipment include: The multivariate data in the multivariate time series dataset is standardized to obtain standardized multivariate data. Calculate the covariance matrix for the standardized multivariate data; The covariance matrix is ​​decomposed according to the eigenvalue decomposition formula, and eigenvalues ​​and eigenvectors are extracted. Take the largest eigenvalue and the largest eigenvector of the largest eigenvalue, and calculate the univariate health index based on the largest eigenvalue and the largest eigenvector.

3. The method as described in claim 1, characterized in that, The time series decomposition algorithm is the STL (Seasonal and Trend decomposition using Loess) algorithm. The steps of decomposing the univariate health indicator into three components using the time series decomposition algorithm include: Determine the disassembly cycle of the industrial equipment; The trend component and the seasonal component are calculated through an inner loop iteration; Calculate the residual components and update the robustness weights based on the trend components and seasonality components; Based on the robustness weights, determine whether to continue the inner loop iteration. If the determination result is not to continue the inner loop iteration, then output the decomposed trend component, the seasonal component, and the residual component.

4. The method as described in claim 3, characterized in that, The step of calculating the trend component and the seasonal component through inner loop iteration includes: Obtain the initial trend component and calculate the detrended sequence; The detrended sequence is decomposed into multiple sub-sequences according to the decomposition period; The multiple subsequences are smoothed by local weighted regression smoothing (Loess) to obtain multiple smoothed subsequences; The multiple subsequences are subjected to two moving averages, and the moving average results after the two moving averages are subjected to a second local weighted regression smoothing to obtain the seasonal component; The seasonal component is subtracted from the detrended sequence to obtain the deseasonal sequence, and the trend component is obtained by smoothing the deseasonal sequence using local weighted regression.

5. The method as described in claim 1, characterized in that, The step of performing a quality assessment on the trend component to obtain a quality score corresponding to the trend component includes: Calculate the monotonicity index, signal-to-noise ratio, and smoothness of the trend component; A quality score is calculated based on the monotonicity index, the signal-to-noise ratio, and the smoothness.

6. The method as described in claim 1, characterized in that, The physical constraint optimization process includes at least one of monotonicity constraints, degradation rate boundary constraints, and smoothness constraints. The step of adaptively applying physical constraint optimization to the trend components based on the quality score to obtain the optimized trend component sequence includes: The weight parameters of the monotonicity constraint, the degradation rate boundary constraint, and the smoothness constraint are calculated based on the quality score. The lower the quality score, the stronger the constraint is. The physical constraint optimization process is achieved by solving the following optimization problem: ,in, Weights constrained by monotonicity. For the rate boundary constraint weights, To smooth the constraint weights, For monotonicity penalty terms, For the rate boundary penalty term, To smooth out the penalty term, To rate the quality, For the original trend component, The optimized trend component, For a certain moment.

7. The method as described in claim 6, characterized in that, Based on the quality score, the monotonicity constraint weight, the rate boundary constraint weight, and the smoothness constraint weight are adaptively adjusted according to the following rules: in, For maximum monotonicity weight, The maximum rate boundary weights, For maximum smoothing weights, , , These are the attenuation coefficients, Rate the quality.

8. The method as described in claim 1, characterized in that, The step of inputting the optimized trend component sequence of a preset length prior to the current moment into the prediction model to predict the predicted trend value over a future period includes: Get the current time value and the length of the history window; Multiple data points are extracted from the optimized trend component sequence based on the time value and the historical window length to form an input sequence; The input sequence is normalized and then input into the prediction model.

9. The method as described in claim 1, characterized in that, The step of determining the failure threshold includes: Acquire historical data of industrial equipment throughout its entire lifecycle, from start of operation to failure, and extract the failure trend component at the moment of failure from the historical data; The failure threshold is calculated using an arithmetic average algorithm based on the failure trend components.

10. The method as described in claim 9, characterized in that, The step of comparing the predicted trend value with the failure threshold and calculating the remaining life of the industrial equipment based on the comparison result includes: Find the first crossing point, wherein the first crossing point is the time point at which the first of the plurality of predicted trend values ​​reaches or exceeds the failure threshold; The remaining lifespan of industrial equipment is calculated using a preset formula, wherein the preset formula is: ,in, This is the predicted trend value of the step preceding the k-th step (i.e., the k-1-th step).

Images