Modeling method for predicting influence of industrial robot introduction on employment based on xgboost model

Abstract

The invention relates to the technical field of employment prediction, in particular to a prediction modeling method for influence of industrial robot introduction on employment based on an xgboost model, and the method comprises the following steps: 1, data preprocessing: constructing a high-dimensional feature space through k-NN interpolation, Robust Scaler normalization and comprehensive time sequence feature engineering; step 2, data inspection and data clustering: PCA anomaly detection and DTW time sequence K-Means clustering are used to provide guarantee for data quality and heterogeneity understanding; step 3, establishing a prediction model: establishing an XGBoost prediction model based on a time sequence; and 4, establishing an early warning framework and conclusion based on the prediction model. According to the method, the robot density of each province is successfully predicted by using the XGBoost model, the model shows relatively good prediction performance on a test set, and the accuracy rate can reach 88.34%; according to the method, the space-time dynamic, the regional mode and the future trend adopted by the Chinese provincial robot can be deeply analyzed, and the potential employment influence can be evaluated.

Description

Technical Field

[0001] This invention relates to the field of employment forecasting technology, and in particular to a predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model. Background Technology

[0002] As a classic application of artificial intelligence, industrial robots have not only driven an efficiency revolution in traditional industries but also profoundly impacted employment structures. However, the release of technological dividends is accompanied by structural contradictions. The large-scale replacement of traditional labor by robots has, on the one hand, indirectly increased employment in R&D and management; on the other hand, it has directly reduced production employment. How to balance this substitution-creation impact on employment is a major issue for the technological upgrading of the manufacturing industry, and further research is needed into the relationship between industrial robots and manufacturing employment. The relationship between industrial robot investment and employment exhibits significant regional differences, and modeling methods for predictive modeling research are still lacking. Summary of the Invention

[0003] The purpose of this invention is to address the shortcomings of existing technologies by proposing a predictive modeling method for the impact of industrial robot adoption on employment based on the XGBoost model. This method can deeply analyze the spatiotemporal dynamics, regional patterns, and future trends of robot adoption at the provincial level in China, and assess its potential employment impact.

[0004] To achieve the above objectives, the present invention adopts the following technical solution:

[0005] A predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model includes the following steps:

[0006] Step 1, Data Preprocessing: A high-dimensional feature space is constructed through k-NN interpolation, RobustScaler normalization, and comprehensive temporal feature engineering;

[0007] Step 2, Data Validation and Clustering: PCA anomaly detection and DTW time series K-Means clustering are used to ensure data quality and understanding of heterogeneity;

[0008] Step 3: Establish a prediction model: Establish an XGBoost prediction model based on time series data;

[0009] Step 4: Establish an early warning framework and conclusions based on the prediction model.

[0010] Preferably, step 1 includes the following method:

[0011] Step 1.1: Data Preprocessing

[0012] To ensure that each region in the data table has complete annual data records within the time span, the k-nearest neighbor (k-NN) imputation method is adopted. This method is based on the idea of ​​clustering similar things, and finds the k nearest neighbor samples that are most similar to the sample containing missing values ​​in other known features, and uses the information of these neighbor samples to estimate the missing values.

[0013] (1) The distance metric uses cosine similarity to measure samples, assuming it is the feature vector x of sample i and sample j. i and x j The degree of similarity between vectors, with cosine similarity focusing on vector direction rather than absolute magnitude, is suitable for handling data with features of different magnitudes.

[0014]

[0015] Similarity(x i ,x j ) represents the cosine similarity between sample i and sample j, used to measure the degree of similarity between the feature vectors of the two samples. i x j Let x be the feature vector of sample i and sample j respectively. i .x j It is the dot product of the feature vectors of sample i and sample j, ||x i ||,||x j || L2 norm x of the eigenvectors of sample i and sample j ik and x jk The values ​​of samples i and j on the k-th feature. n is the total number of features in the samples (i.e., the dimension of the feature vector).

[0016] (2) Nearest neighbor selection: After calculating the similarity between samples, select the k samples that are closest to the target sample, i.e., the samples with the highest similarity, as the nearest neighbor set N. K After testing and evaluation, k=7 was chosen as the optimal number of nearest neighbors.

[0017] N k (x new )={(x j ,y j )|j∈I k}

[0018] N k (x new ) represents the target sample x with missing values. new The set consisting of the k nearest neighbors. new It is the target sample with missing values. j The feature vector of the j-th sample in the nearest neighbor set. jThe variable value of the j-th sample in the nearest neighbor set can be understood as the feature information of that sample. k The index set containing these k nearest neighbor samples is used to mark which samples are the nearest neighbors of the target sample.

[0019] (3) Imputation strategy for missing values ​​in numerical data The filling is performed using the arithmetic mean of the feature values ​​corresponding to its k nearest neighbors, calculated as follows:

[0020]

[0021] The estimated value of the missing feature in the target sample, i.e., the final result used to fill in the missing values. k is the number of nearest neighbor samples selected. N k The set of the k nearest neighbors of the target sample. x' j The value of the j-th sample in the nearest neighbor set for the missing feature.

[0022] By implementing k-NN interpolation, gaps in the dataset are effectively filled, the integrity of time series data is maintained, and a more reliable data foundation is provided for subsequent analysis.

[0023] Step 1.2, Data Normalization

[0024] The dataset contains various indicators such as employment numbers and robot density across different industries, with significant differences in units and numerical ranges. To eliminate these dimensional differences and ensure that the features are on a comparable scale, this invention employs the RobustScaler method for data normalization. RobustScaler is a normalization technique insensitive to outliers. It scales the data by subtracting the median of the features and then dividing by the interquartile range (IQR), as shown in the following formula:

[0025]

[0026] X represents the original feature data to be processed. scaled The scaled feature data after RobustScaler normalization. Median(X) is the median of all sample values ​​for this feature. IQR(X) represents the interquartile range of this feature.

[0027] The interquartile range (IQR) is defined as follows:

[0028] IQR(X) = Q3(X) - Q1(X)

[0029] Q3(X) represents the upper quartile (third quartile) of the feature (the value at the 75th percentile after sorting the feature data from smallest to largest); Q4(X) represents the lower quartile (first quartile) of the feature (the value at the 25th percentile after sorting the feature data from smallest to largest).

[0030] Derivative feature construction and feature engineering

[0031] To more deeply capture time-series relationships in the data, smooth short-term fluctuations, and extract dynamic patterns, a series of derived features were constructed based on the preprocessed base features. These features aim to enhance the ability of subsequent models to capture complex time dependencies. Feature engineering in this phase includes the following aspects; all operations are performed in groups based on data sorted by province and year to ensure the correctness of time-series operations.

[0032] (1) Smoothing

