Hopfield neural network based on clustering algorithms for solving green vehicle routing problem

As a result of the rapidly increasing distribution network, the toxic gases emitted by the vehicles to the environment have also increased, thus posing a threat to health. This study deals with the problem of determining green vehicle routes aiming to minimize CO 2 emissions to meet customers' demand in a supermarket chain that distributes fresh and dried products. A new method based on clustering algorithms and Hopfield Neural Network is proposed to solve the problem. We first divide the large-size green vehicle routing problem into clusters using the K-Means and K-Medoids algorithms, and then the routing problem for each cluster is found using the Hopfield Neural Network, which minimizes CO 2 emissions. A real-life example is carried out to illustrate the performance and applicability of the proposed method. The research concludes that the proposed approach produces very encroaching results.


Introduction
Since the beginning of time, logistics has been a critical link for the worldwide circulation of products. Due to the significant impact of the COVID-19 pandemic on the world economy, logistics has become an essential role in preventing supply chain disruptions. Companies and organizations need to make various decisions to optimize and effectively manage their processes in distribution systems. One of the most critical decisions in the distribution process is the vehicle routing decisions made from the facilities to the customers. The Vehicle Routing Problem (VRP), first scientifically defined by Dantzig and Ramser in 1959, is the problem of determining minimum cost routes with a fleet of vehicles to meet customer demands. The VRP, which is one of the problems in the NP-hard class, has many different types according to its constraints (Dantzig & Ramser, 1959). In recent years, studies that consider an environmental effect and minimize costs in distribution systems have become widespread, and the concept of green logistics has developed. The green vehicle routing problem (GVRP), on the other hand, can be defined as a VRP that takes energy consumption into account in general. Reducing fuel consumption in transport activities will significantly impact greenhouse gas emissions (Erdoğan & Miller-Hooks, 2012). The GVRPs are generally considered as two ways. The first is aimed to minimize fuel consumption, CO2 emissions, load amount, travel time in the general vehicle routing problem. The second is to minimize vehicle routes using either conventional and/or alternative fuels within a vehicle fleet.
There are different solution methods developed to solve GVRPs in the literature. The methods can be classified as an exact, a heuristic, and a metaheuristic. Since the VRP is in the NP-hard problem class, it is challenging and sometimes impossible to reach an optimal solution using exact methods if the size of the problem overgrows. For this reason, heuristic/meta-heuristic methods can be chosen to solve these problems, which usually have the potential to provide the best values (near-optimal) within reasonable solution time constraints (Lenstra & Rinnooy Kan, 1981). This study uses an artificial neural network (ANN) among the metaheuristic methods.
ANNs are complex systems consisting of the connection of artificial neurons, similar to neurons, which are the basic unit of the human brain, with different topology and network models. Instead of giving a step-by-step method in the program, the neural network generates its own internal rules that make the association and organizes them by comparing the results with examples (Rojas, 1996). In general, ANNs are examined in two main groups: feedforward and feedback. For the feedforward networks, neurons are organized through input, hidden, and output layers. The neurons in each layer are associated with neurons in the next layer. However, the layers have no connections between each other. The feedback networks are also called iterative networks. The most important feature of these neural networks is the feedback connection, parallelism, collective analog mode, and dynamic memory. Each neuron collects the inputs of hundreds and thousands of other neurons to provide cascading output (Jain et al., 1996;Lippmann, 1987). Hopfields Neural Network (HNN), one of the feedback networks, is among the first heuristic approaches developed to solve combinatorial optimization problems in ANNs. HNN was first proposed in 1982 by JJ Hopfield to solve optimization problems by reducing programming complexity and quickly reaching steady states. Its high computation speed provides the most suitable solution for neural networks to calculate the shortest path in a much shorter time. The Hopfield model is the energy minimizing network and is helpful as a content-addressable memory for solving combinatorial optimization problems (Hopfield, 1982). For these reasons, HNN is preferred in this study.
A few studies on HNN focus on solving the Traveling Salesman Problem (TSP), the most basic form of vehicle routing problems, can be found in the literature. Hopfield & Tank (1985) firstly used the HNN to solve the TSP. This technique minimizes an energy function with some terms and parameters.
Some studies focused on developing the energy function, ensuring that the model produces feasible solutions ( Van den Bout & Miller, 1988;Li et al., 2016). Selecting the initial parameter values is used more optimally (Kamgar-Parsi & Kamgar-Parsi, 1992;Lai & Coghill, 1992;Wu & Yang, 2004, Tan et al., 2005 to solve the TSP using HNN. Luo (2019) also improved Hopfield's energy function and proposed a genetic algorithm (GA) based on HNN for TSP. He used his proposed method on sample problems with 8 and 10 cities and compared the results with the classical HNN.
To obtain more effective solutions, some researchers hybridized HNN with some methods. Jolai & Ghanbari (2010) integrated Z-score and logarithmic data transformation techniques with HNN to solve the TSP. They compared their results on a sample problem considering ten cities. Oshima et al. (2015) proposed two methods to solve TSP by integrating the random cut and neuronal death method into HNN. Zhong et al. (2017) created a dynamic step-based continuous HNN and applied it on sample problems as 10 and 20 cities.
Abdel-Moetty (2010) applied HNN and the Ant colony algorithm (ACA) for the TSP using 15 cities and compared these results. Sarwar & Bhatti (2012) compared HNN and Lin and Kernighan heuristic for TSP by considering sample size as 5, 10, 20, and 50 customers. This paper proposes a hierarchical approach that consists of a clustering algorithm and HNN to solve the GVRP. We first divide the large-size GVRP into clusters using the K-Means and K-Medoids algorithms, and then vehicle routes are determined for each cluster using the HNN to minimize CO2 emissions. The proposed approach is also applied to a real-life example related to a distribution CVRP meeting the weekly demands of a supermarket chain. The HNN results compared with cluster first route second approach proposed by Comert et al. (2018) using Friedman test.
Studies in the literature show that HNN effectively solves small-size TSP problems. It confirms that the literature survey has given up to 50 customers. In order to effectively solve real-life problems with more than 50 customers, a large-size problem can be transformed into a small-size problem using clustering algorithms. Therefore, our study's first contribution is to deal with HNN's large-scale problems based on K-means and K-medoids clustering algorithms. Additionally, this brief survey shows a significant gap in the literature related to applying the HNN for VRPs. The GVRP is a more comprehensive type of TSP created by adding specific constraints to it. As a result, it requires more processing than the TSP as a solution is difficult. Considering the speed and computational power features of neural network approaches in applications, we contribute to the literature by using the HNN for the large size GVRP.
The remainder of this paper is organized as follows. Section 2 presents the GVRP. The clustering algorithms used in this study are discussed in section 3. A general overview of HNNs and HNNs method for GVRP is mentioned in Section 4. We present real-life examples and the proposed method results in Section 5. Finally, conclusions and discussions of future research directions are presented in section 6.

Green Vehicle Routing Problem (GVRP)
The GVRP aims to find feasible routes for delivering products from one depot to customers by minimizing the CO2 emissions with vehicles considering a certain capacity. The GVRP has many vehicles to meet the demands of customers. Vehicles that carries customer demands should not exceed the capacity constraint. The route of each vehicle starts and ends at the depot. Each customer must be visited only once by a single-vehicle. The problem aims to minimize CO2 emissions.
The GVRP is defined in a complete graph G = (N, A) where N = {0, 1, …, n} is the set of customers and "0" represents the depot. A corresponds to the edges connecting vertices of N. The notation is provided in Table 1. : Demand of customers in the subset S : CO2 emissions between customer i and j Decision Variable =Binary decision variable, taking a value of 1 if vehicle k goes from customer i to customer j and 0 otherwise.
We modified Toth & Vigo's (2002) mathematical model using the CO2 emission model presented by Nie and Li (2013). The mathematical formulation of the GVRP is as follows: Objective Function: Constraints: Constraint (1) is an objective function that aims to minimize the overall CO2 emission of the route. Constraint (2) and (3) guarantee that each customer is visited exactly once. Constraint (4) provides the flow constraint that the same vehicle enters and leaves the same depot in a route. Constraint (5) states that the maximum loading capacity of the vehicles is satisfied. Constraint (6) guarantees that each vehicle is used exactly once. Constraint (7) is used to eliminate sub-tours.

Solution Methodology
We propose a clustering algorithm and the HNN-based approach to solve large GVRP in a reasonable time since GVRP generalizes the TSP, which is NP-hard, and VRPs are generally much more difficult to resolve than TSPs of the same size. In our approach, we first divide the large-size GVRP into clusters using the K-Means and K-Medoids algorithms, and then vehicle routes are determined for each cluster using the HNN to minimize CO2 emissions.

Clustering analysis
Cluster analysis divides data into groups according to the determined similarity criteria. The clustering algorithm aims to minimize the similarity among groups and maximize the similarity of the within groups. Clustering techniques are divided into hierarchical and non-hierarchical clustering techniques. In this study, we use K-Means and K-Medoids clustering algorithms which are a member of non-hierarchical clustering algorithms.

K-Means algorithm
K-Means, one of the oldest clustering algorithms, was developed by J. B. MacQueen in 1967(MacQueen, 1967. The assignment mechanism of K-Means, one of the most widely used unsupervised learning methods, allows each data to belong to only one cluster. This method's main idea is that the central point represents the cluster (Han & Kamber, 2001). The Error Sum of Squares criterion (SSE) is most used to evaluate the K-Means clustering algorithm. The clustering result with the lowest SSE value gives the best results. The sum of the squares of the distances of the objects from the cluster's center points is calculated with Eq. (8) (Tan et al., 2006).
This criterion allows clusters to be determined densely and separately from each other. The algorithm tries to identify k groups that will reduce the error sum of squares function. The K-Means algorithm divides the data set consisting of n data into k clusters with the k parameter given to the algorithm by the user. Cluster similarity is measured by the mean value of the objects in the cluster, which is the gravity center of the cluster (Xu & Wunsch, 2005). The pseudocode for the K-Means algorithm is given in Algorithm 1.
Algorithm 1. The K-Means Algorithm 1: Begin 2: Determine the cluster centroids 3: Do 4: Calculate the distance between specified point and centroids 5: Assign all objects to the group that has the closest centroid 6: Recalculate the positions of the k centroids 7: Until convergence criteria are met 8: End begin

K-Medoids algorithm
Kaufman & Rousseeuw (1990) developed the K-Medoids clustering algorithm to eliminate the sensitivity of the K-Means clustering algorithm to noise and outline. There are many different variants of the K-Medoids algorithm. The first proposed K-Medoids algorithm is PAM (Partitioning Around Medoids). Although representative objects are often called centrotypes in the clustering literature, they are called medoids in the PAM algorithm. In the PAM algorithm, k randomly selected representative objects (Medoids) are chosen as the cluster's center. After k medoids are determined, k clusters are created by assigning each object to the closest medoids. Each medoid is replaced with a non-representative object in the following steps and shifted until improved clustering quality. Thus, the possibility of some data moving to the cluster center is eliminated (Kaufman & Rousseeuw, 1990). The pseudocode for the K-Medoids algorithm is given in Algorithm 2.
Algorithm 2. The K-Medoids Algorithm 1: Begin 2: Determine k cluster number 3: Choose k objects as initial medoids 4: Do 5: Assign the rest of the clusters to the closest medoids x 6: Calculate objective function 7: Randomly select y object that is not medoids 8: If changing x and y minimizes the objective function 9: Change x and y 10: End if 11: Until convergence criteria are met 12: End begin

Hopfield neural networks
Hopfield neural network structures are single-layer and fully connected neural network structures used mainly for autoassociation, developed by the physicist John Hopfield in 1982 (Hopfield, 1982).
All nodes in the HNN are fully interconnected. Each node is an input to every other node in the network. Connections from each node to itself can be thought of as a zero-weighted connection. HNN is a single-layer ANN with feedback and generally accepts bipolar (-1 or +1) or binary (0 or 1) inputs. Fig. 1 shows the HNN architecture (Rao, 1995).

Fig. 1. HNN architecture
HNN is divided into discrete and continuous according to their output values. Although the working styles are the same, the activation function is different. The sign function is used in discrete HNN. However, the hyperbolic tangent function is used in continuous HNN. While discrete HNN is used for self-association, continuous HNN is among the first heuristic approaches developed for solving combinatorial optimization problems for ANN.
In HNN, each variable in the problem is represented by a node. The connections between variables are also characterized by their weight values. Positive weight values show mutually supportive relationships between variables, while negative weight values represent inhibitory relationships. Optimization of the neural network is the process of minimizing the objective function. In HNN, the weight remains unchanged throughout the process, and the current situation of the neuron is related to the last moment. For this reason, the energy function satisfies the stability of HNN. The energy function is as follows (Equation (9)): where N is the number of neurons, is the weight of interconnect between the output of the neuron j and the input of the neuron i, is the output of the neuron i, is the external input signal of the neuron j (Hopfield & Tank, 1985).

Hopfield neural network-based method for GVRP
The proposed method's steps for the GVRP are given as follows: Step 1: The distance matrix among the customers' pairs is initialized.
Step 2: An initial transposition matrix of size NxN is defined to create the connection weights of the network. N represents the number of customers.
The transposition matrix of size NxN consists of 0 and 1 elements. Each row and column represent a particular city and position, respectively. The element "1" in the i th position of row A indicates that the visit order of city A is i. The transposition matrix of the four cities routing problem representing HNN and the route defined as the matrix is shown in Fig.2. The first row of matrix V in Fig. 2 represents city A. Vector (0, 0, 1, 0) means that city A is the third city to visit. The second vector (0, 1, 0, 0) indicates that city B will be visited in the second position, the third vector (0, 0, 0, 1) indicates that city C will be visited in the fourth position and the last vector (1, 0, 0, 0) indicates that city D will be visited first.
Step 3: The energy function is created as in Equation (10) to store the connection of weights (Hopfield & Tank, 1985 In Equation (10), positive constants are A, B, C, and D. It is crucial to choose constants correctly for the stability of HNN. N is the number of cities to be visited on route. , represents the CO2 emissions from city x to city y. , is the output of the i th neuron in the neuron array corresponding to the x th city. The final state of NxN variables consists of the terms 0 or 1 or very close to 0 and 1. If we analyze the energy function, the first and second terms give row and column constraints, respectively. The first term indicates that no city is visited more than once at the same time. The second term represents that each city is precisely once. The third term takes into care that all cities are visited. The first three constraints of the energy function provide a stable combination matrix. The final sum provides the constraint for the shortest CO2 emissions. This term's value is minimum when the total CO2 emissions are the shortest. The value of D is essential for deciding between the time is taken for convergence and the solution's optimality. A low D value means closer to the optimum solution, but it takes longer for the HNN to converge. If the D value is high, it converges quickly, but the solution deviates from the optimum.
Step 4: The connection weights of the network are calculated by equation (11).
where is the Kronecker's delta, which is defined as: Step 5: The network parameter values are entered, and the energy function is iterated for the number of iterations determined using the hyperbolic tangent function, one of the activation functions, to obtain the stable state of the neural network. The hyperbolic tangent function is given in Equation (12).
Here is the input of the neuron and is an initial value. calculates as Equation (13) (Wilson & Pawley, 1988).
= . + (13) Step 6: As a result of the iterations, in the transposition matrix, the position corresponding to the maximum value of each column is set to 1, the others to 0.
Step 7: It is checked whether the created transposition matrix satisfies the constraints of the GVRP problem. If so, the following steps are continued. If not, go to Step 5.
Step 8: The route created according to the transposition matrix is divided by considering the customer demands and the vehicle capacity constraint, and the vehicles are assigned.

A case study
In this study, the problem is to find the vehicle routes that minimize the total CO2 emissions to meet the store demands of the vehicles leaving the main depot in a supermarket chain operating. The supermarket chain has one main depot and 78 stores. The location of the main depot and stores on the map is shown in Fig. 3. The CO2 emissions of inter-store and stores between the main depot are calculated according to the Comprehensive Emission Model proposed by Nie & Li (2013). This problem includes a homogeneous fleet with a 40 pallets capacity. The vehicle fleet has two types of trucks: fresh and dry products for a distribution process. 78 stores are divided into two as daily and not daily stores. The daily stores demand products every weekday. Whereas not daily stores only demand products three days a week (Monday, Wednesday, and Friday). The number of daily and not daily stores are 66 and 12, respectively. The numbers of the daily and not daily stores are given in Table 2.  2,3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,25,26,28,29,31,32,34,35,37,38,40,41,43,44,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,72,73,75,76,78 Not daily 6,24,27,30,33,36,39,42,45,71,74,77 Days are classified as "busy" and "not busy" days. Distributing all stores is called a "busy day", whereas distributing only not daily stores is called a "not busy day". The problem is divided into four sub-problems based on the distribution days of the products.  Table 3, the demands of the sub-problems in the first week. Before using the proposed HNN to solve the problem, 78 and 66 stores are clustered with K-Means and K-Medoids clustering algorithms. First, the number of clusters (k) is determined for the K-Means and K-Medoids algorithms, as explained in Section 3. The Elbow method is used to determine the k number in this study. The elbow method is a method that deals with the percentage of variance, which is a function of the number of clusters. This method is the most preferred method that shows us the big picture. Within cluster sum of square (WCSS) is calculated by taking the sum of the squared distance of each point from the cluster centroid. The Elbow method states that the point where the amount of change in the WCSS decreases, the elbow point, is the optimum point (Salvador & Chan, 2003). The optimum k value is determined as seven with the Elbow method. The Elbow graph of stores distributed on busy days and not busy days are shown in Fig. 4.

Table 9
Friedman Statistics for a sub-problem 1 (C1) Total N 20 Test Statistic 0.000 Degrees of freedom 1 p-value 1.000 • Comparison in terms of using K-medoids algorithm (C2) The null and alternative hypotheses are determined as follows: H0: There is no significant difference between the performance of the two methods (HNN-C2 and CFRS-C2) with a 95% confidence interval. Ha: There is a significant difference between the performance of the two methods (HNN-C2 and CFRS-C2) with a 95% confidence interval. Table 10 shows the Friedman test results. The Friedman test statistics are highly significant (p≥0.05), so the null hypothesis can be accepted, and the H0 hypothesis is accepted. In other words, there is no significant difference between the two methods (HNN-C2 and CFRS-C2).

Table 10
Friedman Statistics for a sub-problem 1 (C2) Total N 20 Test Statistic 0.800 Degrees of freedom 1 p-value 0.371 The other sub-problem' Friedman test results, are given in Table 11.

Table 11
Friedman Test results of sub-problems in solving VRP has not yet been addressed in the literature. Here, we compared the solution performance of HNN with CFRS proposed by Comert et al. (2018) and concluded that HNN obtained competitive results.

Conclusion
In this study, a newly proposed method based on combining K-Means and K-Medoids algorithms with HNN is developed to solve large-scale problems. The method is adopted for GVRP. Firstly, stores are grouped into clusters using K-means and Kmedoids algorithms. Thus, we enable the solution of large-scale problems. Then, we route each cluster with HNN to minimize CO2 emissions. The proposed HNN-based method is applied to a real-life example to illustrate its effectiveness. This real-life problem is related to a distribution GVRP meeting the weekly demands of a supermarket chain with 78 stores and one depot. The results obtained using HNN and CFRS proposed by Comert et al. (2018) are compared using the Friedman test. The test result shows the applicability and effectiveness of the developed HNN-based method in VRP. Future work may consider electric vehicle routing problems from alternative fuel vehicles with the proposed HNN and clustering algorithms-based method.