WiSM: Windowing Surrogate Model for Evaluation of Curvature-Constrained Tours with Dubins vehicle

by   Jan Drchal, et al.
Czech Technical University in Prague

Dubins tours represent a solution of the Dubins Traveling Salesman Problem (DTSP) that is a variant of the optimization routing problem to determine a curvature-constrained shortest path to visit a set of locations such that the path is feasible for Dubins vehicle, which moves only forward and has a limited turning radius. The DTSP combines the NP-hard combinatorial optimization to determine the optimal sequence of visits to the locations, as in the regular TSP, with the continuous optimization of the heading angles at the locations, where the optimal heading values depend on the sequence of visits and vice versa. We address the computationally challenging DTSP by fast evaluation of the sequence of visits by the proposed Windowing Surrogate Model (WiSM) which estimates the length of the optimal Dubins path connecting a sequence of locations in a Dubins tour. The estimation is sped up by a regression model trained using close to optimum solutions of small Dubins tours that are generalized for large-scale instances of the addressed DTSP utilizing the sliding window technique and a cache for already computed results. The reported results support that the proposed WiSM enables a fast convergence of a relatively simple evolutionary algorithm to high-quality solutions of the DTSP. We show that with an increasing number of locations, our algorithm scales significantly better than other state-of-the-art DTSP solvers.


Path Planning Algorithms for a Car-Like Robot visiting a set of Waypoints with Field of View Constraints

This article considers two variants of a shortest path problem for a car...

Lifelong Multi-Agent Path Finding in Large-Scale Warehouses

Multi-Agent Path Finding (MAPF) is the problem of moving a team of agent...

Using Recursive KMeans and Dijkstra Algorithm to Solve CVRP

Capacitated vehicle routing problem (CVRP) is being one of the most comm...

The Bounded Acceleration Shortest Path problem: complexity and solution algorithms

The purpose of this work is to introduce and characterize the Bounded Ac...

Capacitated Vehicle Routing with Target Geometric Constraints

We investigate the capacitated vehicle routing problem (CVRP) under a ro...

Improved Tabu Search Heuristics for Static Dial-A-Ride Problem: Faster and Better Convergence

Multi-vehicle routing has become increasingly important with the rapid d...

I Introduction

In this paper, we address the Dubins Traveling Salesman Problem (DTSP) by the proposed regression surrogate model for fast evaluation of possible solution candidates given by a sequence of visits to the locations of interest by speeding up the convergence of an evolutionary solver that supports finding high-quality solutions. The addressed DTSP is motivated by data collection and surveillance missions where a vehicle is requested to visit a set of locations as quickly as possible [27, 8]. The problem is a variant of the routing problem [26] to determine the optimal sequence of visits to the given set of locations such that the length of the curvature-constrained path connecting the locations in the sequence is minimal [31]. The curvature-constrained path is requested because of motion constraints of the vehicle that are modeled by Dubins vehicle [7] which moves only forward with the limited turning radius . A solution of the DTSP consists of sequencing part to determine the optimal sequence of visits to the locations and continuous optimization to determine the optimal heading values. The DTSP is at least NP-hard [21] because for , it becomes the Euclidean TSP (ETSP) [1].

In 1957, Dubins showed that for two given locations with the prescribed leaving and arrival heading angles of the vehicle, the optimal curvature-constrained path connecting the locations is one of six possible maneuvers that can be found by a closed-form solution [7]. The generalization to an arbitrary number of locations is known as the Dubins Touring Problem (DTP). The state-of-the-art DTP solvers are based on a sampling of the possible heading angles [25]. In [9], an informed sampling method, which we further denote the Iteratively-Refined Informed Sampling (IRIS), enables finding a solution of the DTP close to optimum, i.e., with the relative optimality gap to the lower bound less than 1%, in tens of seconds for problems with up to 100 locations, which outperforms the uniform sampling [25].

The DTSP solvers reported in the literature can be roughly classified into three groups according to 

[23, 33]. The first group represents decoupled

approaches, in which the sequence of visits is determined independently from the headings, e.g., using a solution of the related ETSP. Then, headings are determined in the second step using a heuristic Alternating Algorithm (AA) 

[31], the Local Iterative Optimization (LIO) [33], or by the IRIS proposed in [9]. The methods of the second group involve a joint optimization of both the sequence and the headings using Evolutionary Algorithm (EA) as it is reported in [13, 36]. The last group includes approaches based on transformation of the original problem into a purely combinatorial routing problem. In [27] and [28], the heading angles are sampled, and the discretized instance is transformed to the Generalized TSP that is further transformed into the Asymmetric TSP that can be solved by the heuristic LKH solver [22], or Concorde solver [5] for the optimal solution with respect to the selected sampling. Similar sampling-based strategy is utilized to find a tight lower bound of the DTSP in [24].

