Prediction of cold chain loading environment for agricultural products based on K-medoids-LSTM-XGBoost ensemble model

Data fusion using K-medoids algorithm

The K-medoids algorithm is an improvement upon the K-means algorithm. The K-means algorithm is susceptible to the influence of noise and outlier data, which can significantly distort the data distribution if extreme values or outliers are present. In contrast, the K-medoids algorithm employed in this experiment addresses this sensitivity issue by mitigating the impact of isolated points and noisy data on cluster assignment. It exhibits greater robustness compared to the K-means algorithm. The primary distinction between these two algorithms lies in their approach to determining cluster centroids. The K-means algorithm calculates the average value of the current cluster as the centroid, whereas the K-medoids algorithm selects a specific point within the cluster to serve as the centroid.

The steps involved in the K-medoids algorithm are as follows:

• 1)

  Randomly select k data points from the sample set as initial medoids.

• 2)

  Compute the distances from non-medoid data points in the sample set to the medoids.

• 3)

  Assign non-medoid data points to the cluster with the minimum distance to the medoid.

• 4)

  Randomly select a non-medoid data point $O_r$, calculate the sum of distances between this data point and all other data points D, if the calculated D is smaller than the original D, replace the original medoid with $O_r$.

• 5)

  Repeat steps (2), (3), and (4) until the difference between the calculated D of two adjacent iterations remains unchanged and each data point in the sample set has been traversed, indicating that the medoids no longer change; otherwise, continue with step 4.

• 6)

  Algorithm terminates, and the final medoids are outputted.

Long short-term memory neural networks

Long short-term memory neural networks is a specialized type of RNN that possesses unique storage and forgetfulness functions, enabling it to address the drawback of limited short-term memory in traditional RNNs and learn long-term dependencies.

LSTM consists of three gates: the forget gate ( $f_t$), the input gate ( $i_t$), and the output gate ( $O_t$). These gates control the cell state within LSTM by selectively incorporating or discarding information. The biases for the forget gate, input gate, and output gate are represented as $b_f$, $b_i$ and $b_o$, respectively. Furthermore, $b_C$ denotes the bias for the candidate vector. The weights associated with the forget gate, input gate, output gate, and candidate vector are denoted as $W_f$, $W_i$, $W_o$ and $W_C$, respectively. The activation functions used are the sigmoid function (σ) and the hyperbolic tangent function (tanh). $x_t$ represents the input at time step t, while $h_t$ and $h_{t-1}$ represent the outputs of the model at time steps t and t − 1, respectively. $C_t$ and ${\tilde{\mathrm{C}}}_{\mathrm{t}}$ correspond to the candidate vector at time step t and the updated value of the candidate vector at time step t, respectively. Figure 3 illustrates the internal structure diagram of a temporal LSTM network.

(1) Forget gate

(1) $$\begin{array}{c}{f}_t=\sigma \left(W_f\cdot \left[h_{t-1},x_t\right]+b_f\right)\end{array}$$

(2) Input gate

(2) $$\begin{array}{c}{i}_t=\sigma \left(W_i\cdot \left[h_{t-1},x_t\right]+b_i\right)\end{array}$$

(3) $$\begin{array}{c}{\tilde{C}}_t=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(W_i\cdot \left[h_{t-1},x_t\right]+b_i\right)\end{array}$$

(3) Unit

(4) $$\begin{array}{c}{C}_t=f_t\ast C_{t-1}+i_t\ast {\tilde{\mathrm{C}}}_{\mathrm{t}}\end{array}$$

(4) Output gate

(5) $$\begin{array}{c}{O}_t=\sigma \left(W_o\cdot \left[h_{t-1},x_t\right]+b_o\right)\end{array}$$

(6) $$\begin{array}{c}{h}_t=O_t\tanh \left(C_t\right)\end{array}$$

XGBoost algorithm

XGBoost, short for eXtreme Gradient Boosting, is an algorithm developed by Professor Tianqi Chen as an enhancement to the Gradient Boosting algorithm. Both XGBoost and gradient boosting decision trees (GBDT) utilize the concept of boosting, but their main difference lies in the definition of the objective function. While GBDT optimizes the objective function by expanding it to the first order using a Taylor series approximation, XGBoost expands the objective function to the second order while also incorporating additional regularization terms. This approach allows XGBoost to retain more relevant information from the objective function while controlling model complexity to prevent overfitting.

The XGBoost algorithm begins by generating a tree model at each iteration t. The tree model to be trained in the t^th iteration is denoted as $f_t\left(x\right)$.

