Novel Graph Model for Solving Collision-Free Multiple-Vehicle Traveling Salesman Problem Using Ant Colony Optimization

: In this paper, a novel graph model to ﬁgure Collision-Free Multiple Traveling Salesman Problem (CFMTSP) is proposed. In this problem, a group of vehicles start from di ﬀ erent nodes in an undirected graph and must visit each node in the graph, following the well-known Traveling Salesman Problem (TSP) fashion without any collision. This paper’s main objective is to obtain free-collision routes for each vehicle while minimizing the traveling time of the slowest vehicle. This problem can be approached by applying speed to each vehicle, and a novel augmented graph model can perform it. This approach accommodates not only the position of nodes and inter-node distances, but also the speed of all the vehicles is proposed. The proposed augmented graph should be able to be used to perform optimal trajectories, i.e., routes and speeds, for all vehicles. An ant colony optimization (ACO) algorithm is used on the proposed augmented graph. Simulations show that the algorithm can satisfy the main objective. Considered factors, such as limitation of the mission successfulness, i.e., the inter-vehicle arrival time on a node, the number of vehicles, and the numbers of vehicles and edges of the graph are also discussed.


Introduction
Distribution systems involving autonomous vehicles, such as automated guided vehicles (AGV), are an interesting research topic that continuously grows in the operational research area. Routing and scheduling have been dominating issues to explore. Particularly, in applications involving multiple vehicles, vehicle routing problem (VRP) and its variants are well-known approach to solve client services. Such problems are usually generalized by the traveling salesman problem (TSP) [1][2][3][4][5]. Many approaches to solve complexity in TSP are explored, such as min-k-SCCP [6][7][8][9].
Based on our investigation, less of the publications address a problem of asymmetric TSP (ATSP). In this problem, an individual (vehicle, human, etc.) is required to travel from node-to-node on a can visit the same node in an overlapped time. In Reference [38], a variant of TSP, namely Close Enough Traveling Salesmen Problem (CETSP), is presented. This work considers the node more as an area with a predetermined radius rather than the nodes as points. Moreover, the area is allowed to be occupied by only one vehicle at any time. However, similar to Reference [13], the node-to-node travel time in both articles is assumed to be predetermined. This assumption is practically not realistic, because a single (autonomous) vehicle is controlled automatically by a computer-based control system [39].
The control system has to send some speed instructions to the AGV while traveling some paths [35]. Several published algorithms did not consider speed selection, because they was assumed that some constant speeds are applied to the paths [25][26][27].
The organization of this paper is described as follows. Section 2 introduces several terminologies needed to support the analysis and discussions. Section 3 describes the problem statement. Section 4 explains the proposed methods, including the proposed augmented graph and the application of ACO on the graph. Section 5 reveals the simulation results and discussions. Finally, Section 6 concludes the overall work and describes future works.

Preliminaries
Let G(V, E) be an undirected graph, where V = {v i } N v i=1 denotes a set of N v nodes, and E = e i,j is a set of edges connecting two nodes v i and v j , i.e., e i,j = v i , v j i, j ∈ {1, 2, . . . , N v } . In addition, let B = b i,j ∈ {0, 1} N v ×N v be an adjacency matrix of G, i.e., the matrix that describes the connectivity of any pair of nodes in V, where b i,j is assigned to one if v i and v j , i j, are connected and zero elsewhere. Let S = {s l |l ∈ {1, 2, . . . , N s }} be a set of N s speed options s l to be applied to all vehicles at any node in G. Suppose that there exists a group of N u vehicles that are assigned to visit each node in G.

Definition 1.
Route and sub-route. A sub-route from v i to v j is defined as e i,j , and a route is defined as a sequence of sub-routes, i.e., Ξ = e i,j i, j ∈ {1, 2, . . . , N v } that begins from the start node v i, start ∈ V to the end node v i,end ∈ V. The sub-route and route that are traveled by the k-th vehicle, k ∈ {1, . . . , N u }, are denoted as e k i,j ∈ E and Ξ k , respectively.