The existing decoupled approaches such as [25, 9] use only a single sequence found as a solution of the Euclidean TSP without considering the minimum turning radius of the Dubins vehicle which limits the quality of the DTSP solutions for cases where cannot be ignored (i.e., dense location distribution w.r.t. ). The joint optimization and transformation methods may provide better solutions at the cost of low scalability (dense sampling leads to large problem instances of the transformed problems) and very high computational requirements because a large search-space needs to be covered by the global optimizer such as the evolutionary algorithm.

The herein presented approach is a novel DTSP solver that combines features of both the decoupled approaches and joint optimization. The main idea of the proposed method is that a relatively simple EA optimizes the sequence of visits separately from the continuous optimization part, and the sequence quality (i.e., fitness) of each individual is assessed by means of length of the corresponding DTP solution (which also determines the heading angles). The information about headings is not stored in the genome itself, but rather computed during every evaluation of its fitness. It reduces the complexity of the search-space and makes it strictly combinatorial as the sequence of visits is a permutation; however, it also makes the evaluation of each individual more expensive. Thus, the critical part of our method is the need to evaluate possibly many candidate sequences. Even though the solution of the DTP by the IRIS [9] is significantly faster than the uniform sampling [25]; it is still too computationally demanding.

Similarly to [35], where the authors speed up the expensive evaluation of the objective function using a surrogate model, we propose the Windowing Surrogate Model (WiSM) leveraging on the high-quality solutions of the IRIS-based DTP solver. The main contribution of the proposed WiSM is considered in the windowing decomposition developed from a lower bound which allows training of the model on small fixed-size DTP instances, while the evaluation of the solution length can be performed on DTP instances of arbitrary size. It makes WiSM easily implementable by most types of regression models.

Based on the presented results, WiSM provides solutions competitive to a very dense sampling of the heading values in the DTP, but with significantly lower computational requirements. The paper is mainly focused on the surrogate model based on the windowing decomposition and not on details and tuning of the sequence optimization; however, we show that the combination of WiSM with a simple EA (further denoted as WiSM-EA) outperforms other DTSP solvers for a wide range of instance sizes.

The paper is organized as follows. The addressed DTSP is defined in the next section. A brief overview of the IRIS, the high-quality DTP solver used in learning of the proposed WiSM, is presented in Section III. The model is employed in the Evolutionary Algorithm summarized in Section IV. The surrogate model itself with the windowing decomposition is proposed in Section V. Empirical evaluation of the WiSM-EA and its comparison with existing approaches for the DTSP are in Section VI. Concluding remarks are given in Section VII.

Ii Problem Statement

The proposed solution of the DTSP is based on an evaluation of the candidate sequences to visit the given set of target locations , using a solution of the DTP where the path connecting the locations has to respect the motion constraints of Dubins vehicle [7] that is moving only forward with a constant forward velocity and a minimum turning radius . The state of the vehicle can be expressed as , where is the position of the vehicle and is its heading angle. Thus the state is from the Special Euclidean group, . The motion of the vehicle can be described as


where is the bounded control input. For simplicity and w.l.o.g., we further consider and .

A solution of the DTSP can be expressed as a permutation for that defines a sequence of visits to the locations and particular heading values corresponding to each particular location.

The DTSP stands to find a sequence of visits to with the respective heading values such that the length of the Dubins path connecting the locations is minimal. The problem can be considered as combinatorial optimization over all possible sequences and -variable continuous optimization of the heading values .

Problem 1 (Dubins Traveling Salesman Problem - DTSP)

where denotes the analytically computed length of the optimal Dubins maneuver [7] between and , and for is defined to simplify the notation in (2).

Finding a solution of the DTSP is an optimization problem with a combinatorial part over and continuous part over . If is given, the problem becomes strictly continuous optimization with variables in . The problem is then called the Dubins Touring Problem (DTP), and the solution cost for a particular sequence is denoted .

Problem 2 (Dubins Touring Problem - DTP)

Due to the combinatorial nature of the TSP, many candidate sequences need to be evaluated. Thus, the evaluation of (3) should be fast enough to provide competitive results to existing DTSP approaches [23]. In Section V, we propose a surrogate approximator of (3) trained using high-quality solutions of the DTP found by the IRIS method [9] briefly described in the following section.

Iii Background – Dubins Touring Problem (DTP)

A solution of the DTP is needed to evaluate the cost (3), and thus the final solution of the DTSP, but the solution is also needed to train the proposed surrogate model for fast evaluation of possible candidate sequences. Finding optimal heading angles for the DTP with dense locations is a challenging problem because the length of the Dubins tour connecting a sequence of locations is a piecewise continuous function [12]. Sampling-based approaches are thus utilized to sample the domains of the continuous variables into finite sets of discrete values of possible heading angles. A solution can be then found as the shortest path in a graph representing the discretized instance of the DTP as follows.

Let have samples of the heading angle per each location which form possible states per each . Then each pair of the states corresponding to two consecutive locations in the sequence can be connected by the optimal curvature-constrained path determined as the optimal Dubins maneuver [7]. Having the states and their connections, a search graph can be created, and the shortest path for the given sequence of locations and defined sampling can be found in by a forward search based on dynamic programming [9]. Such a path is a feasible solution to the DTP.

