Industrial robot condition monitoring and predictive maintenance method

Abstract

This invention relates to the field of industrial robot maintenance technology, specifically to a method for industrial robot condition monitoring and predictive maintenance. The method includes: acquiring and fusing multidimensional operational and environmental data; calculating the Mahalanobis distance to the healthy centroid using kernel principal component analysis projection as a first state covariate; constructing a Weibull baseline risk function by classifying historical task types; calculating the current kurtosis and summing it with a forgetting factor to obtain a second time-varying covariate; constructing a proportional hazards model; combining the current task baseline function and the two types of covariates to output the remaining lifetime probability distribution and expectation; defining a utility function that includes spare parts, downtime losses, and expected failure risk costs; optimizing the maintenance window to determine the maintenance timing that maximizes utility and generating instructions. In other words, the solution of this invention can improve the accuracy of condition assessment and lifetime prediction, achieve optimal allocation of maintenance resources, and ensure production continuity.

Description

Technical Field

[0001] This invention relates to the field of industrial robot maintenance technology. More specifically, this invention relates to methods for industrial robot condition monitoring and predictive maintenance. Background Technology

[0002] Industrial robots and related logistics facilities are core equipment in intelligent manufacturing production lines, and their reliability and stability are crucial for ensuring production efficiency and reducing operating costs. In modern production scenarios, the application scope of industrial robots is no longer limited to processing equipment on the main production line, but also encompasses a large number of logistics and auxiliary equipment serving the production line, such as robots responsible for feeding materials, packaging robots, marking robots, conveyor lines that seamlessly connect with production line processes, and various pre-processing lines. For the above-mentioned equipment, existing maintenance strategies mainly rely on periodic maintenance or post-failure maintenance. However, periodic maintenance often fails to match the actual health condition of the robot, and this blindness often leads to under-maintenance or over-maintenance; while post-failure maintenance means addressing the problem only after a failure occurs, which inevitably causes unplanned downtime, leading to serious production interruptions and economic losses. Therefore, to overcome the limitations of the above-mentioned traditional maintenance strategies, the industry urgently needs to develop advanced condition monitoring and predictive maintenance technologies, which aim to proactively intervene before failures occur by monitoring the robot's health status in real time and accurately predicting its remaining service life.

[0003] To achieve these goals, data-driven predictive maintenance methods have become the mainstream technological direction. This approach typically collects multi-source data such as vibration, current, and temperature from sensors installed on the robot itself or in the working environment. Signal processing techniques are used to extract features representing the degradation state, and machine learning or deep learning models are employed to construct health indicators, thereby predicting the remaining service life. This type of improved technology attempts to replace traditional experience-based judgment with data mining and model building, aiming to obtain more accurate maintenance data.

[0004] However, the aforementioned data-driven predictive maintenance methods still have significant shortcomings in practical applications. First, due to the highly nonlinear nature of the degradation process of industrial robots, many existing health status assessment methods struggle to accurately detect changes in the state space, resulting in insufficient sensitivity of the constructed health indicators to early, minor faults. Second, most remaining service life prediction models fail to adequately consider the variability of robot tasks, ignoring the significant differences in wear rates of robot components caused by different types of tasks. Furthermore, the robot's health status is the result of both current operational deviations and historical accumulated damage, but many models focus only on one aspect, lacking a comprehensive analytical framework that combines instantaneous status with historical degradation trends. In addition, regarding maintenance decisions, most methods simply set a threshold based on the predicted remaining service life to trigger maintenance. This strategy is too coarse, failing to comprehensively weigh the multiple costs of spare parts replacement, downtime losses, and potential failure risks incurred due to delayed maintenance, making it difficult to achieve optimal allocation of maintenance resources. Summary of the Invention

[0005] The purpose of this invention is to propose a method for monitoring the condition and predictive maintenance of industrial robots, in order to solve the problem of the difficulty in achieving optimal allocation of maintenance resources in the prior art; to this end, this invention provides a solution in one aspect.