Definition 2.
Trajectory and sub-trajectory. A sub-trajectory ψ i,j,l is defined as a pair of sub-routes and speed option e i,j , s l , or in other words, the sub-route from the i-th node to the j-th node by applying the l-th speed option. A trajectory is a sequence of sub-trajectories, denoted as Ψ = ψ i,j,l connecting the start node v i,start ∈ V to the end node v i,end ∈ V. The sub-trajectory and trajectory that are traveled by the k-th vehicle are denoted as ψ k i,j,l , and Ψ k , respectively. Definition 3. Arrival time. Arrival time of the k-th vehicle to the i-th node, denoted as t k i , is defined as the time when the vehicle starts to enter the node. Definition 4. Operational time. The operational time of the k-th vehicle on the i-th node, denoted as t op , is defined as the difference between the times the vehicle leaves and enters the node.
i > t op , k 1 k 2 , and during the time interval each vehicle is out in the inevitable collision states defined in References [21,22].

Problem Statement
In this study, we enhanced the problem of the typical TSP (see [6,12]) to collision-free multiple TSP (CFMTSP). In this problem, each vehicle attempts to establish its individual TSP mission. The CFMTSP is described as follows. First, we need to find a complete trajectory for all vehicles, i.e., ψ k i,j,l , such that the following function is minimized: subject to the following: 1.
There is no collision between any vehicles (see Definition 6) at each node.

2.
The start times, t k start , of all vehicles are zero, i.e., t 1 Note that the second constraint is designated to be a collision indicator. If the constraint is violated, then the vehicles have collided with each other. The third constraint emphasizes that there is no delay among the start times of the vehicles. Assumption 1. All vehicles start from different nodes.

Assumption 2.
The operational times, t op , for all nodes in Gare assumed to be constant. Assumption 3. The number of speed options, i.e., N s , is the same for all vehicles.

Assumption 4.
Collisions are considered only at the nodes. The edges are assumed to be sufficiently broad, so that any vehicles passing through the same edge will not collide with each other. Assumption 5. The graph G(V, E) is not a multi-graph, i.e., the number of edges between any two nodes is exactly one.

Graph Model
Solving the problem described in (1) is difficult by using typical graph G, since there is no information about the arrival time of each vehicle at each node. Consequently, the inter-vehicle collision problem is unable to solve. To address such a problem, instead of using typical graph G, we developed a novel structure of the graph, which is called an augmented graph, denoted by G a . The augmented node can be constructed from the typical graph G. Let G a (V a , E a ) be an augmented graph, where V a = ψ i,j,l , as long as b i,j = 1, i, j ∈ {1, . . . , N v } and l ∈ {1, 2, . . . , N s }, is defined as a set of augmented nodes and E a = ξ i 2 ,j 2 ,l 2 i 1 ,j 1 ,l 1 , where ξ i 2 ,j 2 ,l 2 i 1 ,j 1 ,l 1 = ψ i 1 ,j 1 ,l 1 , ψ i 2 ,j 2 ,l 2 is augmented edge, i.e., start and end connected sub-trajectories pairs ψ i 1 ,j 1 ,l 1 and ψ i 2 , j 2 ,l 2 . Note that the notation ξ i 2 ,j 2 ,l 2 i 1 ,j 1 ,l 1 implies that ψ i 2 ,j 2 ,l 2 must be a successor of ψ i 1 ,j 1 ,l 1 . Therefore, it is required that i 2 = j 1 . Figure 1 visualizes the proposed G a , V a , and E a . G a expands the typical G from the node-to-node relation into transition-to-transition relation. In the typical graph G, the edges are weighted by node-to-node distance L i,j while in the augmented graph G a , the augmented edges are weighted by acceleration, whose formulation is conducted using start and final speeds the node-to-node distance (see Equation (5)).   Figure 2 shows that, for a node-to-node trajectory, there are some speed alternatives to apply. Therefore, the augmented edge between any pair of augmented nodes can represent information of acceleration and traveling time, for instance, the transition from , , to , , . The applied initial and target speeds at the endpoint of , , are and , respectively. Therefore, the uniform acceleration , along the augmented edge is formulated as following: where , is the length of the edge e , , and it is plotted to the transition from , , to , , .