[0033] To reduce short-term noise and random fluctuations in the data, thereby revealing potential long-term trends more clearly, a simple moving average (SMA) was applied to the target variable (robot density) and all six selected industry characteristic variables (agriculture, manufacturing, build, transaction, it, education). The time window size (SMOOTHING_WINDOW) was set to 3 periods, allowing for fewer than 3 data points in the initial window (min_periods=1) to avoid generating too many missing values ​​at the beginning of the sequence.

[0034]

[0035] SMA t X is the simple moving average at time t (a smoothed value at time t calculated using historical data within the window). W is the size of the moving window and the number of historical data points used in calculating the moving average. t-i Let ti be the value of variable X (the original data of the i-th historical time within the window).

[0036] (2) Lag characteristics

[0037] To capture the time autocorrelation of variables, i.e., how current values ​​may be influenced by their past values. This is crucial for time series forecasting models. For the target variable, robot density, and all six industry characteristic variables, lagged values ​​for the past one and two periods were created, with a maximum lag order N_LAGS set to 2.

[0038] (3) Difference characteristics

[0039] To eliminate trends in the time series and make it closer to a stationary series, the model can analyze the magnitude of changes in variables rather than their absolute levels. A first-order difference was calculated for the target variable, robot density, which is the value at the current time point minus the value at the previous time point. The difference order is 1, generating the feature column `robot_density_diff1`. According to the code implementation, this processing only performs differencing on the target variable.

[0040] ΔX t =X t -X t-1

[0041] ΔX t： The backward difference of variable X at time t describes the magnitude of the change in variable X at time t relative to time t-1. t Let X be the value of variable X at time t (the raw data of the 0th historical time within the window). t-1 Let X be the value of variable X at time t-1 (the original data of the first historical time in the window). This formula is defined by backward difference.

[0042] (4) Window statistical characteristics

[0043] Capture the local dynamic characteristics of variables over a recent period, such as central tendency, volatility, and range of values; set the rolling window size to three periods, allowing for fewer than three data points in the initial stage of the window; calculate the following statistics: mean, standard deviation, minimum, and maximum.

[0044]

[0045] σ roll,t μ represents the rolling standard deviation at time t, used to measure the degree of fluctuation of a variable within a rolling window ending at time t. roll,t This represents the average value of variable X within the scrolling window (of all X values ​​within the window). t-i (arithmetic mean).

[0046] (5) Handling missing values ​​introduced by derived features

[0047] The aforementioned lag, differencing, and windowing operations naturally generate missing NaN values ​​at the beginning of each province's time series. To ensure data quality for subsequent analysis and modeling, all rows containing any NaN values ​​are removed from the dataset.

[0048] Preferably, step 2 includes the following method:

[0049] Data verification

[0050] (1) PCA Algorithm

[0051] Principal Component Analysis (PCA) is a commonly used multivariate statistical method. Besides its core dimensionality reduction function, it is also frequently used for outlier detection. Its basic idea is that most data points should follow the data structure and pattern defined by the principal components, while outliers may significantly deviate from this pattern. By examining the error in projecting data points from the original space to the principal component space and then reconstructing them back to the original space, "maverick" observations can be effectively identified.

[0052] Using PCA technology, a whole scan is performed on the preprocessed dataset containing rich derived features to identify potentially abnormal province-year data points. Based on the detection results, the overall usability of the data is evaluated, providing a basis for subsequent clustering and predictive modeling.

[0053] (2) Anomaly detection based on reconstruction error

[0054] The reconstruction error method is used to identify outliers. The reconstruction error is calculated as follows: for each data point, the mean square reconstruction error is calculated on all standardized features.

[0055]

[0056] MSE i Let be the mean squared reconstruction error of the i-th data point, used to measure the degree of deviation of the data point after PCA reconstruction (the larger the deviation, the more likely it is an outlier). p represents the total number of features, and j is the index variable for summation, used to iterate through each feature dimension. scaled,ij x represents the standardized value (normalized feature value) of the i-th data point on the j-th feature. reconstructed,ij This represents the PCA reconstruction value of the i-th data point on the j-th feature (the feature value obtained by dimensionality reduction and restoration through PCA).

[0057] Determine the threshold: Select all data points MSE i The 99th percentile of the value is used as the threshold for judging outliers; outlier identification: if the MSE of a data point is... i If the value is greater than the threshold, it is marked as a potential outlier.

[0058] (3) Visualization of abnormal data

[0059] To more intuitively understand the distribution of outliers and the setting of thresholds, visual charts are created;

[0060] (4) Construction of K-means model

[0061] To gain a deeper understanding of the similarities and differences in economic development and industrial structure evolution among my country's provinces, cluster analysis is employed. Given that the data in this invention spans multiple years, simply clustering at a single point in time or summarizing statistics may lose important dynamic information. Therefore, the Time Series K-Means algorithm is chosen. This algorithm is used to cluster time series data, aiming to find K cluster centers, grouping the time series data of each province into X... i Assigned to the only cluster C K In this case, all sample points within a cluster are such that they are all within the cluster center M. K Minimize the sum of squares of the distances;

[0062]