[0006] The industrial robot condition monitoring and predictive maintenance method provided by this invention includes the following steps: Multidimensional operating parameters and environmental sensor data of the industrial robot are acquired and fused into a high-dimensional feature vector through time-domain alignment. Kernel principal component analysis is performed on the high-dimensional feature vector to project the vector of the current state onto a nonlinear principal component space. The Mahalanobis distance between the projection point and the centroid of a preset healthy state cluster is calculated, and the distance value is used as the first state covariate. The robot's historical task sequences are identified and task types are classified. A baseline risk function based on the Weibull distribution is constructed for each task type based on historical data. The kurtosis index of the robot joint drive current is calculated, and the kurtosis index is cumulatively summed with a preset forgetting factor within a sliding time window to obtain a second time-varying covariate. A proportional hazards model is constructed, with the first state covariate and the second time-varying covariate as inputs. The corresponding baseline risk function is selected according to the task type currently being performed by the robot. The probability density function of the robot's remaining service life is output using the proportional hazards model, and the expected value of the remaining service life is calculated. A utility function is defined, including spare parts cost, downtime loss cost, and expected failure risk cost calculated based on the probability density function of the remaining service life. By optimizing within a preset discrete maintenance window, the maintenance timing that maximizes the utility function is determined, and maintenance instructions are generated.

[0007] Preferably, the acquisition of multidimensional operating parameters and working environment sensor data of the industrial robot, and the fusion of these data into a high-dimensional feature vector through time-domain alignment, includes: acquiring angle, angular velocity, and torque data of the robot's six-axis joints at a sampling frequency of 100Hz, as well as three-dimensional acceleration data of the end effector; acquiring temperature and humidity data of the working environment at a sampling frequency of 1Hz; upsampling the temperature and humidity data to 100Hz using linear interpolation, aligning it with the robot's operating parameters on the timestamp, and splicing them together to form a high-dimensional feature vector with a dimension of 23.

[0008] Preferably, the step of performing kernel principal component analysis on the high-dimensional feature vector, projecting the vector of the current state onto the nonlinear principal component space, and calculating the Mahalanobis distance between the projection point and the centroid of the preset healthy state cluster includes: using a Gaussian radial basis function as the kernel function of the kernel principal component analysis; selecting several high-dimensional feature vectors collected in the initial stage of robot operation as a healthy state sample set, calculating the geometric center of the sample set after being projected onto the first K principal component dimensions by kernel principal component analysis as the centroid; and after projecting the high-dimensional feature vector of the current state, calculating the Mahalanobis distance between the projection point and the centroid in the K-dimensional principal component space.

[0009] Preferably, the step of selecting several high-dimensional feature vectors collected during the initial stage of robot operation as a health status sample set includes: selecting 5000 high-dimensional feature vectors collected within the first 100 hours of robot operation as the health status sample set; the first K principal component dimensions are the first 3 principal component dimensions.

[0010] Preferably, the process of identifying the robot's historical task sequence and classifying task types includes: dividing historical data into independent task segments based on the robot's unique task program identifier and start timestamp; extracting the average joint angular velocity and average load torque of each task segment as features; and using the K-means clustering algorithm to cluster the features, classifying the historical tasks into three types: heavy-load handling, light-load assembly, and unloaded positioning.

[0011] Preferably, the step of constructing a baseline risk function based on the Weibull distribution for each task type based on historical data includes: for each task type, collecting historical failure time data of multiple robots of the same model under that task type; and fitting the shape parameters of the Weibull distribution using the maximum likelihood estimation method. and scale parameters According to the formula Construct the baseline risk function, where, This is the baseline risk function.

[0012] Preferably, the calculation of the kurtosis index of the robot joint drive current, and the cumulative summation of the kurtosis index with a preset forgetting factor within a sliding time window to obtain a second time-varying covariate, includes: setting the width and sliding step size of the sliding time window; and calculating the kurtosis index of the drive current signal within each time window. The second time-varying covariate is calculated using an exponentially weighted cumulative summation method, and the calculation formula is as follows: ,in, Let be the second time-varying covariate at the current time t. For window width, The preset forgetting factor, and .

[0013] Preferably, the width of the sliding time window is 1000 data sampling points, the sliding step size is 100 data sampling points, and the forgetting factor... It is 0.98.