Additional Adjacency Matrix
We introduced some supporting matrices to support the proposed algorithm. First of all, we introduced edge matrix, denoted as = , | , ∈ {1,2, … , } , whose dimension is the same with the adjacency matrix . Let , where q is the identifier of an edge whose value can be determined by Algorithm 1, be defined as the edge identifier of each element of , i.e., , that has the value of 1. Therefore, we get the following: where , ∈ {1,2, … , }.  Figure 2 shows that, for a node-to-node trajectory, there are some speed alternatives to apply. Therefore, the augmented edge between any pair of augmented nodes can represent information of acceleration and traveling time, for instance, the transition from ψ 1,2,1 to ψ 2,3,2 . The applied initial and target speeds at the endpoint of ψ 1,2,1 are s 1 and s 2 , respectively. Therefore, the uniform acceleration a 1,2 along the augmented edge is formulated as following: where L i,j is the length of the edge e 1,2 , and it is plotted to the transition from ψ 1,2,1 to ψ 2,3,2 .
Algorithms 2020, 13, x FOR PEER REVIEW 5 of 20  Figure 2 shows that, for a node-to-node trajectory, there are some speed alternatives to apply. Therefore, the augmented edge between any pair of augmented nodes can represent information of acceleration and traveling time, for instance, the transition from , , to , , . The applied initial and target speeds at the endpoint of , , are and , respectively. Therefore, the uniform acceleration , along the augmented edge is formulated as following: where , is the length of the edge e , , and it is plotted to the transition from , , to , , .

Additional Adjacency Matrix
We introduced some supporting matrices to support the proposed algorithm. First of all, we introduced edge matrix, denoted as = , | , ∈ {1,2, … , } , whose dimension is the same with the adjacency matrix . Let , where q is the identifier of an edge whose value can be determined by Algorithm 1, be defined as the edge identifier of each element of , i.e., , that has the value of 1. Therefore, we get the following: where , ∈ {1,2, … , }. Furthermore, the traveling time, t, related to the acceleration in Equation (2) is formulated as follows:

Additional Adjacency Matrix
We introduced some supporting matrices to support the proposed algorithm. First of all, we introduced edge matrix, denoted as e B = e b i,j i, j ∈ {1, 2, . . . , N v } , whose dimension is the same with the adjacency matrix B. Let q e, where q is the identifier of an edge whose value can be determined by Algorithm 1, be defined as the edge identifier of each element of B, i.e., b i,j that has the value of 1. Therefore, we get the following: where i, j ∈ {1, 2, . . . , N v }. For j = 1 to N v 5: If i j and b i,j = 1, then q = q +1 and q e = q. 6: End 7: End 8: Output: all q e .
For instance, suppose that we have the following adjacency matrix of a graph G whose number of nodes is 5, as follows.
Therefore, by using (4), we have the following: .
The functions row ψ (x), start e (x), and end e (x) return the row index belonged to the row in ψ B containing x, the indexes of start and end node in e B, respectively. In the previous example, there are ten edges and three speed options. The number of trajectories is 30 and the dimension of ξ B = Z 30×30 . Until this step, we have the following adjacency matrices: e B, ψ B, and ξ B. Suppose that we are given h ξ; then, by using the map described in ξ B, we obtain p 1 ψ and p 2 ψ, where p 1 and p 2 are the indexes of the row and column of ξ B associated to h ξ. After that, we check the map in ψ B, and we obtain the edge q e and speed option s l for each p 1 ψ and p 2 ψ.