[0063] Time series K-Means can directly process the multidimensional feature sequences of each province over time. To better capture the nonlinear alignment or temporal lag or lead of development pace in different provinces, Dynamic Time Warping (DTW) is chosen as the distance metric for time series K-Means. DTW finds the distance between two time series X = (x 1, x2,K x n ) and Y = (y1, y 2, K y m The optimal normalized path between points () is calculated, and their similarity is determined. This path defines a point-to-point mapping relationship that minimizes the cumulative distance between corresponding points. The DTW distance is usually calculated using dynamic programming, and its core recursive relationship is as follows:

[0064] D(i,j)=cost(x i y j )+min\begincasesD(i-1,j)

[0065] D(i,j-1)

[0066] D(i-1,j-1)\endcases

[0067] Where cost(x) i ,y i () is point x i and y i The distance between them, D(i,j) is the cumulative DTW distance between subsequences X[1...i] and Y=[1...j], and the final d DTW (x,y) is equivalent to D(n,m). DTW allows for non-linear scaling of the time axis and can identify sequences with similar shapes but may be offset in time.

[0068] (5) Determine the optimal number of clusters

[0069] Choosing the appropriate number of clusters K is a crucial step in the K-Means algorithm. Two commonly used methods assist in the decision-making process: the Elbow Method calculates the sum of squared intra-cluster distances for different K values. Inertia measures the compactness of samples within a cluster and typically decreases as K increases; finding the K value corresponding to the "elbow" where the rate of Inertia decreases from steep to gentle is considered a better choice. For each K value, the silhouette coefficient is calculated as the average silhouette coefficient of the clustering results. The silhouette coefficient comprehensively considers the compactness of samples within clusters and the separation of samples between clusters, with a value ranging from -1 to 1. A score closer to 1 indicates a better clustering effect. The K value that maximizes the average silhouette coefficient is selected as the candidate. When calculating the silhouette coefficient, the DTW distance matrix between all provincial time series is pre-calculated and passed to the silhouettescore function to ensure that the evaluation metric is consistent with the distance metric used in the clustering algorithm.

[0070] (6) Optimal clustering results

[0071] Using the determined optimal K value of 4, the final DTW clustering is performed on the multidimensional time series data of the province. The algorithm running parameters include setting n_init to max(1, N_INIT, KMEANS / / 2); max_iter to max(10, MAX_ITER, KMEANS / / 2); and random state to RANDOMSTATE to ensure that the results are reproducible.

[0072] (7) Clustering visualization

[0073] Cluster profile visualization: To visually demonstrate the differences in key economic structural features among different clusters, we first calculate the time average of all provinces within each cluster on the selected feature (agri manu build trans it edu), and then draw grouped bar charts for comparison.

[0074] Preferably, step 3 includes the following method:

[0075] 1. Features and target variables

[0076] Target variable: The prediction target is the robot installation density in each province.

[0077] Input features: The input features used to train the XGBoost model include all the numerical features generated in the first preprocessing stage. Except for the identifier column (province PROVINCECOLUMN and year YEAR COLUMN) and the original target variable (TARGETRENAME) itself, all features are used as input. This includes all derived features generated from the original 6 industry indicators (agri manubuild trans it edu) through smoothing, lag, differencing and rolling statistics, especially the lag term and rolling statistics features, which are designed to help the model capture the dynamic dependencies of time series.

[0078] Time series data partitioning

[0079] To simulate real-world prediction scenarios when evaluating model performance and to avoid data leakage, the dataset is strictly divided in chronological order: training set, validation set, and test set.

[0080] Model Configuration

[0081] 2. Model Training

[0082] The XGBRegressor interface of XGBoost is used to build a regression model. The model is trained using the training set (X train, Y train) and monitored for early stopping using the validation set (X val, Y val). During training, the evaluation metrics (RMSE, MAE) on the validation set are output periodically. Finally, the model state that performs best on the validation set is selected. The trained model is saved to a file for later use.

[0083] 3. Model Prediction

[0084] Using a trained XGBoost model, the robot density for the next year (taking 2022 as an example) is predicted;

[0085] Model Validation

[0086] After the model training is completed, its final performance is evaluated using reserved test set data from 2020 and later, which the model has never seen before.

[0087] Preferably, step 4 includes the following method:

[0088] By correlating the predicted growth in robot density with the changing trends in employment in related industries, we can identify regions where automation is progressing rapidly and employment in related industries may be under pressure, and issue early warning signals so that relevant departments can formulate response strategies in advance.

[0089] Compared with the prior art, the present invention has the following beneficial effects:

[0090] 1. This invention successfully predicted the robot density of each province using the XGBoost model. The model showed good prediction performance on the test set, with an accuracy of up to 88.34%.

[0091] 2. This invention can deeply analyze the spatiotemporal dynamics, regional patterns and future trends of robot adoption at the provincial level in China, and assess its potential employment impact. Attached Figure Description

[0092] Figure 1 This is a visualization of the two principal components of the present invention;

[0093] Figure 2 This is the reconstructed error histogram for this invention;

[0094] Figure 3 To reconstruct the error box plot for this invention;

[0095] Figure 4 This is an example diagram of abnormal and normal points in this invention;

[0096] Figure 5 This is the ED vs DTW diagram of the present invention;

[0097] Figure 6 This is a diagram showing the selection of the number of clusters in this invention;

[0098] Figure 7 This is a bar chart showing the cluster grouping of the present invention;

[0099] Figure 8 This is a graph showing the average characteristic values ​​of the education industry in this invention.

[0100] Figure 9 This is the cluster center diagram (cluster 0, 1) of the present invention;

[0101] Figure 10 This is the cluster center diagram (cluster 2, 3) of the present invention;

[0102] Figure 11 This is a pie chart showing the proportion of each cluster in this invention;

[0103] Figure 12 This is a diagram of the xgboost model of the present invention;

[0104] Figure 13 This is a sorting bar chart for the present invention;

[0105] Figure 14 This is a partial diagram showing the importance of features of the present invention.

[0106] Figure 15 This is a graph showing the predicted values ​​versus the actual values ​​in this invention.

[0107] Figure 16 This is the residual diagram of the present invention. Detailed Implementation

[0108] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings, so that those skilled in the art can better understand the advantages and features of the present invention, thereby making a clearer definition of the scope of protection of the present invention. The embodiments described in this invention are only some embodiments of the present invention, not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0109] Example:

[0110] A predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model includes the following steps:

[0111] Step 1: Data preprocessing: A high-dimensional feature space was constructed through k-NN interpolation, RobustScaler normalization, and comprehensive temporal feature engineering;

[0112] 1. Data Preprocessing

[0113] To ensure that each region in the data table has complete annual data records for the time span, the k-nearest neighbor (k-NN) imputation method is used. This method is based on the idea of ​​clustering similar things, finding the k nearest neighbors that are most similar to the sample containing the missing value in other known features, and using the information from these neighbor samples to estimate the missing value.

[0114] (1) The distance metric uses cosine similarity to measure samples, assuming it is the feature vector x of sample i and sample j. i and x j The degree of similarity between vectors. Cosine similarity focuses on the direction of vectors rather than their absolute magnitude, making it suitable for handling data with features of different magnitudes.

[0115]

[0116] Similarity(x i ,x j ) represents the cosine similarity between sample i and sample j, used to measure the degree of similarity between the feature vectors of the two samples. i x j Let x be the feature vector of sample i and sample j respectively. i .x j It is the dot product of the feature vectors of sample i and sample j, ||x i ||,||x j || L2 norm x of the eigenvectors of sample i and sample j ik and x jkThe values ​​of samples i and j on the k-th feature. n is the total number of features in the samples (i.e., the dimension of the feature vector).

[0117] (2) Nearest neighbor selection: After calculating the similarity between samples, select the k samples that are closest to the target sample, i.e., the samples with the highest similarity, as the nearest neighbor set N. K After testing and evaluation, k=7 was chosen as the optimal number of nearest neighbors.

[0118] N k (x new )={(x j y j )|j∈I k}

[0119] N k (x new ) represents the target sample x with missing values. new The set consisting of the k nearest neighbors. new It is the target sample with missing values. j The feature vector of the j-th sample in the nearest neighbor set. j The variable value of the j-th sample in the nearest neighbor set can be understood as the feature information of that sample. k The index set containing these k nearest neighbor samples is used to mark which samples are the nearest neighbors of the target sample.