(7) $${\hat{y}}_i^{\left(t\right)}=\sum _{k=1}^t{f}_k\left(x_i\right)={\hat{y}}_i^{\left(t-1\right)}+f_t\left(x_i\right)$$

${\hat{y}}_i^{\left(t\right)}$ represents the predicted result of sample i after t iterations, while ${\hat{y}}_i^{\left(t-1\right)}$ represents the predicted result from the previous t − 1 trees. $f_t\left(x_i\right)$ denotes the tree model of the t^th tree. The objective function consists of two components: the error function and the regularization term. The objective function employed in this model and its associated regularization term are as follows:

(8) $$O{B}_{j^{\left(t\right)}}=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i^{\left(t\right)})+\sum _{t=1}^t\mathrm{\Omega }\left(f_t\right)=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i^{\left(t-1\right)}+f_t\left(x_i\right))+\mathrm{\Omega }\left(f_t\right)+\sum _{t=1}^{t-1}\mathrm{\Omega }\left(f_t\right)$$

(9) $$\mathrm{\Omega }\left(f_t\right)=\gamma T+{\displaystyle \frac{1}{2}\lambda \sum _{\mathrm{\Omega }-1}^T{w}_o^2}$$

$O{B}_{j^{\left(t\right)}}$ represents the objective function when constructing the t^th tree. $l(y_i,{\hat{y}}_i^{\left(t\right)})$ denotes the loss function. ${\hat{y}}_i^{\left(t\right)}$ represents the predicted result of sample i after t iterations, while ${\hat{y}}_i^{\left(t-1\right)}$ represents the predicted result from the previous t − 1 trees. $\mathrm{\Omega }\left(f_t\right)$ represents the regularization term of the t^th tree, where a smaller value indicates a lower complexity of the model, resulting in stronger generalization ability and more stable predictions.

The complexity of the tree model consists of two components: the number of leaf nodes T and the scores $w_o$ assigned to each leaf node. γ and λ serve as coefficients for the regularization term, controlling the number of leaf nodes and their scores.

The objective function is optimized using Taylor expansion, resulting in the following expression:

(10) $$O{B}_{j^{\left(t\right)}}\approx \sum _{i=1}^n[l(y_i,{\hat{y}}_i^{\left(t\right)})+g_i{f}_t\left(x_i\right)+{\displaystyle \frac{1}{2}{h}_i{f}_t^2\left(x_i\right)]+\mathrm{\Omega }\left(f_t\right)+constant}$$

Definition:

(11) $$g_i={\partial }_{\hat{y}\left(t-1\right)}l(y_i,{\hat{y}}_i^{\left(t-1\right)})$$

(12) $$h_i={\mathrm{\partial }}_{\hat{y}\left(t-1\right)}^{2}l(y_i,{\hat{y}}_i^{\left(t-1\right)})$$

The XGBoost algorithm utilizes a greedy algorithm to enumerate all possible tree structures and find the optimal structure that is subsequently added to the model. To prevent infinite tree generation, it is necessary to set stopping conditions. Common conditions include:

• 1)

  Setting a hyperparameter, max_depth, which stops tree generation when the tree reaches the maximum allowed depth to prevent overfitting.

• 2)

  Stopping the iteration and tree generation when the gain from a potential split falls below a specific threshold. At this point, the split can be disregarded. The threshold parameter is typically a coefficient of the number of leaf nodes T in the regularization term.

Parameter optimization in particle swarm optimization

The PSO algorithm is a stochastic global optimization technique that utilizes interactions among particles to discover the optimal regions within a complex search space. The algorithm draws inspiration from the migration and collective behavior observed in bird flocking during foraging. In the PSO algorithm, birds are abstracted as massless, dimensionless particles, and a group of particles represents a bird flock. The process of searching for the optimal solution within a defined area corresponds to the foraging behavior of the bird flock, achieved through systematic random movements of the particles. During the search for the optimal solution, particles collaborate and share information, with each particle following the vicinity of the current best particle to search for the optimal solution. In this process, particles continuously compare their own information with that of the best solution, enabling them to maintain an optimal state.

The PSO algorithm is described as follows:

In a D-dimensional region, there exist m particles forming a particle swarm { $x_1,x_2,x_3,\dots ,x_m$}. Each particle is characterized by its velocity $x_i=\left\{x_{i1},x_{i2},\dots ,x_{iD}\right\}$ and position $v_i=\left\{v_{i1},v_{i2},\dots ,v_{iD}\right\}$. The individual best position found by the i^th particle is denoted as $P_i=\{P_{i1},P_{i2},P_{i3}\dots ,P_{iD}\}$, while the overall best position discovered by the entire particle swarm is represented as $P_{gBest}=\left\{P_{gBest1},P_{gBest2},P_{gBest3}\dots ,P_{gBestD}\right\}$.