Ant Colony Optimization
Ant Colony Optimization (ACO) was first introduced in Reference [28]. This algorithm is powerful for routing problems such as traveling salesman problem (TSP), vehicle routing problem (VRP), and their variants [17,27]. The algorithm mimics the behavior of a colony of ants in foraging activities. Suppose that, at the beginning, the colony has no information about the location of the food source. An ant system (AS) consists of a set of artificial ants which perform foraging activities, from their nest to a food source. Some ants begin to move randomly to any direction and deposit a chemical called "pheromone", whose trails will be traced by other ants. This process is iterated such that the optimal route is that with the maximum number of pheromone trails.
In this study, we used an ACO algorithm for finding trajectories for a multiple-vehicle system. The algorithm is different from the typical ones. It involves more than one species of ants whose pheromone the others cannot detect trails of. In the cases of CFMTSP, a species represents a vehicle. We develop an algorithm such that each species performs a collision-free trajectory, i.e., all species can prevent collision with each other. In this study, we assume that any two species are not colliding with each other if the difference of their arrival times exceeds a minimum allowable value.
Suppose that there exist N u ant species. Each species consists of N m ants. The r-th ant of the k-th species, denoted by θ k,r , k ∈ {1, 2, . . . , N u }, r ∈ {1, 2, . . . , N m }, represents the k-th vehicle. Each ant is assumed to be able to recognize only the trails produced by ants of the same species. Therefore, if there exists a large concentration of pheromone in a sub-route, if it comes from different species, then it is less possible for the ant to choose that sub-route as a choice. As the typical ACO algorithms, as one ant passes a sub-route, it leaves pheromone trails along the sub-routes. From now on, the others will smell the trails, and based on the largest amount of pheromones, it will choose the sub-route. Even though there are pheromones produced by other different ant species, the ant cannot recognize them. This behavior is similar to the behavior of some colonies of ants in the real world, that is, they cannot recognize the trail of other different ant species, as reported in Reference [40]. Researchers in that study have discovered that a species of ants, i.e., Lasius nigers (La. nigers), is unable to recognize the pheromone trails produced by Novomessor cockerelli (N. cockerelli) and Linepithema humilis (Li. humilis). The reason is that their pheromone trails are constructed by different chemicals.
In this model, we applied the amount of pheromone applied to the proposed augmented graph G a to augmented edges h ξ. Define n, h τ k , and ∆ h τ k,r as the iteration, the total amount of pheromone left by the k-th ant species on augmented edges h ξ, and the increase of pheromone amount left by each r-th ant of the k-th ant species, respectively. We formulate h τ k as follows.
if a solution is found and the i-th and j-th nodes are the part of the augmented edge passed through by the k-th ant species, and where ρ ( 0, 1] is defined as evaporation rate, if the search fails to find a solution. Furthermore, ∆ h τ k,r is formulated as follows: where d k,r is the total distance traveled by the r-th ant from the k-th species, and Q is a positive constant. We apply the probability of selecting h ξ, i.e., Pr( h ξ) as follows: Note that the selection of trajectories of more than one vehicle leads to a consequence of inter-vehicle collision-checking. Therefore, a procedure for checking the collision is developed.
Before describing the main algorithm, we define a number of variables: v k,r i,start and s k,r start are the initial position and speed, respectively. L k unv , L k e , L k ψ , L k ξ , and L k ξ,sel are the lists of unvisited nodes, collection of edges, collection of trajectories, and collection of augmented-edges, the sequence of selected augmented-edges from the start to end nodes, respectively. L i vis is a collection for each node i ∈ {1, 2, . . . , N v } that stores information of the k-th vehicle that has visited the node and its associated arrival time, t k i . Ψ k best is the best trajectory performed by the k-th species. In Line 1 of Algorithm 3, the input is the graph G, initial positions v k,r i,start , and speed s k,r start of each r-ant of the k-th ant species. Note that r also represents the index of iteration. Consequently, N m represents the number of iterations. In Lines 3 and 4, all the required adjacency matrices are constructed, and L k unv for all ant species is set to be empty. The searching process starts from Line 5. For each iteration, trajectories are constructed for each species. In Line 7, some required initializations are performed, such as the initial amount of pheromones. The values are set randomly small to prevent division-by-zero at the beginning at the process. The next processes are focused on identifying the augmented edge of the current occupied node, v k curr . The function CreateEdgeList(v k curr ) is purposed to extract all edges q e in G whose start nodes are v k curr . In this process, the edge matrix e B is used. The resulted q e is then pushed to L k e . The next step is to extract the sub-trajectories, p ψ, whose edges are q e. However, two conditions make the extraction fail. First, it is possible that the end node of the edge q e, i.e., end e ( q e)), was visited. Therefore, the availability of the end node of the edge must be checked (Line 13). Second, even though end e ( q e)) is available, if there exists another ant occupying the node such that the second constraint is violated, then the process is continued to the next edge. This process is revealed in Lines 14-16.
If the end node of the edge q e has not been visited yet and has no collision issue, then the function CreateTrajectoriesList(L k e , q e) is executed (Line 18). The CreateTrajectoriesList(L k e , q e) uses trajectory adjacency matrix ψ B as reference to find the correct trajectories that are spanned by the q-th edge in L k e . If the end node of the edge q e has been visited previously, the process will check the other edges. If there is no edge available, then it can be concluded that a complete trajectory is failed to found. Therefore, the process is continued to new iteration (Lines 29-32). In Line 30, the amount of pheromone trails is reduced by calling ReducePheromoneTrailsAmount(L k ξ , h τ k ), which applies Equation (12).
In Lines 19-21, the function CreateAugmentedEdgesList (L k ψ , p ψ) is used to extract the augmented edges that is spanned by the p-th trajectory in L k ψ . This function uses augmented-edge adjacency matrix ξ B as a reference to find the correct augmented edge that is spanned by the p-th trajectory in L k ψ .
The result is h ξ and is pushed to L k ξ , and we select the augmented edge, whose probability, calculated using Equation (11), is the largest. The selected augmented edge is then pushed into L k ξ,sel (Line 22). In addition, the maximum traveling time, t k i,max , is calculated in Line 23. Here, the final value of t k i,max is the traveling time of the slowest ant species. Moreover, the member of L k unv that is the same to end e ( q e) is removed (Line 24), since the end e ( q e) is selected to be visited.
After all ants from all species complete their routes, the pheromone trails on all augmented edges are updated in Line 31 by using Equation (11). By using the current pheromone trails, the probability of all augmented edges in L k ξ , k ∈ {1, 2, . . . , N u } is calculated by using CalculateProbability(L k ξ,sel , h τ k ) in Line 37. The global minimum, t k best , is then calculated in Line 38. Finally, the output is ψ k best with t k i,max = t k best , for all k ∈ {1, 2, . . . , N u }.