[0120] (3) Imputation strategy for missing values ​​in numerical data The filling is performed using the arithmetic mean of the feature values ​​corresponding to its k nearest neighbors, calculated as follows:

[0121]

[0122] The estimated value of the missing feature in the target sample, i.e., the final result used to fill in the missing values. k is the number of nearest neighbor samples selected. N k The set of the k nearest neighbors of the target sample. x' j The value of the j-th sample in the nearest neighbor set for the missing feature.

[0123] By implementing k-NN interpolation, gaps in the dataset are effectively filled, maintaining the integrity of the time series data and providing a more reliable data foundation for subsequent analysis.

[0124] 2. Data normalization

[0125] The dataset contains various indicators such as employment numbers and robot density across different industries, and these indicators vary greatly in scale and numerical range. To eliminate these dimensional differences and make the features comparable, this invention employs the RobustScaler method for data normalization.

[0126] RobustScaler is a normalization technique that is insensitive to outliers. It scales by subtracting the median of the feature and then dividing by the interquartile range (IQR).

[0127]

[0128] X represents the original feature data to be processed. scaled The scaled feature data after RobustScaler normalization. Median(X) is the median of all sample values ​​for this feature. IQR(X) represents the interquartile range of this feature.

[0129] The interquartile range (IQR) is defined as follows:

[0130] IQR(X) = Q3(X) - Q1(X)

[0131] Q3(X) represents the upper quartile (third quartile) of the feature (the value at the 75th percentile after sorting the feature data from smallest to largest); Q4(X) represents the lower quartile (first quartile) of the feature (the value at the 25th percentile after sorting the feature data from smallest to largest).

[0132] Derivative feature construction and feature engineering

[0133] To more deeply capture time-series relationships in the data, smooth short-term fluctuations, and extract dynamic patterns, this invention further constructs a series of derived features based on the preprocessed basic features. These features aim to enhance the ability of subsequent models to capture complex time dependencies. The feature engineering in this stage mainly includes the following aspects; all operations are performed in groups based on data sorted by province and year to ensure the correctness of time-series operations.

[0134] (1) Smoothing

[0135] To reduce short-term noise and random fluctuations in the data, thereby revealing potential long-term trends more clearly, a simple moving average (SMA) was applied to the target variable (robot density) and all six selected industry characteristic variables (agriculture, manufacturing, build, transaction, it, education). The time window size (SMOOTHING_WINDOW) was set to 3 periods. During calculation, fewer than 3 data points were allowed in the initial stage of the window (min_periods=1) to avoid generating too many missing values ​​at the beginning of the sequence.

[0136]

[0137] SMA t X is the simple moving average at time t (a smoothed value at time t calculated using historical data within the window). W is the size of the moving window and the number of historical data points used in calculating the moving average. t-i Let ti be the value of variable X (the original data of the i-th historical time within the window).

[0138] (2) Lag characteristics

[0139] To capture the time autocorrelation of variables—that is, how current values ​​may be influenced by their past values—is crucial for time series forecasting models. Lag values ​​for the target variable, robot density, and all six industry-specific variables were created for the past 1 and 2 periods, respectively. The maximum lag order N_LAGS was set to 2.

[0140] (3) Difference characteristics

[0141] To eliminate trends in the time series data and make it closer to a stationary series, the model can analyze the magnitude of changes in variables rather than their absolute levels. A first-order difference was calculated for the target variable, robot density, which is the value at the current time point minus the value at the previous time point. The difference order is 1. The feature column `robot_density_diff1` was generated. According to the code implementation, this processing only differencing the target variable.

[0142] ΔX t =X t -X t-1

[0143] ΔX t： The backward difference of variable X at time t describes the magnitude of the change in variable X at time t relative to time t-1. t Let X be the value of variable X at time t (the raw data of the 0th historical time within the window). t-1 Let X be the value of variable X at time t-1 (the original data of the first historical time in the window). This formula is defined by backward difference.

[0144] (4) Window statistical characteristics

[0145] This feature captures the local dynamic characteristics of variables over a recent period, such as central tendency, volatility, and range. The rolling window size is set to 3 periods, allowing for fewer than 3 data points at the beginning of the window. Calculated statistics include: mean, standard deviation, minimum, and maximum.

[0146]

[0147] σ roll,t μ represents the rolling standard deviation at time t, used to measure the degree of fluctuation of a variable within a rolling window ending at time t. roll,t This represents the average value of variable X within the scrolling window (of all X values ​​within the window). t-i (arithmetic mean).

[0148] (5) Handling missing values ​​introduced by derived features

[0149] The aforementioned lag, differencing, and windowing operations naturally introduce missing NaN values ​​at the beginning of each province's time series. To ensure data quality for subsequent analysis and modeling, all rows containing any NaN values ​​are removed from the dataset.

[0150] Step 2: Data Validation and Clustering: PCA anomaly detection and DTW time series K-Means clustering were used to ensure data quality and understanding of heterogeneity.

[0151] 1. Data validation

[0152] (1) PCA Algorithm

[0153] Principal Component Analysis (PCA) is a commonly used multivariate statistical method. Besides its core dimensionality reduction function, it is also frequently used for outlier detection. Its basic idea is that most data points should follow the data structure and pattern defined by the principal components, while outliers may significantly deviate from this pattern. By examining the error in projecting data points from the original space to the principal component space and then reconstructing them back to the original space, "maverick" observations can be effectively identified.

[0154] PCA technology is used to perform a comprehensive scan of the preprocessed dataset containing rich derived features, identifying potentially abnormal province-year data points. Based on the detection results, the overall usability of the data is evaluated, providing a basis for subsequent clustering and predictive modeling. Figure 1 As shown.

[0155] (2) Anomaly detection based on reconstruction error

[0156] The reconstruction error method is used to identify outliers. The reconstruction error is calculated as follows: for each data point, the mean squared reconstruction error is calculated across all standardized features.

[0157]

[0158] MSE i Let be the mean squared reconstruction error of the i-th data point, used to measure the degree of deviation of the data point after PCA reconstruction (the larger the deviation, the more likely it is to be an outlier). p represents the total number of features, and j is the index variable for summation, used to iterate through each feature dimension. scaled,ij x represents the standardized value (normalized feature value) of the i-th data point on the j-th feature. reconstructed,ij This represents the PCA reconstruction value of the i-th data point on the j-th feature (the feature value obtained by dimensionality reduction and restoration through PCA).

[0159] Determine the threshold: Select all data points MSE i The 99th percentile of the value is used as the threshold for judging outliers. Outlier identification: If the MSE of a data point... i If the value is greater than the threshold, it is marked as a potential outlier.

[0160] (3) Visualization of abnormal data

[0161] To more intuitively understand the distribution of outliers and the setting of thresholds, the present invention has created the following visualization charts.