[0014] Preferably, the definition includes a utility function comprising spare parts cost, downtime loss cost, and expected failure risk cost calculated based on the probability density function of remaining useful life, including: constructing a comprehensive utility function: ,in, To the time The cumulative production revenue generated by the robot's normal operation For spare parts costs, To cover planned downtime losses; For the expected failure risk cost, ,in, Costs of serious failures caused by unplanned downtime. Let be the probability density function of the remaining useful life.

[0015] Preferably, the step of determining the maintenance time that maximizes the utility function by optimizing within a preset discrete maintenance window includes: dividing a preset future time period into multiple discrete maintenance windows; calculating the comprehensive utility function value when performing maintenance in each maintenance window; comparing all calculated comprehensive utility function values, and selecting the window corresponding to the maximum value as the optimal maintenance time.

[0016] The beneficial effects of this invention are as follows: By combining nonlinear kernel principal component analysis and Mahalanobis distance, this invention can identify subtle deviations in the operating state of industrial robots, thus improving the accuracy of state assessment. By distinguishing the robot's task types and establishing corresponding baseline risk functions, and simultaneously inputting a time-varying covariate that reflects historical accumulated damage, a remaining life prediction model that comprehensively considers the current state, historical damage, and operating conditions is constructed, improving the accuracy of remaining life prediction. By establishing a comprehensive utility function that includes spare parts, downtime losses, and expected failure risks calculated based on the life probability distribution, and optimizing this function, the optimal maintenance timing is determined, achieving optimal allocation of maintenance resources and ensuring production continuity. Attached Figure Description

[0017] Figure 1 The flowchart illustrating the steps of the industrial robot condition monitoring and predictive maintenance method in this embodiment is shown in the diagram. Detailed Implementation

[0018] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0019] like Figure 1 As shown, the industrial robot condition monitoring and predictive maintenance method in this embodiment includes the following steps: Step S1: Obtain multidimensional operating parameters and working environment sensing data of the industrial robot, and fuse them into a high-dimensional feature vector after time-domain alignment; perform kernel principal component analysis on the high-dimensional feature vector, project the vector of the current state onto the nonlinear principal component space, calculate the Mahalanobis distance between the projection point and the centroid of the preset healthy state cluster, and use the distance value as the first state covariate.

[0020] Specifically, encoders, current sensors, torque sensors deployed on the robot's six joints, and accelerometers installed on the end effector are used to collect real-time data on the angle, angular velocity, drive current, output torque, and end effector vibration signals of each joint. Temperature and humidity data of the robot's working unit are obtained through environmental sensors. The above multi-source heterogeneous data are resampled to a uniform 100Hz frequency using linear interpolation to achieve time-domain alignment. The aligned 38-dimensional data (covering joint states, end effector vibration, and environmental parameters, etc.) are then concatenated to form a high-dimensional feature vector.

[0021] Continuous operational data of the robot during typical tasks, whether brand new or immediately after maintenance, are collected to construct a health state sample set. Radial basis functions are used as the kernel function to perform kernel principal component analysis (KPCA) on this sample set, extracting the top K nonlinear principal components with a cumulative variance contribution rate exceeding 95%, thus constructing a health state space. The geometric center of the projection points of all healthy samples within this space is calculated as the centroid of the health state cluster. During the monitoring phase, the newly acquired high-dimensional feature vectors are projected onto the health state space using the established KPCA model. The Mahalanobis distance between the projection points and the centroid of the health state cluster is calculated; this distance quantifies the degree to which the current state deviates from the normal state and is defined as the first state covariate.

[0022] In an optional embodiment, multidimensional operating parameters of the industrial robot and sensor data of the working environment are acquired and fused into a high-dimensional feature vector after time-domain alignment. This includes: acquiring angle, angular velocity, and torque data of the robot's six-axis joints and three-dimensional acceleration data of the end effector at a sampling frequency of 100Hz; acquiring temperature and humidity data of the working environment at a sampling frequency of 1Hz; upsampling the temperature and humidity data to 100Hz using linear interpolation, aligning it with the robot's operating parameters on the timestamp, and splicing them together to form a high-dimensional feature vector with a dimension of 23.