High-quality solutions, however, require dense sampling [25], which can be computationally demanding. The IRIS method [9] decreases the computational requirements by an iterative refinement of the heading intervals based on the tight lower bound solution from [25].

The DTSP solver presented in this paper requires a high number of sequence evaluations, which makes even the IRIS approach [9] too computationally demanding. Therefore, we propose to employ a surrogate model for fast evaluation of Dubins tour lengths to speed up the convergence towards high-quality DTSP solutions using generated sequences by the Evolutionary Algorithm described in the following section.

Iv Evolutionary Algorithm for the DTSP

In this section, we present the Evolutionary Algorithm (EA) for solving the DTSP. Contrary to existing evolutionary approaches to the DTSP, such as [36] where sequences and headings are encoded together, the proposed evolutionary approach is utilized only for generating proper sequences , and the heading values are determined as a solution of the corresponding DTP. The particular cost function can be determined by IRIS for the DTP [9], but also by its surrogate model estimation proposed in Section V. A generational EA with a simple elitism, where an individual of the population is represented by a permutation of target locations (path representation [20]) is used, and it is summarized in Algorithm 1.

Input :  – a set of the given target locations; – population size; – tournament size;

– mutation probability;

– elite size; S -- termination time;
Output : DTSP solution – (, , )
1 initializePopulation(, );
2 evaluatePopulation(, );
3 while CPU  do
4       ;
5       for  to  do
6             tournamentSelection(, );
7             if random()  then
8                   mutation();
10            else
11                   tournamentSelection(, );
12                   crossover(, );
14            ;
16      evaluatePopulation(, );
17       best();
18       ;
19       for  to  do  ;
20       sort();
21       for  to  do  ;
24 ;
25 return ;
Algorithm 1 Evolutionary Algorithm for the DTSP

The initial population is filled by random permutations in the initializePopulation function representing sequences of visits to the given locations. The fitness of each individual in the population is evaluated using the solution (or its estimation) of the corresponding DTP by the evaluatePopulation function. The main loop (Line 1, Algorithm 1) iterates over the generations until the overall running time CPU reaches the termination time S.

The evolution is performed as follows. A new population of individuals is generated either by mutation with the probability of , or crossover otherwise (Lines 11, Algorithm 1). The mutation operator is the Simple Inversion Mutation (SIM) [15, 20] which reverses a random sub-sequence. The crossover method implements the well-known Order crossover (OX) [6] which copies a random sub-sequence of the first parent and then adds locations from the other parent preserving their order but leaving out those already used. We specifically use the OX1 variety as described in [20], although we generate only a single offspring. Both mutation and crossover employ tournament selection with the tournament size (tournamentSelection) to select the parent(s).

The new population is then evaluated (Line 1, Algorithm 1) and the best solution of the both current and new populations is selected as . Finally, least viable individuals of are replaced by the copies of the best so far solution introducing the elitism (Lines 11, Algorithm 1).

The algorithm terminates when the stop condition is met. The best sequence is extracted from the population (Line 1, Algorithm 1) and the corresponding headings are determined (Line 1, Algorithm 1) calling the DTP solver [9] to provide a feasible solution. On the other hand, sampling-based approaches would be too computationally demanding for evaluation of of all population individuals in evaluatePopulation, and therefore, we propose a surrogate model to get a considerable speedup.

V Proposed Windowing Surrogate Model (WiSM) for a Fast Estimation of the DTP Solution Cost

The proposed Windowing Surrogate Model (WiSM) is designed to approximate with a surrogate function to significantly speed up the evaluation of the DTP instance costs. The WiSM is evaluated based on overlapping windows in a convolutional manner, and the partial results are averaged to get the final cost estimation.

In the studied DTSP, the number of locations

can vary depending on the particular problem instance. Thus, we need to address the limitation of the regression models (e.g., neural networks, random forests, etc.) that are mostly limited to fixed sized inputs. While this issue can be directly approached by recurrent neural networks that allow processing inputs of variable-length as sequences (e.g., LSTM 

[14] or GRU [4]), a collection of training data would be non-trivial as we would need to collect examples of DTP instances having a varying number of locations with a little guarantee on how the network would generalize on sequences of the unseen lengths. To overcome the problem of the variable-size input, we propose a decomposition of the cost function to fixed-size subproblems based on the sliding window technique.

V-a Sliding Windowing based Cost Estimation of the DTP

The proposed idea is to evaluate the cost of Dubins tour for small sub-sequences of the target locations limited by the window size . The total tour cost can be then estimated as an aggregate of the particular costs.

For a specific sequence we can define the particular cost of the subtour of Dubins optimal tour for the fixed -size window as

Fig. 1: Sliding window principle with overlapping windows of the size . Each window covers 2 trips, i.e., 3 locations.

We propose to utilize a sliding window that starts at each target location as it enables to address arbitrarily-length instances for which the number of locations is not necessarily divisible by the utilized window size . The whole Dubins tour is divided into sub-sequences of the windows for each connecting successive locations as it is depicted in Fig. 1. Then the optimal cost for the specific and known optimal heading angles can be computed from the window costs exactly as