Results and Discussions
We tested the performance of the proposed algorithm by conducting simulations involving graphs with various numbers of nodes, i.e., 10, 15, and 20 nodes, as shown in Figures 3-5, respectively. For each graph, four variations of connectivity were simulated to show the success of finding solutions. We denoted the graph, together with the connectivity of variations, by using the code "CONFIG A-B", where A is the number of nodes and B is the label of a variation. For instance, the first variation of the graphs with 10, 15, and 20 nodes were denoted as CONFIG 10-1, CONFIG 15-1, and CONFIG 20-1, respectively. Note that all edges were bi-directional.
Three significant aspects were evaluated. The first aspect was the influence of minimum allowable arrival time on the solution's existence and the convergence of the solutions. The second was the correlation between the average degrees of all nodes, the number of vehicles, and the number of nodes to the success of finding solutions. The last aspect was the accuracy of the resulted minimum traveling time of the slowest vehicle, according to the variation of the evaporating constant and its effect on the convergence of the search results.
For evaluating the first and third aspects, we established simulations for three vehicles, where the 1st, 2nd, and 3rd vehicles started from Nodes 4, 1, and 3, respectively. Each vehicle had four speed options, i.e., 0.1, 0.5, 1, and 1.5 m/s. Meanwhile, for evaluating the second aspect, we applied two until seven vehicles on each connectivity configuration. Vehicles 1 until 7 were started from Nodes 4, 1, 3, 5, 10, 7, and 8, respectively. For one simulation, we applied 3000 iterations of searching the optimal solution.
Algorithms 2020, 13, x FOR PEER REVIEW 10 of 20 final value of , is the traveling time of the slowest ant species. Moreover, the member of that is the same to ( ) is removed (Line 24), since the ( ) is selected to be visited.
After all ants from all species complete their routes, the pheromone trails on all augmented edges are updated in Line 31 by using Equation (11). By using the current pheromone trails, the probability of all augmented edges in , ∈ {1,2, … , } is calculated by using CalculateProbability( , , ) in Line 37. The global minimum, , is then calculated in Line 38. Finally, the output is with t , = , for all ∈ {1,2, … , }.