[0023] Specifically, the data acquisition system acquires the robot's core operating parameters in real time at a frequency of 100 times per second, including the angles, angular velocities, and torques of the robot's six joints (a total of 18 data streams), as well as the accelerations in the X, Y, and Z directions measured by sensors installed on the end effector (a total of 3 data streams). The above high-frequency data provides a foundation for the refined analysis of the robot's behavior.

[0024] Simultaneously, temperature and humidity, two environmental parameters of the work area, are collected at a frequency of 1Hz. To address the issue of inconsistent data sampling frequencies, linear interpolation is used to upsample the low-frequency temperature and humidity data. For example, if the temperature is measured at 25℃ in the first second and 26℃ in the second second, the temperature values ​​at 99 intermediate moments between the two are calculated using the difference, thus increasing the data frequency to 100Hz. Finally, at each aligned timestamp, the 21-dimensional robot operating parameters are concatenated with the interpolated 2-dimensional environmental data to construct a unified high-dimensional feature vector of dimension 23 for subsequent state analysis.

[0025] In an optional embodiment, kernel principal component analysis is performed on the high-dimensional feature vectors to project the vectors of the current state onto a nonlinear principal component space, and the Mahalanobis distance between the projected points and the centroids of a preset healthy state cluster is calculated. This includes: using a Gaussian radial basis function as the kernel function of the kernel principal component analysis; selecting several high-dimensional feature vectors collected during the initial stage of robot operation as a healthy state sample set, and calculating the geometric center of the sample set after being projected onto the first K principal component dimensions by kernel principal component analysis as the centroid; after projecting the high-dimensional feature vectors of the current state, the Mahalanobis distance between the projected points and the centroids in the K-dimensional principal component space is calculated.

[0026] Several high-dimensional feature vectors collected during the initial stage of robot operation are selected as the health status sample set, including: 5000 high-dimensional feature vectors collected within the first 100 hours of robot operation are selected as the health status sample set; the first K principal component dimensions are the first 3 principal component dimensions.

[0027] Specifically, the process of establishing a health status baseline model is as follows: 5000 23-dimensional feature vectors collected during the first 100 hours of operation after a new industrial robot is put into use (or after a major overhaul) are selected to form a health status sample set. Kernel Principal Component Analysis (KPCA) is performed on this sample set, with a Gaussian radial basis function selected as the kernel function, and the kernel parameters are... The value is set to 1.5. Through KPCA transformation, the original 23-dimensional feature space is mapped and reduced to a three-dimensional nonlinear principal component space composed of the first three principal components. This space effectively preserves and reflects the inherent nonlinear structural characteristics of the data. The geometric center of the projection points of these 5000 samples in the three-dimensional space is calculated and denoted as the centroid of the health state cluster. The coordinates of this centroid represent the robot's baseline health state.

[0028] During the real-time monitoring phase of continuous robot operation, for each newly acquired 23-dimensional feature vector, it is transformed using the established KPCA baseline model parameters and projected onto the three-dimensional principal component space to obtain a projection point representing the current state. Subsequently, the Mahalanobis distance between this projection point and the center of the healthy state cluster is calculated. The Mahalanobis distance comprehensively considers the covariance structure and distribution shape of the healthy data cluster, and its value can accurately quantify the degree of deviation of the current state from the baseline healthy state. The calculated Mahalanobis distance is defined as the first state covariance; the larger the value, the further the robot's operating state deviates from the healthy baseline, i.e., the higher the degree of degradation.

[0029] Step S2: Identify the robot's historical task sequence and classify the task types. Based on the historical data, construct a baseline risk function for each task type based on the Weibull distribution. Calculate the kurtosis index of the robot's joint drive current and perform cumulative summation of the kurtosis index with a preset forgetting factor within a sliding time window to obtain the second time-varying covariate.