The optimal heading angles are apparently not available, and thus (5) provides only an intuition behind the proposed approach for fast cost estimation of the Dubins tour.

The cost of the specific window is estimated independently on the final solution and finding the cost of the window sub-sequence is considered as the Open DTP [9], which can be defined as the continuous optimization problem


The Open DTP is further denoted as to emphasize it does not involve a closed tour and the heading angles at both end locations are unconstrained. Note that the unlike in (2), the length of the closing Dubins maneuver is omitted from the definition:


Comparing the particular DTP cost of the -th window defined by a sub-sequence and the corresponding cost , it can be shown that


The proof by contradiction is based on the fact that is a relaxed version of . More specifically, unlike for , the boundary angles and of are not constrained. Thus and can take advantage of arbitrarily chosen heading angles at the endpoints.

Fig. 2: A solution of the closed DTP (dashed line) with one of open DTP subtours (thick solid line).

An example of the Closed DTP with locations with a corresponding instance of the Open DTP over four locations , and is depicted in Fig. 2. Note that the optimal tours between and differ for the open and closed cases. The open version is shorter which corresponds with (8).

The cost of the whole sequence can be estimated as using windows with the size as


Based on (8), the estimation is a lower bound of the optimal cost and we get


In the following, we expect the lower bound (9) is tight enough even for relatively small size and the minimization w.r.t. can be employed as a proxy to the direct minimization of . Thus, we can evaluate DTP costs for instances of a fixed size and use (9) to get the cost estimate for an instance of arbitrary size .

The computational complexity of approximation employing the sampling methods such as IRIS [9] can be bounded by , where is the number of heading values per each location. In practice, the window size is expected to be constant and significantly lower than both and ( and ). Thus, the computational complexity can be bounded by , which is a noticeable improvement in comparison to for the evaluation of . While the performance improvement is appealing, it turns out, that in practice, the sampling methods are still too slow to allow for efficient evolutionary search. Nevertheless, the decomposition, based on a repeated evaluation of fixed-sized sub-tours, enables approximation of the relatively slow sampling-based methods by a wide selection of regression models.

The utilized sliding window approach is a commonly used general technique to process variable length data; however, to the best of the authors’ knowledge, it has not been used for tour length estimates yet. In the particular case of the proposed WiSM, the specific combination of the sliding window with the selected aggregation scheme gave a clear interpretation (10), which is only possible by incorporating domain-specific knowledge of the addressed DTSP using results of the DTP.

V-B Windowing Surrogate Model (WiSM) Approximation of the Cost Estimation

The key idea of the proposed WiSM-based approximation of the DTP solution cost is based on a fast evaluation of the partial costs of by a surrogate regression model. In existing deployments of the surrogate models in evolutionary algorithms, the model is a proxy to the real fitness that is trained online while EA runs. The switching between the real and surrogate fitness can be realized by more or less complex model management [17, 18]. On the other hand, the real-time response of the proposed DTSP solver can benefit from pre-trained surrogate model [16], and therefore, the model is learned offline using a large training dataset based on high-quality solutions of the open-loop .

The WiSM can be implemented using virtually any type of regression model, but we consider the Multi-Layer Perceptron (MLP) 


in this paper. It is because the feed-forward neural network supports fast evaluation of all partial costs

in a single batch. Thus, it can take advantage of highly optimized matrix multiplication routines111We can evaluate the whole population of the evolutionary algorithm in a single batch, which has been utilized for the herein reported empirical results..

Each input vector fed to the surrogate model is constructed as

while the output constitutes single real value corresponding to . Notice, the number of locations is because the cost of Dubins tour is defined for two locations and more, see (4). Hence, the MLP architecture has

inputs, several hidden layers of non-linear units such as ReLU 


, and a single linear output neuron which is common for regression tasks.

Vi Results

The proposed Windowing Surrogate Model (WiSM) has been evaluated in combination with the evolutionary algorithm presented in Section IV and the combined method is further denoted WiSM-EA. The solved DTSP instances are of various sizes and different densities of the locations.

The performance of the WiSM-EA is compared with the existing DTSP approaches. Namely, the Alternating Algorithm (AA) [31], the Local Iterative Optimization (LIO) [33], and the Sampling-based Algorithm (SA) [29]. The SA transforms the problem to the Generalized TSP that is further transformed into Asymmetric TSP solved by the LKH solver [22].

As a baseline, we take the decoupled approach from [9] where the DTSP solution takes the location sequence from the solution of the corresponding Euclidean TSP (obtained by the Concorde solver [5]). The headings are consequently determined by the IRIS method. Besides, we also initially considered the Memetic algorithm [36]; however, the achieved results are not competitive with the selected approaches regarding the solution quality and computational requirements. Therefore it is not included in the reported results.

The DTSP methods were compared using instances, randomly generated similarly to [9]. The benchmark set consists of 10 instances per each possible pair , where denotes the number of target locations and is the density of the locations. The targets of each instance were uniformly sampled from a square region with a side , i.e., .