Results and Discussions
We tested the performance of the proposed algorithm by conducting simulations involving graphs with various numbers of nodes, i.e., 10, 15, and 20 nodes, as shown in Figure 3-5, respectively. For each graph, four variations of connectivity were simulated to show the success of finding solutions. We denoted the graph, together with the connectivity of variations, by using the code "CONFIG A-B", where A is the number of nodes and B is the label of a variation. For instance, the first variation of the graphs with 10, 15, and 20 nodes were denoted as CONFIG 10-1, CONFIG 15-1, and CONFIG 20-1, respectively. Note that all edges were bi-directional.
Three significant aspects were evaluated. The first aspect was the influence of minimum allowable arrival time on the solution's existence and the convergence of the solutions. The second was the correlation between the average degrees of all nodes, the number of vehicles, and the number of nodes to the success of finding solutions. The last aspect was the accuracy of the resulted minimum traveling time of the slowest vehicle, according to the variation of the evaporating constant and its effect on the convergence of the search results.
For evaluating the first and third aspects, we established simulations for three vehicles, where the 1st, 2nd, and 3rd vehicles started from Nodes 4, 1, and 3, respectively. Each vehicle had four speed options, i.e., 0.1, 0.5, 1, and 1.5 m/s. Meanwhile, for evaluating the second aspect, we applied two until seven vehicles on each connectivity configuration. Vehicles 1 until 7 were started from Nodes 4, 1, 3, 5, 10, 7, and 8, respectively. For one simulation, we applied 3000 iterations of searching the optimal solution.

The Effect of Minimum Allowable Arrival Time Difference
The first simulation set was established to evaluate the relation between the minimum allowable arrival time difference, , and the solution's existence. The purpose of the simulations is to confirm the hypotheses. In such simulations, we used = 1 and = 1 and evaporate rate = 0.95. We applied five values of , i.e., 10, 50, 100, 150, and 300 s. Table 1; Table 2 show the results of 10 trials for varying allowable minimum occupation time and varying evaporate constant , respectively. In each table, as shown in Column (1), for each occupation time, we established five trials. Column (1) reveals the values of t . Column (2) is the success rate of finding a complete solution. Columns (3) and (4) show the average and standard

The Effect of Minimum Allowable Arrival Time Difference
The first simulation set was established to evaluate the relation between the minimum allowable arrival time difference, , and the solution's existence. The purpose of the simulations is to confirm the hypotheses. In such simulations, we used = 1 and = 1 and evaporate rate = 0.95. We applied five values of , i.e., 10, 50, 100, 150, and 300 s. Table 1; Table 2 show the results of 10 trials for varying allowable minimum occupation time and varying evaporate constant , respectively. In each table, as shown in Column (1), for each occupation time, we established five trials. Column (1) reveals the values of t . Column (2) is the success rate of finding a complete solution. Columns (3) and (4) show the average and standard

The Effect of Minimum Allowable Arrival Time Difference t occ
The first simulation set was established to evaluate the relation between the minimum allowable arrival time difference, t occ , and the solution's existence. The purpose of the simulations is to confirm the hypotheses. In such simulations, we used α = 1 and β = 1 and evaporate rate ρ = 0.95. We applied five values of t occ , i.e., 10, 50, 100, 150, and 300 s. Tables 1 and 2 show the results of 10 trials for varying allowable minimum occupation time t occ and varying evaporate constant ρ, respectively. In each table, as shown in Column (1), for each occupation time, we established five trials. Column (1) reveals the values of t occ . Column (2) is the success rate of finding a complete solution. Columns (3) and (4) show the average and standard deviation of the maximum traveling times, respectively. Column (5) shows the minimum value of maximal traveling times. The application of t occ is evaluated for its effect on the percentage of successful trials from 10 trials, i.e., n success , average traveling time of the slowest vehicle, T avg , the standard deviation of T avg , i.e., σ T , and the minimum traveling time that is ever found from the ten trials, i.e., T min . For this purpose, we use a statistical correlation technique. The correlations between t occ and n success , T avg , σ T , and T min are −0.99, 0.83, 0.75, and 0.88, respectively. In can be concluded that there is a strong negative relationship between t occ , = and n success and almost-strong positive relation between t occ and the other three variables. In addition, the tendency of failure can be identified by a parameter σ T /t occ . We found from the simulation that σ T /t occ is directly proportional to the probability of success in finding a solution, i.e., n success .

Successfulness of Finding a Solution
Since the graph is not fully connected, it is possible that there exists a situation such that any single solution is failed to find as the effect of the number of nodes, the average degree of all nodes increase, and the number of vehicles. The degree of a node is defined as the number of edges that connect to the node. We established simulations for 10, 15, and 20 nodes under t occ = 10 s.