[0030] Specifically, program call records and load setting information are first extracted from the robot controller log, and combined with historical joint trajectory data obtained through inverse kinematics. K-means clustering is then used to perform unsupervised classification of these historical task segments, thus categorizing the robot's operational modes into three typical task types: heavy-load handling, high-speed precision welding, and small-area picking and placing. For each identified task type, full lifecycle operational data or failure time data from multiple robots of the same model performing that task are collected. Maximum likelihood estimation is then used to fit the shape and scale parameters of the Weibull distribution, thereby constructing a dedicated baseline risk function for each task type.

[0031] Next, the drive current signals of each joint were acquired in real time at a sampling rate of 100Hz, and a sliding time window with a width of 100 data sampling points (i.e., a duration of 1 second) was set. At each sampling moment, the kurtosis index of the current signal within the current time window was calculated to capture the signal impact characteristics caused by early mechanical damage. To quantify the cumulative effect of damage and reflect the degree of historical wear, an exponentially weighted moving average method was used to process the kurtosis value sequence. Specifically, a forgetting factor was set to 0.99, and the cumulative damage value at the current moment was equal to the cumulative value at the previous moment multiplied by the forgetting factor, plus the kurtosis value calculated at the current moment. This calculated cumulative value is used as the second time-varying covariate to dynamically correct the remaining life prediction results.

[0032] In an optional embodiment, identifying the robot's historical task sequence and classifying task types includes: segmenting historical data into independent task segments based on the robot's unique task program identifier and start timestamp; extracting the average joint angular velocity and average load torque of each task segment as features; and using the K-means clustering algorithm to cluster the features, classifying the historical tasks into three types: heavy-load handling, light-load assembly, and unloaded positioning.

[0033] Based on historical data, a baseline risk function based on the Weibull distribution is constructed for each task type. This includes: collecting historical failure time data from multiple robots of the same model under each task type; and fitting the shape parameters of the Weibull distribution using maximum likelihood estimation. and scale parameters According to the formula Construct the baseline risk function, where, This is the baseline risk function.

[0034] Specifically, the historical operation logs of the robot controller are analyzed, and the continuous operation data stream is divided into several independent task execution segments based on the task program number and task start timestamp. For each task segment, two key feature indicators are extracted: the average angular velocity and average load torque of all joints. Subsequently, the K-means clustering algorithm is applied to classify all historical task segments into three typical categories based on the above two feature indicators: for example, tasks with high average speed and average torque are classified as "heavy-load handling"; tasks with low speed and torque are classified as "light-load assembly"; and tasks with high speed but low torque are classified as "no-load positioning".

[0035] After classifying the tasks, for each specific task type (taking "heavy-load handling" as an example), we collect runtime data for all historical tasks belonging to that type up to the point of failure. Using this set of failure time data, we fit the Weibull distribution parameters, i.e., the shape parameters, that best characterize the distribution of this set of data using the maximum likelihood estimation method. and scale parameters For example, calculations may yield results for heavy-duty handling tasks. The value is 2.5. The value is 8000 hours. Finally, the estimated parameters are substituted into the standard formula of the Weibull risk function to construct a baseline risk function specific to this task type. It quantitatively describes the probability density or risk rate of inherent failures of equipment under ideal conditions in this specific task due to the passage of time.

[0036] In an optional embodiment, the kurtosis index of the robot joint drive current is calculated, and the kurtosis index is cumulatively summed with a preset forgetting factor within a sliding time window to obtain a second time-varying covariate, including: setting the width and sliding step size of the sliding time window; and calculating the kurtosis index of the drive current signal within each time window. The second time-varying covariate is calculated using an exponentially weighted cumulative summation method, and the calculation formula is as follows: ,in, Let be the second time-varying covariate at the current time t. For window width, The preset forgetting factor, and .

[0037] The sliding time window has a width of 1000 data sampling points, and the sliding step size is 100 data sampling points. The forgetting factor... It is 0.98.

[0038] Specifically, a sliding time window with a width of 1000 data points (approximately corresponding to a duration of 10 seconds) is set to analyze the joint drive current signal. This time window slides forward in steps of 100 data sampling points (i.e., 1 second), thereby achieving continuous monitoring of the current signal. Within each time window, the kurtosis index of the current signal is calculated. As a high-order statistic, kurtosis is highly sensitive to transient impacts or abnormal spikes in the signal, and therefore can effectively detect weak abnormal current fluctuations caused by early mechanical wear, gear pitting, or precursory faults.