Due to unknown optimal solutions of the DTSP instances, the results are evaluated using the normalized cost defined as


where is the cost provided by the baseline.

The headings for the baseline and also for the WiSM-EA (Line 1, Algorithm 1) were determined using IRIS [9] with the maximum number of samples set to .

The parameters used by the evolutionary algorithm were: the population size , tournament size , mutation probability , and elitist size .

The WiSM-EA was implemented in Julia language [2], run with -O3 and --check-bounds=no options and the computational requirements were further significantly decreased by caching model approximations using as a key. The IRIS method (called from Line 1, Algorithm 1) and also all the other algorithms (AA, LIO, SA) were implemented in C++ and compiled by the gcc compiler with -O3 and -march=native flags, which also holds for the LKH. In all the cases, the algorithms were run on a single core of the Intel Xeon Gold 6130 @ processor.

Vi-a Learning Setup of the Proposed WiSM

The choice of an appropriate window size of the WiSM deals is a trade-off between the model complexity, accuracy, and computational time needed to evaluate the model. Increased window size leads to the increased number of the model inputs and parameters that also increases demands on the size of the training dataset. Besides, a prediction for a complex model is computationally demanding which might be critical for the evolutionary search because it can limit the total number of evaluations achievable under the particular computational time limit S (Line 1, Algorithm 1).

Empirical evaluation has been performed for in the same way as described in this section, but data are not shown for brevity. Based on that, we found out that provides overall best results, and therefore, this window size is considered for the rest of the presented results. The achieved average computational time per a single window for is and with and without the caching, respectively.

The empirically selected neural network architecture was the MLP as described in Section V-B having three hidden layers of 256 ReLU units each. The training dataset consisted of samples of coordinates and the corresponding target cost value was computed using a solution of the open DTP found by the IRIS method with

. The individual coordinates of the learning datasets were independently drawn from the normal distribution

which introduced a small fraction of longer distances between the target locations improving generalization overall considered densities of the benchmark instances222It would take roughly 21 days to generate the full dataset for on a single core, and thus 64 cores were used to speed up the process.. The distribution was also selected empirically, where the process was bootstrapped by observing the distributions of locations in sub tours of the random DTSP instances. We found WiSM-EA to be reasonably robust w.r.t. to the training distributions; nevertheless, future research should focus on training data generation methods.

The model training was done as follows. We split the dataset to training part () and the validation part (

). The loss function to train the surrogate model was the Mean Squared Error (MSE), which is a common choice for the regression. No regularization method such as L1, L2, or drop out 

[32] was used as overfitting was not a problem due to the size of the training data. We employed Adam [19] using recommended parameter values (the learning rate , and ). The training was terminated using early stopping

when there was no improvement in the validation loss for ten successive epochs. The single WiSM employed for all

pairs in the following experiments needed 134 epochs to converge, achieving MSE of on the validation set333The target values were standardized as usual when training neural networks..

Vi-B WiSM-EA Performance Evaluation

The evaluation results of the proposed WiSM-EA and its comparison with the other DTSP solvers is depicted in Figures 3, 6 and 9. Each data point represents the mean value of the real required computational time and the mean value of the corresponding normalized cost computed over 100 runs of the particular algorithm, i.e., ten runs per each of the ten random instances for each problem setup

. The boxes in the presented plots delimit lower and upper quartiles. Note that, in general,

increases with the decreasing density of as the problem becomes more similar to its underlying Euclidean TSP.

The WiSM and SA methods were run using multiple settings in order to study the tradeoff between the solution quality and computational requirements. In particular, for the WiSM-EA, we used eight values of the time limit S (evaluated by the stopCondition(S) method) ranging from 1 second to 600 seconds, denoting the algorithm configuration as WiSM-EAS. The total required computational time , as in its final stage WiSM-EA calls IRIS to compute the headings for the sequence found by the evolutionary search (Line 1, Algorithm 1). The SA was run for four sample sizes which are referred to as SA.

The benefit of the surrogate model is supported by a comparison with two more methods denoted the WIRIS-EA and IRIS-EA. The WIRIS-EA is the WISM-EA (including the caching) with a sole exception of using IRIS instead of its neural network surrogate model to compute the partial costs . The IRIS-EA simply combines the evolutionary algorithm presented in Section IV with IRIS approximation of the closed DTP of the complete sequence , i.e., .

Due to computational requirements of IRIS, its sampling precision was reduced to for both WIRIS-EA and IRIS-EA. However, in many cases, we were unable to get results competitive to other methods for lower values of the S. Nevertheless, the value of was selected empirically based on the overall best tradeoff between the DTP length approximation precision and the number of evaluations achievable in S. Notice that the final IRIS call (Line 1, Algorithm 1) was still executed using to show impact of the found sequences by the EA.