[0162] First, this invention plotted the distribution of the reconstruction error MSE to understand its overall statistical characteristics. Histogram: This shows the frequency distribution of the reconstruction error for all data points, such as... Figure 2 As shown.

[0163] A PCA reconstruction error distribution histogram. The horizontal axis represents the magnitude of the reconstruction error (MSE), and the vertical axis represents the frequency. The graph shows that the error for most data points is concentrated at a low level, while the error for a few points is significantly higher, exhibiting a long-tailed distribution. The red dashed line marks the 99th percentile threshold.

[0164] Box plots provide another view of the reconstruction error distribution, making it easier to observe the median, interquartile range, and extreme values, such as... Figure 3 As shown.

[0165] PCA reconstruction error box plot. The boxes represent the interquartile range of the data, the median line indicates the median, and the whiskers extend to most of the data range. Points outside the boxes, especially those far beyond the whiskers, are generally considered potential outliers. The red dashed line also indicates the 99th percentile threshold, clearly showing the degree of deviation of points exceeding the threshold from the main data.

[0166] Secondly, in order to observe the location of outliers while preserving the main data structure, this invention projects the data onto a three-dimensional space composed of the first three principal components and uses Plotly Express to generate an interactive 3D scatter plot, such as... Figure 4 As shown.

[0167] The distribution of data points in the first three principal component spaces (PC1, PC2, PC3). Each point represents an observation from a province-year. The color and size of the points can be encoded according to their reconstruction error; the darker the color or the larger the point, the higher the reconstruction error and the more likely it is to be an outlier.

[0168] (4) Construction of K-means model

[0169] To gain a deeper understanding of the similarities and differences in economic development and industrial structure evolution among my country's provinces, cluster analysis is employed. Given that the data in this invention spans multiple years, simply clustering at a single point in time or summarizing statistics may lose important dynamic information. Therefore, the Time Series K-Means algorithm is chosen. This algorithm is used to cluster time series data, aiming to find K cluster centers, grouping the time series data of each province into X... i Assigned to the only cluster C K In this case, all sample points within a cluster are such that they are all within the cluster center M. K Minimize the sum of squares of the distances:

[0170]

[0171] J DTW C represents the objective function value of the K-Means time series dataset, which is the sum of squared DTW distances from all samples within a cluster to their corresponding cluster centers. The goal is to minimize this value to achieve compactness of samples within each cluster. K is the number of clusters. That is, the samples should be divided into K categories. k X represents the sample set corresponding to the k-th cluster, containing all time series samples assigned to the k-th cluster. i Let M be the i-th time series sample within the k-th cluster, corresponding to a multidimensional feature sequence of a province over time. k d represents the center of the k-th cluster, and is the "average" time series of all time series samples within that cluster, serving as the representative of the cluster. DTW (X i M k ) is the sample X i With cluster center M k The DTW distance between them.

[0172] Time series K-Means can directly process the multidimensional feature sequences of each province over time. To better capture the nonlinear alignment or temporal lag or lead of development pace in different provinces, Dynamic Time Warping (DTW) is chosen as the distance metric for time series K-Means. DTW finds the distance between two time series X = (x1, x2, K x...) n ) and Y = (y1, y2, Ky) m The optimal normalized path between points () is calculated, and their similarity is computed. This path defines a point-to-point mapping relationship that minimizes the cumulative distance between corresponding points. DTW distance is typically calculated using dynamic programming, with the core recursive relationship as follows:

[0173] D(i,j)=cost(x i y j )+min\begincasesD(i-1,j)

[0174] D(i,j-1)

[0175] D(i-1,j-1)\endcases

[0176] D ij It is the cumulative distance between the i-th element of time series X and the j-th element of time series Y, where cost(x) is the cumulative distance between the i-th element of time series X and the j-th element of time series Y. i ,y i The distance between the i-th element of sequence X and the j-th element of sequence Y is used to measure the difference between the two individual data points. Begincases and endcases represent the start and end of the iteration, respectively.

[0177] Where cost(x) i ,y i () is point x i and y i The distance between them, D(i,j) is the cumulative DTW distance between subsequences X[1...i] and Y=[1...j]. The final d DTW (x,y) is equivalent to D(n,m). DTW allows for non-linear scaling of the time axis and can identify sequences with similar shapes but potentially different temporal shifts, such as... Figure 5 As shown.

[0178] The goal here is to use the Time Series K-Means algorithm to group provinces with similar development trajectories and patterns into the same cluster, and to reveal the diversity and differences in regional development in my country by analyzing the characteristics of each cluster.

[0179] The core features used for clustering are attributes across six dimensions: agri, manu, build, trans, it, and edu. This invention focuses on the trends of these key industry indicators over time.

[0180] This invention is implemented using Time Series K-Means provided by the tslearn library. The key configuration includes the Time Series K-Means clustering algorithm. The distance metric is set to dtw, explicitly using dynamic time warping distance. This enables the clustering to identify morphologically similar but potentially time-distorted development patterns.

[0181] (5) Determine the optimal number of clusters

[0182] Choosing the appropriate number of clusters K is a crucial step in the K-Means algorithm. Two commonly used methods assist in the decision-making process: the elbow method, which calculates the sum of squared intra-cluster distances for different K values; and the inertia method, which measures the compactness of samples within a cluster and typically decreases as K increases. This invention seeks the K value corresponding to the "elbow" where the rate of inertia decreases from steep to gradual, considered a superior choice. For each K value, the silhouette coefficient is calculated as the average silhouette coefficient of the clustering results. The silhouette coefficient comprehensively considers the compactness of intra-cluster samples and the separation of samples between clusters, with a value ranging from -1 to 1. A score closer to 1 indicates better clustering performance. This invention selects the K value that maximizes the average silhouette coefficient as the candidate. When calculating the silhouette coefficient, this invention utilizes a pre-calculated DTW distance matrix between all provincial time series and passes it to the silhouettescore function to ensure consistency between the evaluation metric and the distance metric used in the clustering algorithm.

[0183] This invention plotted curves of Inertia and the profile coefficient as a function of the K value, as shown in the figure. Figure 6 :

[0184] The left figure shows the elbow rule curve, illustrating how the sum of squared intra-cluster distances changes with increasing K value; the right figure shows the profile coefficient curve, illustrating how the average profile coefficient changes with increasing K value.

[0185] According to the elbow rule diagram, a clear elbow was observed at K=4, and the profile coefficient reached its maximum value of 0.630 when K=2. Therefore, K=4 was finally selected as the number of clusters for this clustering.

[0186] (6) Optimal clustering results

[0187] Using a determined optimal K value of 4, this invention performs the final DTW clustering on the multidimensional time series data of provinces. The algorithm's running parameters include setting n_init to max(1, N_INIT, KMEANS / / 2); max_iter to max(10, MAX_ITER, KMEANS / / 2); and random state to RANDOMSTATE to ensure the reproducibility of the results.