During the t^th iteration, each i^th particle updates its velocity and position using the following formulas, as shown in Eqs. (13) and (14):

(13) $$v_{id}^{t+1}=\omega v_{id}^t+c_1{r}_1\left(P_{id}-x_{id}^t\right)+c_2{r}_2\left(P_{gBestd}-x_{id}^t\right)$$

(14) $$x_{id}^{t+1}=x_{id}^t+v_{id}^{t+1}$$

In the above equations, i ranges from 1 to m, and d ranges from 1 to D. The inertia weight ω represents the degree to which the particle retains memory of its previous velocity. A higher ω value enhances the particle’s global search capability, while a lower ω value emphasizes its local search ability. The learning factors, $c_1$ and $c_2$, correspond to the particle’s self-cognition ability and social sharing ability, respectively. These factors assist particles in converging towards local optimal points. The variables $r_1$ and $r_2$ denote random numbers. The algorithm terminates either when the maximum iteration count is reached or when the desired precision is achieved.

The flowchart of the PSO algorithm is illustrated in Fig. 4.

A combined predictive model: K-medoids-LSTM-XGBoost

The ensemble prediction model integrates information from multiple individual prediction models and combines their results using specific methods. By extracting and integrating this information, the ensemble model leverages the strengths of different models, thereby effectively enhancing the accuracy of prediction results and bolstering the persuasiveness of experimental outcomes.

In this experiment, a weighted approach is employed to combine the two models, as demonstrated by the following formula:

(15) $$f_t=w_1{f}_{1t}+w_2{f}_{2t},t=1,2,3,\dots ,n$$

Here, $f_{1t}$ denotes the prediction results derived from the LSTM model, while $f_{2t}$ represents the prediction results obtained from the XGBoost model.

The weight coefficients, denoted as $w_i$, are determined according to the following equations:

(16) $$w_1=\{\begin{array}{c}{\displaystyle \frac{{\alpha }_1}{{\alpha }_1+{\alpha }_2},{\alpha }_1\ge {\alpha }_2}\\ {\displaystyle \frac{{\alpha }_2}{{\alpha }_1+{\alpha }_2},{\alpha }_1<{\alpha }_2}\end{array}$$

(17) $$w_2=\{\begin{array}{c}{\displaystyle \frac{{\alpha }_1}{{\alpha }_1+{\alpha }_2},{\alpha }_1\ge {\alpha }_2}\\ {\displaystyle \frac{{\alpha }_2}{{\alpha }_1+{\alpha }_2},{\alpha }_1<{\alpha }_2}\end{array}$$

The values α₁ and α₂ represent the prediction errors of the LSTM and XGBoost models, respectively. The model with smaller error and higher accuracy will have a relatively higher proportion in the ensemble model, whereas a model with larger error will have a lower proportion.

The specific steps are as follows:

• 1)

  Collect temperature data from sensors in the remote cold chain vehicle.

• 2)

  Preprocess the data, including cleaning the data, removing outliers, filling missing data, and normalizing the data.

• 3)

  Apply the K-medoids algorithm to fuse the data and improve its accuracy.

• 4)

  Divide the fused data into training and testing sets in proportion.

• 5)

  Initialize the LSTM neural network model and optimize its parameters using the PSO algorithm. Then, input the testing set data into the trained LSTM model for prediction, obtaining the predicted data.

• 6)

  Initialize the XGBoost model and import the training set data. Perform grid search to obtain the optimal parameters for the XGBoost model. Then, input the testing set data into the trained XGBoost model for prediction, obtaining the predicted data.

• 7)

  Combine the two sets of models using the weighted fusion method described above to generate the final prediction results.

• 8)

  Introduce evaluation metrics for model assessment and conduct validation and analysis of the prediction results for the two individual models and the combined model. Additionally, to enhance the persuasiveness of the results, this study will compare, validate, and analyze them using models such as random forest and gate recurrent unit (GRU).

• 9)

  Figure 5 illustrates the flowchart depicting the prediction process of the cold chain loading environment.

Experimental preparation

K-medoids data fusion results

Figures 6 and 7 present the original data for ambient temperatures of 23 °C and 26 °C, respectively. These figures indicate periodic temperature variations with a period exceeding 190 to 200 min.