Both the methods significantly underperform the WiSM-EA with the windowing WIRIS-EA being a better alternative. The reason is evident from aggregated results over all densities for , where the WiSM-EA achieves a significantly higher rate of evaluations of the DTP solution cost per second in comparison to much slower WIRIS-EA and IRIS-EA with and evaluations per second, respectively.

Examples of the found solutions are shown in Fig. 16 and a detailed report of the evaluation follows.

Fig. 3: Normalized cost and real computational requirements for DTSP instances with the density and target location. Several values for the WiSM-EA indicate different time limit and for SA the particular number of the samples. Note the logarithmic scale of the time axis.
(a) high density instances, ,
(b) low density instances, ,
Fig. 6: Performance of the evaluated DTSP solvers in high () and low () density instances with target locations. The achieved normalized cost is shown according to the required computational time or a given time limit (for the WiSM-EA). For , the iWiSM-EA refers to the proposed WiSM-EA with the initialization of the population by the solution of the Euclidean TSP provided by the Concorde solver.
(a) small instances ,
(b) large instances ,
Fig. 9: Performance of the evaluated DTSP solvers in small () and large () instances with the density .
(a) WiSM-EA, ,
(b) WiSM-EA, ,
(c) WiSM-EA, ,
(d) Baseline, ,
(e) Baseline, ,
(f) Baseline, ,
Fig. 16: Solutions of the DTSP with target locations and density found by the proposed WiSM-EA (upper) and baseline (bottom). The target locations are depicted as small blue disks, and the found Dubins tour is in the red. Each column represents the same single instance, and thus the solutions can be directly compared. The final cost is computed using the IRIS with the maximum resolution .

The achieved results for the DTSP instances with the density and locations are shown in Fig. 3, where the Pareto front is formed by three methods only: the AA, baseline, and the proposed WiSM-EA. Taking only the shortest limit into consideration, all three methods give results in units of seconds. In particular, the AA in , the baseline in , and the WiSM-EA in . The WiSM-EA achieves a significantly better than the AA with and the baseline with . Besides, the normalized cost further improves to in approximately ten minutes (WiSM-EA600). The only rival in the solution cost is the SA method with the best value achieved by the SA8 in . The solution cost provided by the SA should theoretically decrease with the increasing number of samples; however, the observed behavior is different, and the opposite trend can be noticed for the SA16 and SA32. This behavior can be explained by too large instances of the underlying Generalized TSP and heuristic solution of the transformed TSP instance. Similar behavior was observed for all tested instances having .

Evaluation results for the DTSP instances with locations and high and low densities are depicted in Fig. 6. For the high density instances with , the results are similar to (compare Fig. 3 and Fig. (a)a). The AA reaches the mean , while the proposed WiSM-EA achieves significantly better results ranging from (WiSM-EA1) to (WiSM-EA600). LIO does not perform well, while the SA provides competitive results for the time limit around .

For the low-density instances with , the DTSP becomes closer to the corresponding Euclidean TSP because locations are relatively far from each other and the minimum turning radius constraint does not significantly affect the length of the final Dubins tour. It is also indicated by the results for the baseline with the sequence determined by the ETSP solver. The proposed WiSM-EA outperforms all other solvers, but the baseline is outperformed only for the time limits . Based on that, we employed a solution of the corresponding Euclidean TSP provided by the Concorde solver [5] in the initialization of the whole initial generation in the initializePopulation method (Line 1, Algorithm 1) of the proposed EA instead of a random permutation to improve the performance. The modified algorithm is denoted iWiSM-EA in Fig. (b)b. Although the iWiSM-EA provides outstanding performance in low-density instances, it performs similarly to the WiSM-EA for and . For large instances, we observed slightly higher final costs, which can be explained by the zero diversity of the initial population making the optimization prone to local minima (data not shown for readability). We conclude that the initialization of the iWiSM-EA approach is beneficial for low-density instances only.

Computational performance for instances with small () and high () number of locations is shown in Fig. 9. For small instances, the Pareto front is formed by the AA, SA, and WiSM-EA solvers. Note that WIRIS-EA achieves the best result for while its performance deteriorates for longer time limits which can be explained by the reduced precision using . Solutions found by the WiSM-EA1 have the mean cost that is significantly lower than for the AA with , the SA4 with , but also the SA8 with with the mean computational time . The AA is the fastest approach with the mean and the SA4 needs . For large instances with , the proposed WiSM-EA dominates the other methods. The WiSM-EA10 provides the mean in and reaches in approximately ten minutes.

Vi-C Discussion and Possible Future Work

Based on the reported results, the proposed approximator WiSM is a vital approach that enables the solution of the DTSP instances by a relatively simple evolutionary algorithm. The developed WiSM-EA scales significantly better than the other evaluated methods in both the problem size and density . For large instances, it dominates the other methods in the computational requirements and the quality of the found solutions. From a practical point of view, the WiSM-EA gives the best results in units of seconds for small and medium-sized instances, while for the large instances with hundreds of locations, it needs only low tens of seconds. To improve the rate of convergence for low-density instances of the DTSP, we suggest using the ETSP initialization of the WiSM-EA.