[0039] To construct a health indicator that can both smooth out random fluctuations and sensitively reflect recent deterioration trends, an exponentially weighted cumulative algorithm with a forgetting factor is adopted.

[0040] Preferably, the forgetting factor The value is set to 0.98. This parameter indicates that when calculating the second time-varying covariate at the current moment, the kurtosis value calculated recently is given a higher weight, while the influence of earlier historical data on the current indicator decays exponentially over time. Through the weighted accumulation mechanism, the resulting covariate can not only respond promptly to newly emerging fault impact signals, but also effectively suppress the interference of isolated noise points by utilizing the smoothing effect of historical information, thereby accurately characterizing the cumulative damage state of the equipment.

[0041] Step S3: Construct a proportional hazards model, with the first state covariate and the second time-varying covariate as inputs. Select the corresponding baseline hazard function according to the task type currently being performed by the robot. Use the proportional hazards model to output the probability density function of the robot's remaining service life and calculate the expected value of the remaining service life.

[0042] Specifically, a proportional hazards model is constructed for remaining life prediction. The risk function of this model is defined as the product of the baseline risk function for a specific task type and the covariate influencing factors. The covariate influencing factors are in exponential form, and their exponential terms are composed of a linear combination of a first state covariate (reflecting the deviation from the current state) and a second time-varying covariate (reflecting historical cumulative damage).

[0043] During real-time prediction, the system executes the following steps: First, based on the task type currently identified by the robot, such as heavy-load handling, the corresponding baseline risk function is called from the pre-built function library; second, the first and second covariates calculated in real time are substituted into the model to calculate the robot's current instantaneous risk rate function; third, based on this risk rate function, an integral operation is performed to derive the survival probability function of the device at various future time points, and the derivative of the survival probability function is further calculated to obtain the probability density function of the remaining service life; finally, by integrating this probability density function to solve for the mathematical expectation, the expected value of the robot's remaining service life is calculated.

[0044] Step S4: Define a utility function that includes spare parts cost, downtime loss cost, and expected failure risk cost calculated based on the probability density function of remaining service life; determine the maintenance timing that maximizes the utility function by optimizing within a preset discrete maintenance window, and generate a maintenance instruction.

[0045] Specifically, a utility function is defined that includes spare parts cost, downtime loss cost, and expected failure risk cost calculated based on the probability density function of remaining useful life, including: constructing a comprehensive utility function: ,in, To the time The cumulative production revenue generated by the robot's normal operation For spare parts costs, To cover planned downtime losses; For the expected failure risk cost, ,in, Costs of serious failures caused by unplanned downtime. Let be the probability density function of the remaining useful life.

[0046] To determine the optimal timing for maintenance, a comprehensive utility model is established to maximize the net benefits that can be obtained by performing maintenance at different future time points. The net benefit is defined as the difference between the revenue generated from continued production and the costs associated with each maintenance activity.

[0047] The specific components are as follows: First, there is the accumulated continuous production revenue as operating time increases. For example, let's say the robot can generate 5000 yuan in value for every hour it operates normally. Secondly, the later the maintenance plan is scheduled (i.e., the later the timing...), the greater the potential profit. The larger the value of the robot, the longer its accumulated operating time before maintenance, resulting in higher total production revenue. Offsetting this revenue are maintenance-related costs, which consist of three parts. The first part is the fixed spare parts cost, a constant independent of time. For example, the hardware cost of replacing a critical bearing is fixed at 5000 yuan. The second part is the planned downtime loss cost, referring to the production interruption loss caused by proactively scheduling maintenance work. For example, assuming a standard maintenance job takes 4 hours, and the downtime loss per hour is 3000 yuan, then this part of the cost is fixed at 12000 yuan, which is also a cost determined before the decision. The third part is the expected failure risk cost, which varies with the maintenance time point. The variable cost, which increases monotonically over time, represents the planned maintenance time. The anticipated losses from unplanned robot malfunctions prior to arrival are calculated as follows: The cost of a serious malfunction due to unplanned downtime is set at a high value (e.g., 80,000 yuan), and this value is multiplied by the time... The cumulative probability of previous equipment failures. This cumulative probability is calculated using the remaining lifetime probability density function predicted by a proportional hazards model within the range... This is obtained by integrating the data. The later the maintenance plan is implemented, the higher the cumulative probability of failure, and the greater the risk cost associated with this aspect.