Unlike the K-means algorithm, the K-medoids algorithm does not require the computation of means as centroids. Instead, it selects existing points as medoids during each iteration. This characteristic enables the K-medoids algorithm to perform computations rapidly and efficiently. By employing this method for data fusion on the collected temperature time series, it becomes possible to effectively extract the desired feature information, thereby enhancing the accuracy of the prediction results. Figure 8 illustrates the fusion results of the original data.

LSTM neural network prediction results

After employing the K-medoids algorithm for data fusion, the fused data is subsequently subjected to normalization and divided into distinct training and testing sets, following a predetermined ratio. In this experiment, to showcase the universality of the model, the conventional approach to data set partitioning is eschewed. Instead, the training set encompasses data recorded at temperatures of 23 °C and 26 °C, while the testing set comprises data recorded at 28 °C. Figure 9 presents the original data recorded at 28 °C.

Based on Fig. 10, the hidden layer predominantly comprises a recurrent neural network (RNN) constructed using LSTM and Dropout. The intermediate RepeatVector layer serves to establish a connection between the encoder and decoder components. To mitigate convergence challenges associated with deep neural networks, the hidden layer employs a modest configuration of only two layers. The LSTM layer performs computations and predictions, whereas the Dropout layer safeguards against overfitting, thereby enhancing the LSTM model’s generalization capabilities.

The output layer will utilize the TimeDistributed layer to process the outputs from the LSTM hidden layer, reducing dimensionality and predicting the final outcome.

Once the LSTM model is constructed, the PSO algorithm is employed to optimize its parameters, aiming to enhance the model’s accuracy. The results of the LSTM model are depicted in Fig. 11.

XGBoost model predictions

Before running the XGBoost model, it is crucial to determine three types of parameters: general parameters, booster parameters, and task parameters. General parameters control the overall functions and are commonly used for tree models or linear models. Booster parameters can control the booster type (tree/regression) at each step. Once the general parameters are set, the corresponding booster parameter types are determined. Booster parameters can influence the model’s effectiveness and computational cost. In essence, tuning the XGBoost model largely involves adjusting the booster parameters. Task parameters control the performance of the training objective, such as classification or regression tasks, including binary or multi-class classification.

These three sets of parameters significantly impact the algorithm’s performance. After determining the general and task parameters, the next step involves parameter tuning, specifically the booster parameters. Common booster parameters include learning rate, minimum loss reduction for splitting (gamma), and maximum depth of trees. As the prediction results are directly influenced by the maximum depth of trees, optimizing this parameter is the first step in model tuning. The optimization process begins by assigning initial values to the parameters other than the maximum depth of trees, typically using common or default values. Through multiple iterations and adjustments, the optimal combination of parameters is obtained.

Once the maximum depth of trees is determined, a traversal search method is used to find the best combination of tree depth with other parameters. After the traversal search, the optimal parameter combination obtained is presented in Table 2.

After acquiring the optimal parameters, the dataset is inputted into the XGBoost model for prediction. The temperature prediction values and actual values derived from the XGBoost model demonstrate their respective trends, as depicted in Fig. 12.

Prediction results of the combined K-medoids-LSTM-XGBoost model

By employing the variable weight combination method described above, the LSTM model undergoes fine-tuning using the PSO algorithm, while the XGBoost model is tuned through an exhaustive search approach. Following this, the two sets of models are combined to form an ensemble prediction model. Finally, an analysis and prediction are conducted under identical prediction conditions.

Figure 13 illustrates the prediction results obtained from the K-medoids-LSTM-XGBoost ensemble forecasting model.

Figure 13 showcases the prediction performance of the K-medoids-LSTM-XGBoost ensemble forecasting model, revealing a favorable overall performance. While certain prediction errors are evident in specific aspects, such as peaks and valleys, these inconsistencies could be attributed to the limited dataset quantity and the relatively large prediction interval. Such factors may lead to insufficient training. In practical applications, as more data is continuously collected from each farm’s cold chain loading, the predictive accuracy of the model is expected to improve significantly.

Comparison and analysis

Model evaluation metrics provide a intuitive way to assess the accuracy of prediction results and serve as a means to evaluate the quality of predictions. Therefore, in order to evaluate the accuracy of the model’s predictions, this study will employ three metrics to evaluate the predicted results $y(i)$ against the true values of the dataset.

R-squared.

(19) $$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

Mean absolute error (MAE).

(20) $$MAE={\displaystyle \frac{1}{n}\sum _{i=1}^n\left|{\hat{y}}_i-y_i\right|}$$

Root mean square error (RMSE).