[0188] (7) Clustering visualization

[0189] To understand the development patterns behind each cluster, this invention provides a visual analysis of the clustering results.

[0190] Cluster Profile Visualization: To intuitively demonstrate the differences in key economic structural features among different clusters, this invention first calculates the time average of all provinces within each cluster on the selected feature (agri manu build trans it edu), and then plots grouped bar charts for comparison. For example... Figure 7 As shown.

[0191] The chart compares the average values ​​of the X-axis across six key features on the Y-axis for each cluster. Each cluster has a set of bars, with different colored bars representing different features. This figure illustrates the typical differences in industrial structure among different clusters. Figure 8 Take the average characteristics of the education industry as an example.

[0192] This static profile comparison provides a foundation for the subsequent analysis of the time-series dynamic evolution of cluster centers in this invention.

[0193] Next, this invention analyzes the time series plots of the cluster centers of each cluster on various features. The cluster centers represent the average or typical pattern of the time series of all provinces within that cluster. Figure 9 , Figure 10 It can be reflected.

[0194] Figure 9 , Figure 10 The time-series patterns of the cluster centers of four clusters across six features are shown. Each subplot represents a typical development trajectory of a cluster on a single feature.

[0195] By analyzing these center curves, this invention can summarize the dynamic evolution characteristics of each cluster in terms of robot density development.

[0196] Cluster 0: This cluster is characterized by an early start in robot density but slow initial growth. Its standardized value remained at a low level for a long period, but from approximately 2010 to 2012, it began to show a sustained and robust acceleration in growth, eventually reaching a density value significantly above the average level. This represents a steady catch-up development model.

[0197] Cluster 1: This cluster exhibited the longest period of low-level stagnation. Before 2015, its standardized robot density grew extremely slowly, remaining almost at its initial low level. Despite starting the latest, its growth momentum stabilized later, and by 2021 its standardized density value had reached a relatively high level, even slightly higher than Cluster 0. This can be seen as a late-blooming growth pattern.

[0198] Cluster 2: This cluster initially had a relatively high standardized robot density. Its growth trajectory was similar to other clusters in the early stages, but after 2017 it exhibited the most rapid acceleration, with a sharply increasing slope, ultimately reaching the highest standardized density level among all clusters. This represents a high-speed, leading development model.

[0199] Cluster 3: Similar to Cluster 0, this cluster experienced a later start and slower growth in robot density in its early stages. Growth accelerated around 2012, reaching a high peak in 2018-2019. This cluster exhibits a unique pattern, showing a significant slowdown or even decline in growth near 2020, followed by a slight recovery. This may reflect specific economic fluctuations, policy adjustments, or statistical changes, representing a fluctuating growth pattern. Figure 1 A pie chart showing the proportion of each cluster is provided.

[0200] Step 3: Establish a prediction model: Establish an XGBoost prediction model based on time series.

[0201] Building upon the aforementioned data preprocessing, anomaly detection, and cluster analysis, the core objective here is to utilize machine learning models to capture the complex time-series relationship between regional economic development and robot density, and thereby predict future development trends. By predicting future robot density, this invention can further assess the automation transformation status of specific industries in various regions and provide early warning information for potentially affected employees, thus providing data support for policy-making and labor market adjustments.

[0202] Considering the rich features and potential nonlinear relationships in the data, this invention chooses to use the XGBoost (Extreme Gradient Boosting) model, an efficient, flexible and accurate gradient boosting tree algorithm. It performs well in prediction tasks involving structured data and time-series related features. It can effectively handle high-dimensional features and capture complex interaction effects, and has a built-in regularization term to prevent overfitting. Figure 12 This reflects the basic logic of the xgboost model.

[0203] This section will detail the construction, training, and evaluation process of the XGBoost model, demonstrate how to use the trained model to predict future robot density, and finally explore how to apply the prediction results to automation transformation assessment and employment early warning.

[0204] 1. Features and target variables

[0205] Target variable: The prediction target of this invention is the robot installation density in each province.

[0206] Input Features: The input features used to train the XGBoost model include all numerical features generated in the first preprocessing stage. All features are used as input except for the identifier columns (Province COLUMN and Year COLUMN) and the original target variable (TARGET RENAME) itself. This includes all derived features generated from the original six industry indicators (agrimanu build trans it edu) through smoothing, lag, differencing, and rolling statistics calculations, particularly the lag and rolling statistics features, designed to help the model capture the dynamic dependencies of time series data.

[0207] Time series data partitioning

[0208] To simulate real-world prediction scenarios when evaluating model performance and to avoid data leakage, this invention strictly divides the dataset according to chronological order: The training set contains data from all provinces with years less than 2019 (less than the VALIDATION SPLIT YEAR), used for learning model parameters. The validation set contains data from all provinces with years equal to 2019 (equal to VALIDATION SPLITY YEAR), used for hyperparameter tuning and performance monitoring during model training. The test set contains data from all provinces with years greater than or equal to 2020 (greater than or equal to TESTSPLITY YEAR), used to evaluate the model's final generalization ability on unseen data after training. This division ensures that the model predicts future data based on past data, conforming to the fundamental principles of time series forecasting.

[0209] Model Configuration

[0210] The key hyperparameter settings for the model are as follows:

[0211] For an objective:reg:squarederror regression task, optimize the mean squared error.

[0212] eval metric:[rmse,mae] uses root mean square error and mean absolute error as evaluation metrics;

[0213] eta (learning rate): 0.02. Smaller learning rates typically require more trees but help improve model accuracy.

[0214] max depth: 5 limits the maximum depth of a single tree to prevent overfitting;

[0215] subsample:0.7 In each iteration, 70% of the training samples are randomly selected for tree construction;

[0216] colsample bytree:0.7 randomly extracts 70% of the features each time the tree is built;

[0217] Early stopping rounds:50 If the performance on the validation set does not improve for 50 consecutive rounds, training is stopped early to prevent overfitting and the best-performing model is selected.

[0218] 2. Model Training

[0219] The XGBRegressor interface of XGBoost is used to build the regression model. The model is trained using the training set (X train, Y train) and monitored for early stopping using the validation set (X val, Y val). During training, evaluation metrics (RMSE, MAE) on the validation set are output periodically. Finally, the model state that performs best on the validation set is selected. The trained model is saved to a file for later use.

[0220] 3. Model Prediction

[0221] Using a trained XGBoost model, this invention predicts robot density for the next year (taking 2022 as an example). To more intuitively illustrate these predictions, this invention generates a sorted bar chart, as shown below. Figure 13 As shown.

[0222] Model Validation

