An Improved Chimp-Inspired Optimization Algorithm for Large-Scale Spherical Vehicle Routing Problem with Time Windows

The vehicle routing problem with time windows (VRPTW) is a classical optimization problem. There have been many related studies in recent years. At present, many studies have generally analyzed this problem on the two-dimensional plane, and few studies have explored it on spherical surfaces. In order to carry out research related to the distribution of goods by unmanned vehicles and unmanned aerial vehicles, this study carries out research based on the situation of a three-dimensional sphere and proposes a three-dimensional spherical VRPTW model. All of the customer nodes in this problem were mapped to the three-dimensional sphere. The chimp optimization algorithm is an excellent intelligent optimization algorithm proposed recently, which has been successfully applied to solve various practical problems and has achieved good results. The chimp optimization algorithm (ChOA) is characterized by its excellent ability to balance exploration and exploitation in the optimization process so that the algorithm can search the solution space adaptively, which is closely related to its outstanding adaptive factors. However, the performance of the chimp optimization algorithm in solving discrete optimization problems still needs to be improved. Firstly, the convergence speed of the algorithm is fast at first, but it becomes slower and slower as the number of iterations increases. Therefore, this paper introduces the multiple-population strategy, genetic operators, and local search methods into the algorithm to improve its overall exploration ability and convergence speed so that the algorithm can quickly find solutions with higher accuracy. Secondly, the algorithm is not suitable for discrete problems. In conclusion, this paper proposes an improved chimp optimization algorithm (MG-ChOA) and applies it to solve the spherical VRPTW model. Finally, this paper analyzes the performance of this algorithm in a multi-dimensional way by comparing it with many excellent algorithms available at present. The experimental result shows that the proposed algorithm is effective and superior in solving the discrete problem of spherical VRPTW, and its performance is superior to that of other algorithms.