(21) $$RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}}$$

Comparison of predictive results for baseline models

In this experiment, the raw data underwent data fusion using the K-medoids algorithm. The resulting processed data was then utilized as input for comparative models in the prediction process. To facilitate meaningful comparisons, a variant of the RNN known as the GRU algorithm was employed. The GRU algorithm shares many similarities with RNN and provides a highly comparable approach. Additionally, random forests, an ensemble learning algorithm based on decision trees, was chosen as the second comparative model. This choice allows for a reasonable basis of comparison. The prediction results obtained from the GRU algorithm and Random Forests are displayed in the figure below.

From Figs. 11, 12, 14A, and 14B, it can be observed that the prediction outcomes of the GRU and Random Forest models exhibit a similar overall trend as the true curve. However, they lack accuracy in capturing the finer details compared to the LSTM and XGBoost algorithms. Specifically, discrepancies are evident in identifying peaks and valleys. Therefore, in the context of agricultural cold chain loading, the LSTM and XGBoost algorithms demonstrate relative advantages in single-model prediction. Consequently, this study selects the LSTM and XGBoost models as part of the ensemble model.

Comparison of predictive results for ensemble models

The predicted results of the K-medoids-GRU-Random Forest ensemble prediction model are shown in Fig. 15A.

Figures 15A and 15B clearly demonstrate that the K-medoids-LSTM-XGBoost ensemble prediction model significantly outperforms the K-medoids-GRU-Random Forest ensemble prediction model in terms of accuracy. This superiority is particularly evident in the fitting of the slopes on both sides of the peaks and valleys, where the K-medoids-LSTM-XGBoost ensemble prediction model exhibits a high degree of overlap. This finding strongly supports the superiority of the K-medoids-LSTM-XGBoost ensemble prediction model in the context of agricultural cold chain loading.

Practical applications

Based on the comparison between the predicted figures and the model error analysis presented in Table 3, it is evident that the prediction models proposed in this study successfully achieve temperature forecasting within the cold chain loading environment to varying degrees. Upon evaluating the predictive performance of these models, most of them exhibit a substantial overlap between the predicted values and the true value curve in their respective prediction figures. Notably, each combination model consistently outperforms the individual models. Among the ensemble models, the K-medoids-LSTM-XGBoost ensemble prediction model consistently outperforms the K-medoids-GRU-Random Forest ensemble prediction model in terms of MAE, RMAE, and R-square. This study also incorporates the K-medoids-Attention-BiGRU model for comparison. The model integrates an Attention mechanism and bidirectional gated recurrent units (BiGRU). The inclusion of the Attention mechanism allows the model to dynamically assign weights based on the importance of the input data, thereby enabling more accurate identification of key information within the sequence. BiGRU, with its superior ability to capture both forward and backward dependencies in time series, demonstrates improved performance compared to LSTM networks. As a result, this model is particularly well-suited for long-term time series predictions that require capturing both global and local significance, as well as handling complex data dependencies. However, as shown in Table 3, while the K-medoids-Attention-BiGRU model performs effectively, the K-medoids-LSTM-XGBoost model is more suitable for predicting the temperature environment in agricultural cold chain logistics. The experimental results demonstrate that the K-medoids-LSTM-XGBoost prediction model for temperature in the agricultural cold chain loading environment achieves an RMSE, MAE, and R-square of 2.5343, 5.1906, and 0.8971, respectively, indicating its remarkable accuracy in reflecting the actual cold chain loading environment data.

Based on the predictive results, this experiment preemptively controls the temperature environment within cold chain vehicles in practical applications. Through comparative analysis, it is concluded that the algorithm and methods proposed in this article can effectively predict temperature fluctuations inside cold chain vehicles after activating the refrigeration system under various ambient temperatures during agricultural product loading. This allows for an optimal determination of when to initiate the refrigeration system, thereby minimizing the adverse effects of excessively high temperatures on transported goods during loading, reducing energy consumption, and lowering transportation costs. In practical applications, temperature control and real-time monitoring are the core tasks of the entire experiment. To prevent uncontrollable incidents, users can establish a redundancy system to ensure that temperature is continuously monitored by sensors, even in the event of partial equipment failure. To ensure the system’s reliability, regular maintenance, upgrades, and component replacement are necessary. During data transmission, issues such as network congestion or equipment performance problems may lead to transmission delays or even data loss, ultimately affecting the system’s real-time prediction and cargo safety. Therefore, optimizing the data transmission protocol by reducing packet size and increasing transmission speed can help mitigate the risk of transmission failures.