[0223] After the model was trained, its final performance was evaluated using a reserved test set of data from 2020 and later, which the model had never seen before. The following evaluation metrics were calculated on the test set. To gain a deeper understanding of the model's predictive behavior, the following visualization chart, a feature importance plot, was created to analyze the features that the model considers most important when making predictions.

[0224] Figure 14 This demonstrates the importance of some features.

[0225] The XGBoost model features are ranked in the top 30 for importance. The figure shows the features that contribute the most to the prediction of robot density, such as robot densitylag1, manu roll mean w3, etc. This helps to understand the key drivers affecting changes in robot density.

[0226] Predicted vs. True Value Plot: This plot compares the predicted and true values ​​on the test set using a scatter plot. For example... Figure 15 .

[0227] A scatter plot comparing the predicted values ​​(Y-axis) and the actual values ​​(X-axis) of the test set. The closer the point is to the red y=x diagonal, the more accurate the prediction.

[0228] Residual plot: This plot shows the relationship between the predicted residuals and the predicted values. For example... Figure 16 .

[0229] A scatter plot of test set residuals (Y-axis) versus predicted values ​​(X-axis) is ideal. The residuals should be randomly distributed near the zero line (red dashed line) without any obvious pattern. This can be used to check whether the model has heteroscedasticity or other systematic problems.

[0230] Step 4: Establish an early warning framework and conclusions based on the predictive model.

[0231] This study correlates predicted robot density growth with employment trends in related industries. It identifies regions with rapid automation and potential pressure on employment in related sectors, issuing early warnings to enable relevant departments to develop timely response strategies.

[0232] This framework provides a way to translate model predictions into support for practical management decisions. Specific evaluation criteria and early warning trigger conditions need to be further refined and quantified based on actual policy needs and more comprehensive data.

[0233] The XGBoost model was successfully used to predict robot density in various provinces. The model demonstrated good predictive performance (88.34% accuracy) on the test set. The forecasting process for the next year was demonstrated, and significant regional differentiation trends were observed based on the forecast results. Importantly, this invention proposes an analytical framework for assessing automation transformation and providing early warning of employment based on the forecast results. This aims to apply the predictive model to practical problems, providing a data-driven perspective for understanding the potential impact of automation on regional economies and employment. Although the forecast itself relies on simplifying assumptions, and the early warning framework needs further refinement, this invention demonstrates a complete technical path from data analysis to predictive early warning, which has certain reference value.

[0234] In summary, the predictive model of this invention can deeply analyze the spatiotemporal dynamics, regional patterns and future trends of robot adoption at the provincial level in China, and assess its potential employment impact. It has a high accuracy rate, with a measured accuracy rate of 88.34%.

[0235] The descriptions and practices disclosed in this invention are readily apparent and understandable to those skilled in the art, and various modifications and refinements can be made without departing from the principles of this invention. Therefore, any modifications or improvements made without departing from the spirit of this invention should also be considered within the scope of protection of this invention.

Claims

1. A predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model, characterized in that, Includes the following steps: Step 1, Data Preprocessing: A high-dimensional feature space is constructed through k-NN interpolation, RobustScaler normalization, and comprehensive temporal feature engineering; Step 2, Data Validation and Clustering: PCA anomaly detection and DTW time series K-Means clustering are used to ensure data quality and understanding of heterogeneity; Step 3: Establish a prediction model: Establish an XGBoost prediction model based on time series data; Step 4: Establish an early warning framework and conclusions based on the prediction model.