Simulations for 10 Nodes
For the ten nodes cases, we applied four different configurations of connectivity, as shown in Table 3, i.e., CONFIGs 10-1, 10-2, 10-3, and 10-4. CONFIG 10-1 is the graph with the connectivity visualized in Figure 3. CONFIG 10-2 is constructed from CONFIG 10-1 with the elimination of the edge connecting Nodes 5 to 8. CONFIG 10-3 is constructed from CONFIG 10-3, with the elimination of the edge connecting Nodes 5 to 9. CONFIG 10-4 is constructed from the CONFIG 10-1, with the addition of the edges connecting Nodes 1 to 4, 3 to 6, and 4 to 7. Table 3 describes the degree of each node. The averages of the degree of nodes for CONFIG 10-1, 10-2, 10-3, and 10-4 are 3.20, 3.00, 2.80, and 4.00.
We investigate the relationship between the nodes' average degrees and the number of vehicles to the possibility of successfulness of finding a solution. Here, we introduce a parameter φ defined as follows: where r avg is defined as the average degree of all nodes. From Table 3, it can be seen that the solution can be found for φ ≤ 0.19. It can be concluded that, from the 10 trials, at the range of the φ, there exists at least one successful trial. Table 3. Simulation results for 10 nodes. Node  CONFIG 10-1  CONFIG 10-2  CONFIG 10-3  CONFIG 10-4   1  3  3  3  4  2  4  4  4  4  3  3  3  3  4  4  2  2  2  4  5  5  4  3  5  6  3  3  3  4  7  3  3  3  5  8  3  2  2  3

Simulations for 15 Nodes
In this simulation, we added five nodes, as follows. Nodes 11,12,13,14,and 15 Figure 4. CONFIG 15-2 is constructed by adding an edge connecting Nodes 1 and 5 on the CONFIG 15-1. The CONFIG 15-3 configuration is performed by adding edges connecting Nodes 1 to 4, 2 to 12, and 7 to 10 to CONFIG 15-2. The CONFIG 15-4 configuration is performed by adding edges connecting nodes 4 to 7 and 3 to 6 to CONFIG 15-3. The degree of each edge in each connectivity configuration is described in Table 4.  3  3  3  3  12  2  2  3  3  13  3  3  3  3  14  2  3  3  3  15  3  3  3  3 Average 3.60 3.73 4.13 4.40 From Table 4, it can be seen that the solution can be found for φ ≤ 0.08, which means that, at the range of the φ, it can be concluded that, from the 10 trials, there always exists at least one successful trial.

Simulations for 20 Nodes
In this simulation, we added five nodes from the CONFIG 15-1, as follows. Nodes 16,17,18,19,and 20 Figure 5. CONFIG 20-2 is constructed by modifying CONFIG 20-1, i.e., adding edges connecting Nodes 1 to 16, 4 to 16, 2 to 20, 6 to 9, 8 to 19, and 12 to 17. CONFIG 20-3 is performed by adding Edges 2 to 17, 3 to 5 to 9, 7 to 15, and 9 to 15 and eliminating Edges 4 to 16, 4 to 10, and 14 to 17 on CONFIG 20-2. CONFIG 20-4 is performed by the addition of Edges 3 to 17, 4 to 16, 4 to 10, and 5 to 14 to CONFIG 20-3. Table 5 reveals the simulation results of such connectivity configurations. It can be shown that the solution can be found for φ ≤ 0.02. It means that, from the 10 trials, it can be concluded that there always exists at least one successful trial at the range of the φ.