Regarding the other methods, the AA is the fastest algorithm, but it always provides and thus worse solutions than the baseline. The SA provides competitive results to the proposed WiSM-EA in medium and high-density instances, but it does not scale with the problem size and the number of the utilized samples, mostly because of the limitations of the underlying solver to the transformed Generalized TSP.

Although it is not the aim of this paper, our initial experiments indicated (data not shown), that WiSM is robust w.r.t. different neural network architectures. Nevertheless, future research should focus on finding possibly more effective regressors. The similar applies to approaches generating representative training data for WiSMs learning.

Regarding the future deployments of the proposed surrogate approximator of the Dubins tour costs, we believe it can also be utilized in other Dubins routing problems such as the Dubins Orienteering Problem [30] where many sequences have to be evaluated. Moreover, the model can be generalized for touring problems with neighborhood areas instead of single locations, where the recently introduced Generalized Dubins Interval Problem [34] can be utilized for the model learning using high-quality solutions, and then for solving the Dubins TSP with Neighborhoods [8]. Finally, the herein presented results motivate for future work on the approximation of costs of curvature-constrained tours for more complex vehicle models than Dubins vehicle, e.g., Bézier curves [10].

Vii Conclusion

We present a novel approach to solve the Dubins Traveling Salesman Problem (DTSP) by a relatively simple Evolutionary Algorithm that is based on the proposed surrogate approximator of the Dubins tour cost called WiSM. Even though collecting enough training data might take considerable time, once the dataset is built and the WiSM is trained, it can provide Dubins tour cost estimates in a very fast rate, which can be exploited by robust global search methods such as Evolutionary Algorithms. The developed WiSM-EA has been evaluated on DTSP instances of varying size and also the density of the locations to be visited. Based on the reported results, the WiSM-EA outperforms the existing state-of-the-art approaches in the quality of the found solutions and also computational requirements. The results demonstrate the proposed method scales with the problem size and density, and thus, it is a suitable heuristic for finding high-quality solutions of curvature-constrained routing problems with Dubins vehicle.


  • [1] D. Applegate, R. Bixby, V. Chvátal, and W. Cook (2007) The traveling salesman problem: a computational study. Princeton University Press. Cited by: §I.
  • [2] J. Bezanzon, S. Karpinski, V. Shah, and A. Edelman (2012) Julia: A fast dynamic language for technical computing. CoRR abs/1209.5145. Cited by: §VI.
  • [3] C. M. Bishop (2006) Pattern recognition and machine learning (information science and statistics). Springer-Verlag. External Links: ISBN 978-0-387-31073-2 Cited by: §V-B.
  • [4] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio (2014) Learning phrase representations using RNN encoder–decoder for statistical machine translation. In

    Conference on Empirical Methods in Natural Language Processing (EMNLP)

    pp. 1724–1734. External Links: Document Cited by: §V.
  • [5] D. Applegate, R. Bixby, V. Chvátal, and W. Cook (2003) Concorde TSP Solver. Note: [cited 22 Aug 2018] External Links: Link Cited by: §I, §VI-B, §VI.
  • [6] L. Davis (1985) Applying adaptive algorithms to epistatic domains. In

    International Joint Conference on Artificial Intelligence (IJCAI)

    pp. 162–164. Cited by: §IV.
  • [7] L. E. Dubins (1957) On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents. American Journal of Mathematics, pp. 497–516. Cited by: §I, §I, §II, §III, Problem 1.
  • [8] J. Faigl, P. Váňa, R. Pěnička, and M. Saska (2019) Unsupervised learning-based flexible framework for surveillance planning with aerial vehicles. Journal of Field Robotics 36 (1), pp. 270–301. External Links: Document Cited by: §I, §VI-C.
  • [9] J. Faigl, P. Váňa, M. Saska, T. Báča, and V. Spurný (2017) On solution of the dubins touring problem. In European Conference on Mobile Robots (ECMR), pp. 1–6. Cited by: §I, §I, §I, §I, §II, §III, §III, §III, §IV, §IV, §V-A, §V-A, §VI, §VI, §VI.
  • [10] J. Faigl and P. Váňa (2018) Surveillance planning with Bézier curves. IEEE Robotics and Automation Letters 3 (2), pp. 750–757. External Links: Document Cited by: §VI-C.
  • [11] X. Glorot, A. Bordes, and Y. Bengio (2011) Deep sparse rectifier neural networks. In International Conference on Artificial Intelligence and Statistics, pp. 315–323. Cited by: §V-B.
  • [12] X. Goaoc, H. Kim, and S. Lazard (2013) Bounded-curvature shortest paths through a sequence of points using convex optimization. SIAM Journal on Computing 42 (2), pp. 662–684. Cited by: §III.
  • [13] D. Guimaraes Macharet, A. Alves Neto, V. Fiuza da Camara Neto, and M. Montenegro Campos (2012) An evolutionary approach for the Dubins’ traveling salesman problem with neighborhoods. In

    Annual Conference on Genetic and Evolutionary Computation

    pp. 377–384. Cited by: §I.
  • [14] S. Hochreiter and J. Schmidhuber (1997) Long short-term memory. Neural computation 9 (8), pp. 1735–1780. Cited by: §V.
  • [15] J. H. Holland (1992) Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control and artificial intelligence. MIT Press. External Links: ISBN 978-0-262-08213-6 Cited by: §IV.
  • [16] S. Horng, S. Lin, L. H. Lee, and C. Chen (2013) Memetic algorithm for real-time combinatorial stochastic simulation optimization problems with performance analysis. IEEE Transactions on Cybernetics 43 (5), pp. 1495–1509. External Links: Document Cited by: §V-B.
  • [17] Y. Jin (2005-01) A comprehensive survey of fitness approximation in evolutionary computation. Soft Computing 9 (1), pp. 3–12. External Links: Document Cited by: §V-B.
  • [18] Y. Jin (2011) Surrogate-assisted evolutionary computation: recent advances and future challenges. Swarm and Evolutionary Computation 1 (2), pp. 61–70. Cited by: §V-B.
  • [19] D. P. Kingma and J. Ba (2015) Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), Cited by: §VI-A.
  • [20] P. Larranaga, C. M. H. Kuijpers, R. H. Murga, I. Inza, and S. Dizdarevic (1999) Genetic algorithms for the travelling salesman problem: a review of representations and operators. Artificial Intelligence Review 13 (2), pp. 129–170. Cited by: §IV, §IV.
  • [21] J. Le Ny, E. Feron, and E. Frazzoli (2012) On the Dubins Traveling Salesman Problem. IEEE Transactions on Automatic Control 57 (1), pp. 265–270. Cited by: §I.
  • [22] K. Helsgaun (2018) LKH solver 2.0.9. Note: [cited 23 Dec 2018] External Links: Link Cited by: §I, §VI.
  • [23] D. G. Macharet and M. F. M. Campos (2018) A survey on routing problems and robotic systems. Robotica 36 (12), pp. 1781–1803. External Links: Document Cited by: §I, §II.
  • [24] S. G. Manyam and S. Rathinam (2018) On tightly bounding the dubins traveling salesman’s optimum. Journal of Dynamic Systems, Measurement, and Control 140 (7), pp. 071013. Cited by: §I.
  • [25] S. G. Manyam, S. Rathinam, D. Casbeer, and E. Garcia (2017) Tightly Bounding the Shortest Dubins Paths Through a Sequence of Points. Journal of Intelligent & Robotic Systems, pp. 1–17. Cited by: §I, §I, §I, §III.
  • [26] M. Niendorf and A. R. Girard (2018-02) Exact and Approximate Stability of Solutions to Traveling Salesman Problems. IEEE Transactions on Cybernetics 48 (2), pp. 583–595. External Links: Document Cited by: §I.
  • [27] P. Oberlin, S. Rathinam, and S. Darbha (2010-12) Today’s traveling salesman problem. Robotics & Automation Magazine, IEEE 17 (4), pp. 70–77. Cited by: §I, §I.
  • [28] K. J. Obermeyer, P. Oberlin, and S. Darbha (2012) Sampling-based path planning for a visual reconnaissance unmanned air vehicle. Journal of Guidance, Control, and Dynamics 35 (2), pp. 619–631. External Links: Document Cited by: §I.
  • [29] K. J. Obermeyer (2009) Path planning for a UAV performing reconnaissance of static ground targets in terrain. In AIAA Guidance, Navigation, and Control Conference, pp. 10–13. Cited by: §VI.
  • [30] R. Pěnička, J. Faigl, P. Váňa, and M. Saska (2017) Dubins orienteering problem. IEEE Robotics and Automation Letters 2 (2), pp. 1210–1217. Cited by: §VI-C.
  • [31] K. Savla, E. Frazzoli, and F. Bullo (2005) On the point-to-point and traveling salesperson problems for Dubins’ vehicle. In American Control Conference, pp. 786–791. Cited by: §I, §I, §VI.
  • [32] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov (2014) Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research 15 (1), pp. 1929–1958. Cited by: §VI-A.
  • [33] P. Váňa and J. Faigl (2015) On the dubins traveling salesman problem with neighborhoods. In IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 4029–4034. Cited by: §I, §VI.
  • [34] P. Váňa and J. Faigl (2018) Optimal Solution of the Generalized Dubins Interval Problem. In Robotics: Science and Systems (RSS), Cited by: §VI-C.
  • [35] Y. Wang, D. Yin, S. Yang, and G. Sun (2018) Global and Local Surrogate-Assisted Differential Evolution for Expensive Constrained Optimization Problems With Inequality Constraints. IEEE Transactions on Cybernetics, pp. 1–15. External Links: Document Cited by: §I.
  • [36] X. Zhang, J. Chen, B. Xin, and Z. Peng (2014) A memetic algorithm for path planning of curvature-constrained UAVs performing surveillance of multiple ground targets. Chinese Journal of Aeronautics 27 (3), pp. 622–633. Cited by: §I, §IV, §VI.