[0048] In an optional embodiment, the maintenance timing that maximizes the utility function is determined by optimizing within a preset discrete maintenance window, including: dividing a preset future time period into multiple discrete maintenance windows; calculating the comprehensive utility function value when performing maintenance in each maintenance window; comparing all calculated comprehensive utility function values, and selecting the window corresponding to the maximum value as the optimal maintenance timing.

[0049] First, a future decision interval is defined and divided into a series of discrete maintenance opportunity points, called maintenance windows. In this embodiment, the next 30 days are set as the decision interval, with "hours" as the smallest time unit. Therefore, this 30-day decision interval corresponds to 720 discrete candidate maintenance windows.

[0050] The optimization process involves evaluating each of the 720 candidate time points. The system starts from the first window and iterates through all 720 windows, performing optimizations for each maintenance window. ( Assuming maintenance is performed at that moment, the corresponding net utility value is calculated using the aforementioned comprehensive utility model. The program calculates the net utility value for the first hour, the net utility value for the second hour, and so on, until the net utility value for the 720th hour is calculated, thus generating a sequence list containing 720 utility values.

[0051] Finally, the system compares and analyzes the list of utility values, searching for and selecting the window that maximizes the utility function as the optimal maintenance time. For example, if the calculation finds that the net utility value of maintenance performed at the 158th hour is the highest among all 720 net utility values, the system determines the 158th hour as the optimal maintenance time. Based on this calculation result, the system automatically generates a maintenance instruction that clearly specifies which robot should be maintained, when (i.e., at the 158th hour), and what maintenance operation should be performed, thereby automating and refining maintenance decisions.

[0052] In the description of this specification, "multiple" means at least two, such as two, three or more, etc., unless otherwise expressly and specifically defined.

[0053] While various embodiments of the invention have been shown and described in this specification, it will be apparent to those skilled in the art that such embodiments are provided by way of example only. Many modifications, alterations, and alternatives will occur to those skilled in the art without departing from the spirit and essence of the invention.

Claims

1. A method for condition monitoring and predictive maintenance of industrial robots, characterized in that, Includes the following steps: Multidimensional operating parameters and working environment sensor data of industrial robots are acquired and fused into high-dimensional feature vectors through time-domain alignment; kernel principal component analysis is performed on the high-dimensional feature vectors, the vector of the current state is projected onto the nonlinear principal component space, the Mahalanobis distance between the projection point and the centroid of the preset healthy state cluster is calculated, and the distance value is used as the first state covariate; Identify the robot's historical task sequences and classify task types, and construct a baseline risk function based on the Weibull distribution for each task type based on historical data; Calculate the kurtosis index of the robot joint drive current, and sum the kurtosis index with a preset forgetting factor within a sliding time window to obtain the second time-varying covariate; A proportional hazards model is constructed, with the first state covariate and the second time-varying covariate as inputs. A corresponding baseline hazard function is selected based on the task type currently being performed by the robot. The proportional hazards model is used to output the probability density function of the robot's remaining service life, and the expected value of the remaining service life is calculated. Define a utility function that includes spare parts cost, downtime loss cost, and expected failure risk cost calculated based on the probability density function of remaining service life; determine the maintenance timing that maximizes the utility function by optimizing within a preset discrete maintenance window, and generate maintenance instructions.

2. The industrial robot condition monitoring and predictive maintenance method according to claim 1, characterized in that, The acquisition of multidimensional operating parameters and working environment sensor data of the industrial robot, fused into a high-dimensional feature vector through time-domain alignment, includes: The robot's six-axis joint angle, angular velocity, and torque data, as well as the three-dimensional acceleration data of the end effector, are collected at a sampling frequency of 100Hz. The temperature and humidity data of the working environment are collected at a sampling frequency of 1Hz. The temperature and humidity data are upsampled to 100Hz using linear interpolation, aligned with the robot's operating parameters on the timestamp, and spliced ​​together to form a high-dimensional feature vector with a dimension of 23.

3. The method for industrial robot condition monitoring and predictive maintenance according to claim 1, characterized in that, The step of performing kernel principal component analysis on the high-dimensional feature vector, projecting the vector of the current state onto the nonlinear principal component space, and calculating the Mahalanobis distance between the projected point and the centroid of the preset healthy state cluster includes: Gaussian radial basis function is used as kernel function for kernel principal component analysis; several high-dimensional feature vectors collected in the initial stage of robot operation are selected as health state sample set, and the geometric center of the sample set after being projected onto the first K principal component dimensions by kernel principal component analysis is calculated as centroid; after projecting the high-dimensional feature vector of the current state, the Mahalanobis distance between the projection point and the centroid in the K-dimensional principal component space is calculated.

4. The industrial robot condition monitoring and predictive maintenance method according to claim 3, characterized in that, The selection of several high-dimensional feature vectors collected during the initial stage of robot operation as a health status sample set includes: Five thousand high-dimensional feature vectors collected during the first 100 hours of robot operation were selected as the health status sample set. The first K principal component dimensions are the first 3 principal component dimensions.

5. The method for industrial robot condition monitoring and predictive maintenance according to claim 1, characterized in that, The process of identifying the robot's historical task sequences and classifying task types includes: Based on the unique identifier and start timestamp of the robot task program, historical data is divided into independent task segments; the average joint angular velocity and average load torque of each task segment are extracted as features; the K-means clustering algorithm is used to cluster the features, and the historical tasks are clustered into three types: heavy-load handling, light-load assembly, and unloaded positioning.

6. The method for industrial robot condition monitoring and predictive maintenance according to claim 1, characterized in that, The construction of a baseline risk function based on the Weibull distribution for each task type based on historical data includes: For each task type, historical failure time data of multiple robots of the same model under that task type were collected; the shape parameters of the Weibull distribution were fitted using the maximum likelihood estimation method. and scale parameters According to the formula Construct the baseline risk function, where, This is the baseline risk function.

7. The industrial robot condition monitoring and predictive maintenance method according to claim 1, characterized in that, The calculation involves determining the kurtosis index of the robot joint drive current and summing the kurtosis index with a preset forgetting factor within a sliding time window to obtain a second time-varying covariate, including: Set the width and step size of the sliding time window; within each time window, calculate the kurtosis index of the drive current signal. The second time-varying covariate is calculated using an exponentially weighted cumulative summation method, and the calculation formula is as follows: ,in, Let be the second time-varying covariate at the current time t. For window width, The preset forgetting factor, and .

8. The method for industrial robot condition monitoring and predictive maintenance according to claim 7, characterized in that, The sliding time window has a width of 1000 data sampling points, and the sliding step size is 100 data sampling points. The forgetting factor... It is 0.

98.

9. The method for industrial robot condition monitoring and predictive maintenance according to claim 1, characterized in that, The definition includes a utility function of spare parts cost, downtime loss cost, and expected failure risk cost calculated based on the probability density function of remaining useful life, including: Construct a comprehensive utility function: ,in, To the time The cumulative production revenue generated by the robot's normal operation For spare parts costs, To cover planned downtime losses; For the expected failure risk cost, ,in, Costs of serious failures caused by unplanned downtime. Let be the probability density function of the remaining useful life.

10. The method for industrial robot condition monitoring and predictive maintenance according to claim 1, characterized in that, The step of determining the maintenance timing that maximizes the utility function by optimizing within a preset discrete maintenance window includes: The future preset time period is divided into multiple discrete maintenance windows; for each maintenance window, the comprehensive utility function value when performing maintenance in that window is calculated; all calculated comprehensive utility function values ​​are compared, and the window corresponding to the maximum value is selected as the optimal maintenance time.

Images