Introduction
The vehicle routing problem (VRP) is a famous path-planning problem that was first proposed by Dantzig and Ramser [1]. It has been widely studied in the field of optimization problems and is a very practical model (Toth et al. [2]). In recent years, many variants of the VRP have been developed, such as the VRP with time windows (VRPTW) (Yu et al. [3]), green VRP (GVRP) (Xu et al. [4]), capacitated VRP (CVRP) (Zhang et al. [5]), multi-depot VRP (MDVRP) (Duan et al. [6]), and heterogeneous VRP (HVRP) (Ghannadpour et al. [7]). To sum up, there are mainly three kinds of classical algorithms to solve VRPs, which are exact algorithms, traditional heuristic algorithms, approximation algorithms, and metaheuristic algorithms. Exact algorithms (such as branch-and-price (Li et al. [8]), sophisticated branch-cut-and-price methods (Pessoa et al. [9]), mixed-integer nonlinear programming algorithms (Xiao et al. [10]), approximate dynamic programming algorithms (Çimen et al. [11], etc.) use mathematical methods to search for the optimal solution. Although they can often find a good solution, there are also problems, such as the single form of the solution, inability to avoid exponential explosion, and consumption of a lot of computing time. The basic idea of traditional heuristic algorithms (such as saving algorithms (Li et al. [8]), improved Dijkstra algorithms (Behnk et al. [12]), etc.) is to start from the current solution, search for a better solution in its neighborhood to replace the current solution, and then continue to search until there are no better solutions (Da Costa et al. [13]). The traditional heuristic algorithm easily falls into the local optimum and cannot easily achieve the global optimum. In addition, approximation algorithms (Das et al. [14], Khachay et al. [15]) with theoretical performance guarantees and approximation schemes have been widely used to solve the problems of covering non-Euclidean settings. Metaheuristic algorithms have excellent performance and constantly perturb near the current solution to search for better solutions. Therefore, this study uses the meta-heuristic algorithm to solve the problem.
However, many studies on the VRPTW are generally carried out on two-dimensional planes. In many research fields, the three-dimensional spherical structure is also of great significance. For example, celestial bodies, particle structures, daily foods, proteins and other nutrients, balls, buildings, and path-planning problems are all problems related to spheres. Therefore, it is of great practical significance to expand the research on the VRPTW from the two-dimensional plane to three-dimensional spheres. Zhang et al. [59] applied the BBMA to the spherical MST problem. To sum up, in order to carry out research related to the distribution of goods by unmanned vehicles and unmanned aerial vehicles, this paper plans the path of a robot through these coordinates in three-dimensional space in order to carry out research on the VRPTW in the spatial dimension.
ChOA has been applied to various practical problems, such as image segmentation for medical diagnosis (Si et al. [60]), clustering analysis (Sharma et al. [61]; Yang et al. [62]), Said-Ball curve degree optimization (Hu et al. [63]), and convolution neural networks (Chen et al. [64]). In addition, Du and Zhou (Du et al. [65]; Du et al. [66]) improved this algorithm and applied it to 3D path-planning problems and color image-enhancement problems. The algorithm divides the population of each generation into four groups, namely attackers, barriers, chasers, and drivers, and they cooperate against prey. Therefore, different groups search different spaces, which enhances the searching ability of ChOA. The adaptive factor of ChOA has a faster convergence speed and can adaptively balance exploration and exploitation, but it also easily falls into the local optimum. In order to obtain a group of better solutions with limited resources and time, this paper proposes an improved ChOA (MG-ChOA) for solving the spherical VRPTW model. The main contributions of this paper are as follows. Firstly, the proposed algorithm combines the ChOA algorithm with the quantum coding, local search, multiple population, and genetic operators to ensure that the algorithm can not only achieve adaptive and rapid convergence, but also find solutions with higher accuracy. Secondly, this paper proposes a three-dimensional spherical VRPTW and applies the proposed algorithm to solve this problem. Finally, by comparing with the running result of popular swarm intelligence algorithms for eight different instances, the effectiveness and superiority of the proposed algorithm in dealing with large-scale combinatorial optimization problems are strongly verified.
The remaining organizational structure of this paper is as follows. Section 2, the Related Work, briefly depicts works related to the model proposed. Section 3 analyzes twodimensional VRPTW and spherical VRPTW models to propose the mathematical model of a spherical VRPTW. The proposed algorithm (MG-ChOA) for a spherical VRPTW, an improved MG-ChOA algorithm based on ChOA, is presented in Section 4. The discussion of the experimental results analyzes the performance of the algorithm in Section 5. The conclusion and future work proposals are presented in Section 6.

Geometric Definition of Sphere
A sphere refers to a set of points in 3D space with equal distance from the center point of the sphere, and radius is the distance from the center point of the sphere to a point on the sphere, as shown in Figure 1a. Therefore, a sphere with radius r can be defined by the following formula.
where x, y, and z are coordinate axes of three-dimensional space, which are used to describe each point.

Definition of Points on the Sphere
The coordinates of points on the sphere can be described in detail with the following equation (Hearn et al. [67]).
The coordinate of each point can be represented by x, y, and z, and they can be expressed by normalized parameters (such as u and v) in [0, 1]. Equations (3)-(5) specify the coordinates of each point on the sphere (Eldem et al. [68]).
where u and v respectively represent the longitude and latitude lines used to calibrate the position, as shown in Figure 1b. Different combinations of u and v describe points of the sphere (Ugur et al. [69]), as shown in Figure 1b. In order to save computing resources and compare the performance of algorithms more conveniently, this study uses a sphere with a radius of one to carry out experiments and discussions.

Geodesics between Two Points on the Sphere
The big circle is a figure formed by the intersection of a plane passing through the sphere's center. The shortest path between two points on the sphere is a certain arc length of the big circle, and the geodesic line is this arc (Lomnitz [70]). According to the description above, the geodesic line between p i (point i) and p j (point j) on the sphere is shown in Figure 1c, which can be described by vector → v i and vector → v j , respectively, and their product is defined as follows.
where θ indicates the angle of two vectors, and the formula of the shortest path can be expressed as follows.
According to Formulae (6)-(8), we get ∧ d p i p j = r arccos( Calculate the distance between two points on a sphere with n points to obtain a symmetric distance matrix Dis with the size of n × n, and the distance matrix can be described as follows.
where ∧ d ij denotes the distance of the geodesic line formed by two points on the sphere and ∧ d ij equals infinity, meaning that the point cannot reach itself.

VRPTW on a 2D Plane
VRPTW describes a path-planning problem of allocating goods from a distribution center to different customers within a specified time, including K vehicles and N customer nodes. The logistics network reasonably plans a series of routes to serve customers according to the transportation demand, and vehicles must leave and finally return to the depot. As the loading capacity of each vehicle is limited, the transportation company must serve customers within a specified time to meet the customers' needs. In addition, if the vehicle arrives at the customer node in advance, it needs to wait for a period of time before starting the service. Figure 2 shows the process of VRPTW. The parameters used in the model are shown in Table 1.
Subject to: ∑ i∈C x ijk = y jk , ∀j ∈ C, ∀k ∈ V Equations (11) and (12) represent objective functions consisting of the travel distance and the number of vehicles. Equation (13) indicates that each customer's demands on a route should not be greater than the maximum capacity of the vehicle. Equation (14) denotes that one customer should only be served once. Equation (15) indicates that the vehicle only serves one node before serving the next node. Equation (16) denotes that vehicles only visit one node after visiting the previous node. Equation (17) indicates that the start and end nodes of each vehicle should return to the depot. Equations (18) and (19) are time window constraints.

Three-Dimensional Spherical VRPTW Model
The 3D spherical VRPTW model maps the customer nodes from the 2D plane model to the 3D space. Therefore, each customer node of i can be defined as c i = (x i , y i , z i ). Therefore, the mathematical model of the three-dimensional spherical VRPTW model can be defined as follows: In the formula above,  The ChOA was first proposed by Khishe and Mosavi [17] in 2020, and it simulates the predatory behavior of chimps. The algorithm divides the population of each generation into four groups, namely attackers, barriers, chasers, and drivers, and they cooperate against prey. The ChOA has excellent adaptive factors that can help it balance exploration and exploitation so as to find better solutions. The mathematical model of ChoA is as follows: Equations (21) and (22)    The model above describes the main flow of the algorithm. In each iteration, the algorithm firstly selects the four best individuals, and then the remaining individuals update their positions based on them. The specific mathematical model of the algorithm is as follows: The random vector → c strengthens (c > 1) or weakens (c < 1) the moving range of prey. When the random vector → a is greater than 1 or less than −1, the algorithm will be in the exploration stage; otherwise, it will be in the exploitation stage. Therefore, the algorithm can adaptively adjust exploration and exploitation to find a better solution. The pseudo code of ChOA is shown in Algorithm 1. for each agent 10.
Update its position based on Equations (26)

The Proposed MG-ChOA for the Spherical VRPTW
ChOA can adaptively adjust exploration and exploitation and its convergence speed is fast, but it still has the shortcomings of limited exploration capacity, easily falling into the local optimum, and not being suitable for discrete problems. Therefore, this paper improves the ChOA algorithm to solve the three-dimensional spherical VRPTW. The proposed algorithm (MG-ChOA) introduces the quantum coding method, multiple-population strategy, genetic operators, and local search strategy to improve the search ability of the algorithm. Using the quantum coding method to initialize the population can increase the population's diversity at the initial stage. Similarly, the multiple-population strategy can increase the population diversity of the algorithm in the iterative process to find better solutions. Genetic operators enhance the exploration ability of the algorithm, and the local search strategy helps the algorithm to search for better solutions in the neighborhood of each solution.

Encoding and Decoding of the Spherical VRPTW
As shown in Figure 3, each code is divided into many small parts by the number of 0 that represents the distribution center, and the remaining numbers represent customer nodes. Each independent small part above represents one path served by one vehicle. Therefore, the path encoding of Figure 3 can be sequentially decoded into sequences of [1,2,3], [4,9,8], and [7,6,5,10].

Initializing Population Using the Quantum Coding
In this study, we used quantum coding to initialize the population. This method enhances the population's diversity at the beginning of the iteration, which is conducive to enabling the algorithm to quickly find a better solution at the initial stage. The smallest unit of quantum computing is the quantum bit (qubit), whose state is the superposition state of 0 or 1. The quantum bit is defined as follows.
where α and β are complex numbers and α 2 and β 2 respectively represent the probability of the state being 0 or 1, and they satisfy the equation of α 2 + β 2 = 1. The specific representation of a qubit is as follows: Inspired by the quantum encoding, individuals of the population can be described as follows: where QC i denotes the ith individual of the population, θ denotes the angle falling in the interval of [0, 2π], and n and d represent the total number and dimensions of each individual, respectively. Therefore, each individual has two candidate schemes, which can be defined as follows: The initialization population of the proposed algorithm consists of n individuals and d dimensions, so the initialization population is to construct a matrix of n × d. According to the rules of quantum coding, we firstly needed to initialize the angle matrix and then obtain the initialization population according to the angle matrix. The angle matrix can be calculated by the following formula: where lb ij and ub ij respectively represent the minimum and maximum values of each of the individuals' dimensions at the problem boundary, and their values are 0 and 2 π, respectively, while rand indicates random numbers falling in the interval of [0, 1]. The initialized angle matrix is as follows.
where QC represents a population of N individuals. Each individual has two positions in the space, and they represent the candidate solution of the problem. Therefore, each individual has two different candidate solutions. To sum up, the quantum population can be defined by the following formula.
According to Equation (35), the distribution of individuals in the initial population is closely related to two different trigonometric functions. The value distribution of sine and cosine functions greatly affects the distribution of the population, which can not only enhance the diversity of the population distribution in the solution space, but also balance the initial exploration and exploitation; therefore, this is helpful for the algorithm to find excellent solutions quickly. In addition, individuals obtained by the formula above are all decimals in the interval of [0, 1]. Therefore, the proposed algorithm converts these decimals into percentages and then multiplies them by the total number of customers, rounds them, and obtains the route code. Finally, we removed and reinserted the duplicated part of the route code to obtain the final initialization population.

The Multiple-Population Strategy for MG-ChOA
The main method of the multiple-population search strategy is to use multiple populations with different parameters to search the target solution together, and it mainly includes the immigration and manual selection operator. When the algorithm iteratively searches for the target solution, the migration operator is used to contact and exchange the optimal solution among various populations, and the optimal solution of each generation is saved to the essence population by the manual selection operator as the basis of algorithm convergence. The whole process above is shown in Figure 4. As shown in the figure, when the population 1-N searches the solution space, the migration operator will replace the worst individual of population i + 1 with the best individual of population i, and the best solution of population N will be used to replace the worst solution of the first population. Therefore, information exchange between populations is achieved.

Genetic Operators
The proposed algorithm introduces genetic operators to help the algorithm update the population to search for excellent individuals in the solution space. Specifically, the proposed algorithm firstly selects some excellent individuals as a new population, then assigns a probability pd to each individual in the new population; if it is greater than or equal to 0.5, the current individual will cross with the next individual. Figure 5 shows an example of the crossover operation, in which the numbers 3-8 represent the selected individuals, and then the algorithm assigns a random value pd to each of them. From the method above, pd of 3 and 5 are numbers greater than 0.5, so they were selected to cross with their next individual. Moreover, mutation is also an integral part of genetic operators. Similarly, the proposed algorithm assigns a probability pmd to each individual in the population, and if it is greater than the mutation probability pm given in this paper, the individual will perform the mutation operation. Therefore, a new route can be obtained by exchanging randomly selected customers on two individuals. The introduction of genetic operators enhances the exploration ability of the proposed algorithm.

Local Search Strategy
The local search strategy is helpful for the proposed method to search in the neighborhood of a specific solution to find a better solution. This method allows customers to move on the path in a certain way, and finally obtains more potential individuals without violating various constraints. The algorithm proposed in this paper uses a mechanism to generate a number of customers to be removed in the search process, and then reinserts the removed customers in a more reasonable location without violating the constraints of the time window and capacity. Therefore, better solutions in the neighborhood will be obtained. Figure 6 shows that a better route was obtained by randomly removing some customers and reinserting them in a more appropriate location.

Experimental Results and Discussion
This paper uses eight different cases to test MG-ChOA's ability to solve spherical VRPTW. All experiments were conducted on the premise that r equals one, and the numbers of customer nodes in these instances were 80, 100, 200, 400, 600, 800, 1000, and 1200. For cases with fewer than 100 customers, the time window and load data were from the Solomon dataset (Solomon [71]), and these data of the remaining cases were generated randomly. Owing to the strong randomness of metaheuristic algorithms, all algorithms were independently run 30 times in each case to obtain the result. The structure of this chapter is as follows. Firstly, Experimental Setup provides the parameter settings of all algorithms and experimental configurations. Secondly, this paper presents the results of the proposed algorithm running for 2D instances. Moreover, this paper compares the search ability and results of the proposed algorithm with other algorithms in low-and high-dimensional cases. Finally, this paper discusses the impact of different improvement methods on the overall performance of the algorithm.

Experimental Setup
The code of all experiments carried out in this paper was compiled in MATLAB. The computer system was configured with an Intel Core Intel (R) Core (TM) i7-9700 CPU, 16 GB RAM, and Windows 10 operating system. In this experiment, each algorithm iterated 300 generations and the population size of all algorithms was set to 100. This paper also compares the performance of MG-ChOA with the GA, ant colony algorithm (ACO), PSO, slime mold algorithm (SMA), firefly algorithm (FA), chimp optimization algorithm (ChOA), and gray wolf optimizer (GWO). The GA was integrated with the local search strategy. Moreover, many excellent improved algorithms (RPSO (Borowska [72]), JADE (Su et al. [73]), L-SHADE (Chen et al. [74]), learning CSO (Borowska [75]), and CMA-ES (Tong et al. [76])) should also be used to comprehensively analyze the performance of the proposed algorithm, and this study selected two of them (RPSO, JADE). In order to verify the effectiveness and superiority of the MG-ChOA algorithm, this paper also comprehensively compares the convergence curve, ANOVA test, fitness value obtained by 30 runs, Wilcoxon rank-sum test (Gibbons et al. [77]; Derrac et al. [78]), effects of different improvement methods, and running results. In addition, the control parameter for each algorithm was set as follows ( Table 2) (Zhang et al. [59]). The pheromone was set to 4, heuristic information was 5, waiting time was 2, time window width was 3, parameter controlling ant movement was 0.5, evaporation rate of pheromone was 0.85, and the constant affecting pheromone updating was 5 PSO Kennedy et al. [32] The inertia weight was 0.2, global learning coefficient was 1, and the self-learning coefficient was 0.7

RPSO Borowska [72]
The inertia weight was 0.6, acceleration constants were c1 = c2 = 1.7, and the number of particles with the worst fitness p was set as 3.

JADE
Borowska [73] Parameters of the algorithm changed adaptively SMA Li et al. [79] The parameter controlling foraging was 0.03 FA Yang [35] The basic value of the attraction coefficient was 0.8, the mutation coefficient was 0.8, and the light absorption coefficient was 0.8 ChOA Khishe et al. [25] Parameters of the algorithm changed adaptively GWO Mirjalili et al. [37] Parameters of the algorithm changed adaptively

Performance Comparison of Algorithms for Two-Dimensional Datasets
In order to analyze the performance of the proposed algorithm in multiple dimensions, this chapter briefly analyzes the results of the proposed algorithm in the two-dimensional plane. This research selected four different Solomon datasets to test the algorithm, which were C101, R102, R201, and RC105. The differences between the algorithm and the most famous results are provided in Table 3. Figure 8 shows the convergence speed of this algorithm. Figure 9 shows the optimal path obtained by the algorithm for each dataset. Table 3 shows that the performance of the algorithm for four different datasets was significantly different from the most famous results, with minimum and maximum values of 0 and 7%. Figure 8 shows the superior convergence performance of the algorithm, which converged to the optimal value quickly for all datasets. In conclusion, the results of four 2D datasets show the effectiveness of the proposed algorithm.

Performance Comparison of Algorithms for Low-Dimensional Instances
This part tests the performance of the algorithm through four instances of different scales. All algorithms were independently run 30 times to obtain the results, and then the optimal value, the worst value, the average value, and the standard deviation of the result obtained were recorded. As shown in Table 4, the font marked in bold indicates the best value of all algorithms. Figure 10 shows the convergence curve after 300 iterations of the algorithms. Figure 11 shows the ANOVA test of the result after 30 runs. Figure 12 shows the result of 30 runs by algorithms, Figure 13 shows the optimal path found by the proposed algorithm, and Table 5 shows the Wilcoxon rank-sum test, which shows the significance of the difference between the proposed algorithm and other algorithms. It is worth noting that the rank metric in Table 4 was obtained by the Friedman statistical test (Zimmerman et al. [80]), and the fitness values shown in Figures 11, 12, 15, 16, 19 and 20 could be calculated according to the following formula, which represents the distance from the coordinate composed of TD and NV to the coordinate origin. In addition, as the convergence curve of TD is more representative, this paper will mainly analyze the convergence speed of algorithms according to the TD.       Fitness Value = (TD) 2 + (NV) 2 It can be seen from the instance with 80 customers that the comprehensive performance of MG-ChOA was better than that of the other algorithms. It obtained the optimal value and the minimum average value, but did not obtain the optimal standard deviation. The convergence curve in Figure 10 displays the excellent convergence ability of the proposed algorithm, which converged to the optimal solution faster than the other algorithms. In addition, GA ranked second overall, and the performance gap between algorithms was very small. The ANOVA test in Figure 11 shows that the result obtained by the proposed algorithm was relatively uniform, which means that it had strong stability and good optimization ability. The running result of Figure 12 shows that the search ability of the proposed method was relatively stable, and results obtained were better than those of other algorithms. Figure 13 shows an optimal path obtained by the proposed algorithm for this instance.
From the instance with 100 customers, we can see that MG-ChOA's comprehensive performance was superior to that of other algorithms, and it obtained superior average and optimal values to all algorithms. ChOA obtained the optimal standard deviation, which indicated that its performance was stable for this instance. The convergence curve in Figure 10 indicates that the proposed algorithm had excellent convergence ability for the instance of 100 customers. It not only found a better solution, but also had the fastest convergence speed. GA ranked second overall, and the solution searched by GA was only second to the first. The ANOVA analysis in Figure 11 shows that the comprehensive performance of the proposed algorithm was relatively good for this instance, and it had stable search capability. The running result of Figure 12 shows that results obtained by the algorithm were generally better than those of the other algorithms. Figure 13 shows the excellent solution for the proposed algorithm for this instance.
The comparison result of the instance with 200 customers indicated that the proposed algorithm had the strongest capability, and not only was the minimum average value obtained, but also the search result was the best among all algorithms. GA ranked first overall-its performance was more stable than that of the proposed algorithm-and PSO achieved the best standard deviation. The convergence curve in Figure 10 shows that the convergence ability of the proposed algorithm was faster and better than that of the other algorithms for this instance. The ANOVA analysis in Figure 11 shows that, although the search ability of the proposed method was not stable, it could find better solutions than the other algorithms. The running result in Figure 12 shows that results obtained by the algorithm were generally better than those of the other algorithms, and the gap was increasingly obvious. Figure 13 shows a good path obtained by the proposed algorithm for this instance.
The comparison result of the instance with 400 customers indicates the superior ability of the proposed algorithm in this instance, and GWO achieved the best standard deviation; in addition, RPSO ranked second overall. Compared with small-scale instances, except for the genetic algorithm, the other algorithms were relatively stable for this instance, but their search ability was greatly reduced. The proposed algorithm and GA both hadbetter search ability and they found better solutions, but the proposed algorithm had better convergence ability and could find good solutions faster than GA. The ANOVA analysis in Figure 11 indicates that the proposed algorithm's performance was very strong, although it was not as stable as other algorithms. It had the best optimization ability, and the gap between GA and other algorithms is getting smaller and smaller. Figure 13 shows a feasible route found by the proposed algorithm for this instance.
To sum up, the comparison result above shows that there was little difference in the capability of the algorithms for small-scale instances. When the number of customers increased, the performance of the algorithm started to show a significant difference. Compared with small-scale instances, differences between algorithms for large-scale instances increased significantly. The search ability of the proposed algorithm was still excellent for large-scale instances, and it had fast convergence and search capability. The performance of GA ranked second in the instance above, which was related to the strong search ability of its genetic operators. Both ChOA and GWO are stable because they have adaptive factors. With the increase in the instance's size, the search ability of other algorithms became weaker and weaker, and the gap between the genetic algorithm and other algorithms became smaller and smaller. In general, the proposed algorithm performed best in the four instances, but the stability of the algorithm needs to be improved, which is related to the richness of improved methods. Finally, Table 5 shows the significant difference between the proposed algorithm and other algorithms in each instance, with a significance level of 0.05. This means that the difference between two groups of data is significant if the p value is less than 0.05. Table 5 shows that the experimental results between the proposed algorithm and other algorithms were significantly different in low-dimensional instances. Therefore, the proposed method is better than other algorithms.

Performance Comparison of Algorithms for High-Dimensional Instances
In this section, the search ability of the algorithms is tested through four instances of different scales, namely, instances with 600, 800, 1000, and 1200 customers. All algorithms were independently run 30 times. Table 6 shows the best value, the worst value, the mean value, and the standard deviation of the result obtained, and the font marked in bold represents the optimal value obtained by all algorithms. Figure 14 compares the convergence curve of the algorithm after iterating for 300 generations. Figure 15 shows the result of the ANOVA analysis of the algorithm after 30 runs, Figure 16 displays the result of these algorithms for 30 runs, Figure 17 shows the optimal path of the proposed algorithm, and Table 7 displays the Wilcoxon rank-sum test, which shows the significance of the difference between the proposed algorithm and other algorithms.     From the instance with 600 customers, we can see that MG-ChOA had a better search ability than the other algorithms. It obtained the best value and the optimal average value, and RPSO obtained the best standard deviation; in addition, RPSO ranked second. RPSO and JADE also performed well and had little difference from each other. ACO, ChOA, and GWO had relatively poor search ability. The convergence curve in Figure 14 shows that the proposed method achieved faster convergence than the other algorithms. The ANOVA analysis in Figure 15 shows that the search ability of the proposed method was relatively stable in all algorithms, and the solutions found were better than those of other methods. The running result in Figure 16 shows that, although the search ability of the proposed algorithm fluctuated greatly, the result obtained was generally better than that of other algorithms. Figure 17 shows an excellent solution for the proposed algorithm for this instance.
From the instance with 800 customers, we can see that MG-ChOA obtained the optimal and superior average values to all algorithms, and its comprehensive performance was the best; moreover, RPSO ranked second overall, and JADE obtained the best standard deviation. The convergence curve in Figure 14 indicates that, for this large instance, the excellent convergence speed of the proposed method enabled the algorithm to constantly search for better solutions; in addition, the convergence speed of the proposed method was much faster than that of other algorithms and the solution found was better. Secondly, the search ability of the chimp optimization algorithm was better than that of the GA, which indicates that its adaptive factor well-balanced the search process of the algorithm. The ANOVA test in Figure 15 shows that the search ability of the proposed method was relatively stable for this instance, and the JADE's performance was the most stable. The running result in Figure 16 shows that results obtained by the proposed algorithm were much better than those of other algorithms, but the performance of most algorithms was unstable. Figure 17 shows an excellent solution found by the proposed algorithm for this instance.
The comparison result of the instance with 1000 customers indicates that the proposed algorithm had the best search ability. It not only found the lowest average value, but also obtained the best result among all algorithms, and the GA ranked second overall. The convergence curve in Figure 14 shows that the proposed algorithm had the best convergence ability for this instance, while the convergence ability of other algorithms became slower and slower; this indicates that the proposed algorithm was still effective for large instances. The ANOVA analysis in Figure 15 reveals that all algorithms were stable, and the proposed algorithm found better solutions than the other algorithms. The running result in Figure 16 shows that the result obtained by the proposed algorithm was generally better than that obtained by the other techniques, and the performance gap between them was obvious. Finally, Figure 17 shows a feasible path obtained by the proposed algorithm for this instance.
The comparison result of the instance with 1200 customers indicates the superior ability of the proposed algorithm for this instance. The performance of other algorithms was relatively stable for this instance, and their search ability was greatly weakened compared with that for small-scale instances. Compared with other algorithms, the proposed algorithm still had a good convergence ability and obvious advantages, and it could find excellent solutions faster. The analysis of the ANOVA test in Figure 15 indicates that these algorithms were relatively stable with little difference from each other. The proposed algorithm still had the best performance, and the gap with other algorithms was the largest. Figure 17 shows a good solution obtained by the algorithm for this instance.
To sum up, the above comparison result shows that the differences in the search ability of these algorithms for large-scale instances became more and more obvious, but the performance of the proposed algorithm was still excellent, with fast convergence ability and good search ability, and the gap was about 30% higher than that for instances with 80 customers. Therefore, this shows that the adaptive factor of the algorithm well-balanced the search ability, and the combination of genetic operators and local search strengthened the exploration and convergence ability of the algorithm; in addition, the multiple-population strategy strengthened the communication between populations, allowing the algorithm to find better solutions faster. Secondly, the performance of the GA generally ranked second, which indicates that its genetic operator had strong search ability. The adaptive factor of the chimp optimization algorithm could well-balance the exploration and exploitation. Therefore, these two methods strengthened the performance and stability of the algorithm. With the increase in the instance's size, the performance of these algorithms became weaker and weaker, and the gap between other algorithms became smaller and smaller. In general, the proposed algorithm performed best for the above eight instances, but the stability of the proposed method needed to be improved, indicating that the improvement method for the proposed method was successful and efficient. Table 7 shows that the experimental results between the proposed algorithm and other algorithms were significantly different for high-dimensional instances. Therefore, the proposed method was obviously better than other algorithms.

Performance Limit Test of MG-ChOA
In order to measure the boundaries of the algorithm proposed in this paper, this section adds two additional instances to verify the limits of the algorithm, which are instances of 1600 customers and 1800 customers. Table 8 records the best value, the worst value, the average value, and the standard deviation of the algorithm for these instances. Next, this section uses the GA for comparison with the proposed algorithm in terms of their convergence speeds. Figure 18 shows the convergence speed of the algorithm. Figures 19 and 20 respectively show the ANOVA test and the result statistics of the algorithms when run 30 times. To sum up, the convergence speed of the algorithm became slower and slower. Compared with the 1200-customer instance, the 1600-customer instance had a significant decline, and the convergence ability of the proposed algorithm on the 1800-customer instance was close to the limit. Moreover, the analysis of the ANOVA test of the algorithm showed that its stability became weaker and weaker. Therefore, there is still much more room to improve the performance of the algorithm. Table 9 shows that the results of these two algorithms were significantly different. Figure 21 shows an excellent solution for the proposed algorithm for these instances

Performance Analysis of Different Strategies
This section analyzes the impacts of different improved methods on the overall performance of the proposed algorithm to further test the effectiveness of these improved methods. Figure 22 shows the search results of various improved methods on the 80-customer and 100-customer instances. In this paper, we propose initializing the population by using quantum coding. The main idea is as follows. Firstly, the population generated by the quantum coding technique is twice as large as the original population, which increases the diversity of the initial population. Secondly, this method mainly includes two trigonometric functions, which generate individuals with good probability distribution, and this method well-balances the initial exploration and exploitation of the algorithm. As can be seen from Figure 22, the method of initializing the population by quantum coding sped up the convergence of the algorithm and caused it to find a feasible solution faster, but the effect was not obvious. Moreover, genetic operators and the local search strategy were introduced into the algorithm. Genetic operators have strong search ability and can broaden the search scope of an algorithm. Their performance is very suitable for solving path problems, while local search algorithms can find better solutions in the neighborhoods of these excellent solutions found in each iteration. Therefore, the combination of the above two strategies can greatly improve the exploitation and exploration ability of algorithms and help them to find solutions with high accuracy. Figure 22 shows that this method made the largest proportion of contributions to the performance of the proposed algorithm, which shows that this method was suitable for this algorithm and enhanced the effectiveness of the proposed algorithm. This paper also introduced a multiple-population strategy and generated two populations. This method achieved communication between populations through migration operators, so the algorithm could find multiple excellent solutions in each generation, which enhanced the exploration ability. Figure 22 shows that the introduction of this method significantly improved the convergence speed of the proposed method, indicating that the proposed method could be feasibly integrated with multiple-population strategies. In addition, this study also conducted experiments to select the best population number, and the results are shown in Table 10, where p is the population's number. Table 10 shows that the optimal result was obtained when the population number was set to 2, so the population number of the method proposed in this paper was 2.

Conclusions and Future Work
The chimp optimization algorithm (ChOA) is a new swarm intelligence algorithm that has excellent search ability and is suitable for solving continuous problems. The characteristic of this algorithm is that the convergence speed is fast during the initial stage of iteration, and the solution's accuracy is high, but the convergence ability is weakened during the later stage of iteration, and it easily falls into the local optimum. Furthermore, the algorithm can adaptively adjust its exploration and exploitation when searching the solution space because the algorithm has well-designed adaptive factors to balance the exploitation and exploration in the process of optimization.
Although the performance of the chimp optimization algorithm itself was superior, it was not suitable for dealing with discrete optimization problems in real life, and the convergence speed of the algorithm became slower and slower, so it easily fell into the local optimum. Therefore, this paper improved the performance of the algorithm according to the shortcomings above, and the multiple-population strategy, genetic operators, and local search strategy were integrated into the algorithm to enhance the overall exploration ability and convergence speed of the proposed method. The multiple-population strategy initializes multiple populations, uses migration operators to exchange information among various populations, and finally selects excellent individuals to enter the next generation through manual selection operators. The combination of genetic operators and local search strategy not only strengthened the overall search ability of the algorithm, but also improved the convergence speed so that the algorithm could find better solutions faster.
In order to verify the effectiveness of the algorithm's improvement, this paper analyzed the performance of the proposed algorithm and several excellent algorithms for instances of different scales. The test results indicated that the proposed method was effective and superior in solving the spherical VRPTW model, and its results were better than those of other algorithms. With the increase in the instance's size, the gap became more obvious. Finally, this paper analyzed the improvement method for this method, and the experimental result showed that the improvement of the proposed algorithm was effective.
However, according to the NFL theorem, the proposed algorithm still has some limitations, such as the search ability of the algorithm not being stable enough and the running time being relatively long. Therefore, the performance of MG-ChOA will continue to be explored and improved through practical applications in the future, and the spherical VRPTW model studied in this paper will also be discussed and studied in combination with green logistics, robot path planning, and other topics.

Conflicts of Interest:
The authors declare no conflict of interest.