Analysis of Accuracy
From the simulations established in Section 5.3, we can analyze the accuracy of the result, i.e., minimum traveling time of the slowest vehicle. The standard deviation σ T can identify the accuracy of the simulation. As shown in Tables 1 and 2, the averages of σ T for t occ > 10 s and t occ > 150 s are 373.48 s and 631.2 s. It indicates that the results in lower t occ are more accurate than the higher ones.
In addition, we analyzed the relationship between the average T avg and σ T by using the same technique as the one used in Section 5.2., i.e., statistical correlation. The results are described as follows. The correlations between T avg and σ T for t occ > 10 s and t occ > 150 s are 0.92 and 0.84, respectively. It indicates that there is a relatively strong positive relation between T avg and σ T : larger T avg tends to lead to a larger σ T , or, in other words, a larger T avg tends to lead to lower accuracy. Figures 6 and 7 show the searching progression for t occ > 10 s and t occ > 150 s, respectively, each with ρ = 0.1 and ρ = 0.5. It can be concluded that the algorithm's ability to converge is better for t occ > 10 s than for t occ > 150 s. It can be concluded from the final values until the 3000th iteration. Moreover, Figures 6 and 7 confirm the accuracy conclusions conducted from the statistical correlation analysis. Accuracy analysis can infer other behavior, that is, the ability to reach T avg . Figure 7 shows that, for t occ > 10 s, the algorithm's ability to reach T avg is better than that for t occ > 150 s. Furthermore, it can be seen that, for t occ > 150 s, even until the 1500th iteration, the algorithm still can update the minimum of maximum times. However, it cannot go further to T avg . The interesting situation occurs for ρ = 0.5 in the same t occ . In this situation, mostly, once a solution is found, the algorithm cannot find another better solution. iteration, the algorithm still can update the minimum of maximum times. However, it cannot go further to . The interesting situation occurs for = 0.5 in the same . In this situation, mostly, once a solution is found, the algorithm cannot find another better solution.

The Near-Optimal Trajectories
We show two samples of near-optimal trajectories obtained from the simulations to verify that any point in the graph G is visited by different vehicles with the arrival time difference between any vehicles exceeds t . Tables 6 and 7 show the optimal trajectories for > 10 s and > 150 s, respectively. To verify that the arrival time differences exceed , we provide Table 8. According to

The Near-Optimal Trajectories
We show two samples of near-optimal trajectories obtained from the simulations to verify that any point in the graph G is visited by different vehicles with the arrival time difference between any vehicles exceeds t occ . Tables 6 and 7 show the optimal trajectories for t occ > 10 s and t occ > 150 s, respectively. To verify that the arrival time differences exceed t occ , we provide Table 8. According to Table 8, for t occ > 10 s, the smallest arrival time difference is 11.7 s, i.e., between Vehicles 2 and 3 at Node 2. For t occ > 150 s, the smallest arrival time difference is 165.4 s, i.e., between Vehicles 2 and 3 at Node 2.

Conclusions
A method to solve the Collision-Free Multiple Traveling Salesman Problem (CFMTSP) applied to multiple agents was proposed. In this problem, each agent must visit all nodes in a provided graph, following the provided edges. The graph is modified into an augmented graph such that the information of nodes' position and speed options can be accommodated, and in turn, the arrival time to each node can be determined. According to the optimization of collision-free for each vehicle, an ACO involving multiple ant species is utilized on the augmented graph. The pheromone trails are not left on the edges (such as in typical graph model), but the augmented edges. As a consequence, the probability of selection is assigned to those augmented edges, as well. As a result, the solution is not the sequence of routes but also the trajectory (routes and selected speeds). The algorithm has to guarantee that the resulted trajectories are collision-free.
Simulations were established with no violation of minimum allowable arrival time difference. In addition, the simulations have shown three behaviors of the resulted solutions. First, the increase of minimum allowable arrival time difference leads to the decrease of the number of successful trials, almost-strongly increases of the average, the standard deviation, and the ever-found minimum traveling times. Second, it can be concluded that there exists a threshold of the ratio φ such that the solution for all vehicles is unable to find. It was demonstrated that the threshold of φ is lower as the number of nodes increases. Moreover, for the same number of nodes, the increase of the average of nodes' degree leads to the number of vehicles that are able to find solutions successfully. Third, it is concluded that that the increase in the average of minimum traveling time tends to decrease the algorithm's accuracy. Third, the evaporate rate in the algorithm has a weak influence on the average, the standard deviation, and the ever-found minimum traveling times.
Future works will focus on some issues, such as the continuation of speed options and the consequences of modifying the augmented graph. Another issue is in the performance of the ACO algorithm, such as the convergence of the optimal trajectories, the decrease of the standard deviation of the slowest agent's traveling time in high minimum allowable arrival time difference.
Author Contributions: A.K.P. contributed to writing the original draft, formal analysis, conceptualization, methodology, and visualization; D.B.S. contributed to validation, review and editing of the manuscript, and supervision; conceptualization, methodology, and other aspects were agreed on by all authors. All authors have read and agreed to the published version of the manuscript.
Funding: This research received internal funding from Universitas Atma Jaya Yogyakarta, Indonesia.