2. The predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model according to claim 1, characterized in that, Step 1 includes the following methods: Step 1.1: Data Preprocessing To ensure that each region in the data table has complete annual data records within the time span, the k-nearest neighbor (k-NN) imputation method is adopted. This method is based on the idea of ​​clustering similar things, and finds the k nearest neighbor samples that are most similar to the sample containing missing values ​​in known features, and uses the information of these neighbor samples to estimate the missing values. (1) The distance metric uses cosine similarity to measure samples, assuming it is the feature vector x of sample i and sample j. i and x j The degree of similarity between vectors, with cosine similarity focusing on vector direction rather than absolute magnitude, is suitable for handling data with features of different magnitudes. Similarity(x i ,x j The cosine similarity (x) between sample i and sample j measures the degree of similarity between the feature vectors of the two samples. i x j Let x be the feature vector of sample i and sample j respectively. i .x j It is the dot product of the feature vectors of sample i and sample j, ||x i ||,||x j || L2 norm x of the eigenvectors of sample i and sample j ik and x jk The values ​​of sample i and sample j on the k-th feature, where n is the total number of features of the sample, i.e. the dimension of the feature vector; (2) Nearest neighbor selection: After calculating the similarity between samples, select the k samples that are closest to the target sample, i.e., the samples with the highest similarity, as the nearest neighbor set N. K After testing and evaluation, k=7 was chosen as the optimal number of nearest neighbors. N k (x new )={(x j ,y j )|j∈I k } N k (x new ) represents the target sample x with missing values. new The set of the k nearest neighbors, x new The target sample with missing values, x j The feature vector y of the j-th sample in the nearest neighbor set j The variable value of the j-th sample in the nearest neighbor set can be understood as the feature information of that sample, I. k An index set containing these k nearest neighbor samples is used to mark which samples are the nearest neighbors of the target sample; (3) Imputation strategy for missing values ​​in numerical data The filling is performed using the arithmetic mean of the feature values ​​corresponding to its k nearest neighbors, calculated as follows: k represents the estimated value of the missing features in the target sample, i.e., the final result used to fill in the missing values; k is the number of nearest neighbor samples selected; N is the number of nearest neighbor samples selected. k The set of the k nearest neighbors of the target sample, x' j The value of the j-th sample in the nearest neighbor set for the missing feature; By implementing k-NN interpolation, gaps in the dataset are effectively filled, the integrity of the time series data is maintained, and a data foundation is provided for subsequent analysis. Step 1.2, Data Normalization Data normalization is performed using the RobustScaler method. RobustScaler is a normalization technique that is insensitive to outliers. It scales the data by subtracting the median of the features and then dividing by the interquartile range (IQR). The calculation formula is as follows: X represents the original feature data to be processed. scaled After scaling the feature data by RobustScaler normalization, Median(X) is the median of all sample values ​​for the feature, and IQR(X) represents the interquartile range of the feature. The interquartile range (IQR) is defined as follows: IQR(X) = Q3(X) - Q1(X) Q3(X) represents the upper quartile of this feature; Q4(X) represents the lower quartile of this feature; Derivative feature construction and feature engineering The feature engineering in this phase includes the following aspects. All operations are performed in groups based on data sorted by province and year to ensure the correctness of time series operations. (1) Smoothing A simple moving average (SMA) was applied to the target variable and all six selected industry characteristic variables. The time window size (SMOOTHING_WINDOW) was set to 3 periods. During the calculation, it was allowed that there were fewer than 3 data points in the early stage of the window (min_periods=1) to avoid generating too many missing values ​​at the beginning of the sequence. SMA t Let X be the simple moving average at time t, W be the size of the moving window, and X be the number of historical data points used in the calculation of the moving average. t-i Let t be the value of variable X at time ti; (2) Lag characteristics For the target variable robot density and all six industry characteristic variables, lag values ​​for the past one period and the past two periods are created respectively, with the maximum lag order N_LAGS set to 2; (3) Difference characteristics The first-order difference of the target variable robot density was calculated, which is the value of the current time point minus the value of the previous time point. The difference order is 1, and the robot_density_diff1 feature column is generated. According to the code implementation, this processing only performs difference on the target variable. ΔX t =X t -X t-1 ΔX t： The backward difference of variable X at time t describes the magnitude of the change in variable X at time t relative to time t-1. t Let X be the value of variable X at time t. t-1 Let X be the value of variable X at time t-1; (4) Window statistical characteristics Capture the local dynamic characteristics of variables over a recent period, such as central tendency, volatility, and range of values; set the rolling window size to three periods, allowing for fewer than three data points in the initial stage of the window; calculate the following statistics: mean, standard deviation, minimum, and maximum. σ roll,t μ represents the rolling standard deviation at time t, used to measure the degree of fluctuation of a variable within a rolling window ending at time t. roll,t This represents the average value of variable X within the scrolling window; (5) Handling missing values ​​introduced by derived features The aforementioned lag, differencing, and windowing operations naturally generate missing NaN values ​​at the beginning of each province's time series. To ensure data quality for subsequent analysis and modeling, all rows containing any NaN values ​​are removed from the dataset.

3. The predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model according to claim 1, characterized in that, Step 2 includes the following methods: Data verification (1) PCA Algorithm Principal component analysis (PCA) is used to perform a holistic scan on the preprocessed dataset containing rich derived features, identify potentially abnormal province-year data points, evaluate the overall usability of the data based on the detection results, and provide a basis for subsequent clustering and predictive modeling. (2) Anomaly detection based on reconstruction error The reconstruction error method is used to identify outliers. The reconstruction error is calculated as follows: for each data point, the mean square reconstruction error is calculated on all standardized features. MSE i Let be the mean squared reconstruction error of the i-th data point, used to measure the degree of deviation of this data point after PCA reconstruction, p represents the total number of features, j is the index variable for summation, used to traverse each feature dimension, and x scaled,ij Let x represent the standardized value of the i-th data point on the j-th feature. reconstructed,ij This represents the PCA reconstruction value of the i-th data point on the j-th feature; Determine the threshold: Select all data points MSE i The 99th percentile of the value is used as the threshold for judging outliers; outlier identification: if the MSE of a data point is... i If the value is greater than the threshold, it is marked as a potential outlier. (3) Visualization of abnormal data To more intuitively understand the distribution of outliers and the setting of thresholds, visual charts are created; (4) Construction of K-means model The K-Means algorithm for time series data was chosen. This algorithm is used to cluster time series data, with the goal of finding K cluster centers. The time series data for each province is X... i Assigned to the only cluster C K In this case, all sample points within a cluster are such that they are all within the cluster center M. K Minimize the sum of squares of the distances; Dynamic Time Warping (DTW) is chosen as the distance metric for time series K-Means. DTW works by finding the distance between two time series X = (x... 1, x2,K x n ) and Y = (y1, y 2, K y m The optimal normalized path between points () is calculated, and their similarity is determined. This path defines a point-to-point mapping relationship that minimizes the cumulative distance between corresponding points. The DTW distance is usually calculated using dynamic programming, and its core recursive relationship is as follows: D(i,j)=cost(x i ,y j )+min\begincasesD(i-1,j) D(i,j-1) D(i-1,j-1)\endcases Where cost(x) i ,y i () is point x i and y i The distance between them, D(i,j) is the cumulative DTW distance between subsequences X[1...i] and Y=[1...j], and the final d DTW (x,y) is equivalent to D(n,m). DTW allows for non-linear scaling of the time axis and can identify sequences with similar shapes but may be offset in time. (5) Determine the optimal number of clusters Two methods are used to assist decision-making: the elbow rule calculates the sum of squared intra-cluster distances for different K values, and Inertia measures the compactness of intra-cluster samples, which decreases as K increases; for each K value, the average silhouette coefficient of the clustering results is calculated, and the K value that makes the average silhouette coefficient the highest is selected as the candidate. When calculating the silhouette coefficient, the DTW distance matrix between the time series of all provinces is pre-calculated and passed to the silhouette score function to ensure that the evaluation index is consistent with the distance metric used by the clustering algorithm. (6) Optimal clustering results Using the determined optimal K value of 4, the final DTW clustering is performed on the multidimensional time series data of the province. The algorithm running parameters include setting n_init to max(1, N_INIT, KMEANS / / 2); max_iter to max(10, MAX_ITER, KMEANS / / 2); and random state to RANDOM STATE to ensure that the results are reproducible. (7) Clustering visualization Cluster profile visualization: In order to intuitively show the differences in key economic structure features among different clusters, we first calculate the time average of all provinces in each cluster on the selected features, and then draw grouped bar charts for comparison.

4. The predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model according to claim 1, characterized in that, Step 3 includes the following methods: (1) Features and target variables Target variable: The prediction target is the robot installation density in each province. Input features: The input features used to train the XGBoost model include all numerical features generated in the first part of the preprocessing stage. Except for the identifier column and the original target variable itself, all features are used as input. This includes all derived features generated from the original 6 industry indicators through smoothing, lag, differencing and rolling statistics, especially the lag term and rolling statistics features, which are designed to help the model capture the dynamic dependencies of time series. Time series data partitioning The dataset is divided strictly according to chronological order: training set, validation set, and test set. Model Configuration (2) Model training The XGBRegressor interface of XGBoost is used to build a regression model. The model is trained using the training set (X train, Y train) and monitored for early stopping using the validation set (X val, Y val). During training, the evaluation metrics (RMSE, MAE) on the validation set are output periodically. Finally, the model state that performs best on the validation set is selected. The trained model is saved to a file for later use. (3) Model prediction Using a trained XGBoost model, the robot density for the next year is predicted. Model Validation After the model training is completed, its final performance is evaluated using reserved test set data from 2020 and later, which the model has never seen before.

5. The predictive modeling method for the impact of industrial robot introduction on employment based on the XGBoost model according to claim 1, characterized in that, Step 4 includes the following methods: By correlating the predicted growth in robot density with the changing trends in employment in related industries, we can identify regions where automation is progressing rapidly and employment in related industries may be under pressure, and issue early warning signals so that relevant departments can formulate response strategies in advance.

Images