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 highquality 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 curvatureconstrained path connecting the locations in the sequence is minimal [31]. The curvatureconstrained 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 NPhard [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 curvatureconstrained path connecting the locations is one of six possible maneuvers that can be found by a closedform solution [7]. The generalization to an arbitrary number of locations is known as the Dubins Touring Problem (DTP). The stateoftheart DTP solvers are based on a sampling of the possible heading angles [25]. In [9], an informed sampling method, which we further denote the IterativelyRefined 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 decoupledapproaches, 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 samplingbased 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 searchspace 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 searchspace 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 highquality solutions of the IRISbased 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 fixedsize 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 WiSMEA) 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 highquality 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 WiSMEA 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
(1) 
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)
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)
(3) 
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 highquality 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]. Samplingbased 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 curvatureconstrained 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.
Highquality 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 highquality 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.
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 1–1, Algorithm 1). The mutation operator is the Simple Inversion Mutation (SIM) [15, 20] which reverses a random subsequence. The crossover method implements the wellknown Order crossover (OX) [6] which copies a random subsequence 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 1–1, 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, samplingbased 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 variablelength as sequences (e.g., LSTM
[14] or GRU [4]), a collection of training data would be nontrivial 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 variablesize input, we propose a decomposition of the cost function to fixedsize subproblems based on the sliding window technique.Va Sliding Windowing based Cost Estimation of the DTP
The proposed idea is to evaluate the cost of Dubins tour for small subsequences 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
(4) 
We propose to utilize a sliding window that starts at each target location as it enables to address arbitrarilylength instances for which the number of locations is not necessarily divisible by the utilized window size . The whole Dubins tour is divided into subsequences 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
(5) 
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 subsequence is considered as the Open DTP [9], which can be defined as the continuous optimization problem
(6) 
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:
(7) 
Comparing the particular DTP cost of the th window defined by a subsequence and the corresponding cost , it can be shown that
(8) 
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.
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
(9) 
Based on (8), the estimation is a lower bound of the optimal cost and we get
(10) 
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 fixedsized subtours, enables approximation of the relatively slow samplingbased 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 domainspecific knowledge of the addressed DTSP using results of the DTP.
VB Windowing Surrogate Model (WiSM) Approximation of the Cost Estimation
The key idea of the proposed WiSMbased 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 realtime response of the proposed DTSP solver can benefit from pretrained surrogate model [16], and therefore, the model is learned offline using a large training dataset based on highquality solutions of the openloop .
The WiSM can be implemented using virtually any type of regression model, but we consider the MultiLayer Perceptron (MLP)
[3]in this paper. It is because the feedforward neural network supports fast evaluation of all partial costs
in a single batch. Thus, it can take advantage of highly optimized matrix multiplication routines^{1}^{1}1We 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 hasinputs, several hidden layers of nonlinear units such as ReLU
[11], 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 WiSMEA. The solved DTSP instances are of various sizes and different densities of the locations.
The performance of the WiSMEA is compared with the existing DTSP approaches. Namely, the Alternating Algorithm (AA) [31], the Local Iterative Optimization (LIO) [33], and the Samplingbased 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
(11) 
where is the cost provided by the baseline.
The headings for the baseline and also for the WiSMEA (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 WiSMEA was implemented in Julia language [2], run with O3 and checkbounds=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.
Via Learning Setup of the Proposed WiSM
The choice of an appropriate window size of the WiSM deals is a tradeoff 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 VB 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 instances^{2}^{2}2It 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 WiSMEA 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 stoppingwhen 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 set^{3}^{3}3The target values were standardized as usual when training neural networks..ViB WiSMEA Performance Evaluation
The evaluation results of the proposed WiSMEA 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 WiSMEA, 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 WiSMEA_{TS}. The total required computational time , as in its final stage WiSMEA 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 WIRISEA and IRISEA. The WIRISEA is the WISMEA (including the caching) with a sole exception of using IRIS instead of its neural network surrogate model to compute the partial costs . The IRISEA 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 WIRISEA and IRISEA. 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 WiSMEA with the windowing WIRISEA being a better alternative. The reason is evident from aggregated results over all densities for , where the WiSMEA achieves a significantly higher rate of evaluations of the DTP solution cost per second in comparison to much slower WIRISEA and IRISEA with and evaluations per second, respectively.
Examples of the found solutions are shown in Fig. 16 and a detailed report of the evaluation follows.
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 WiSMEA. 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 WiSMEA in . The WiSMEA achieves a significantly better than the AA with and the baseline with . Besides, the normalized cost further improves to in approximately ten minutes (WiSMEA_{600}). The only rival in the solution cost is the SA method with the best value achieved by the SA_{8} 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 SA_{16} and SA_{32}. 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 WiSMEA achieves significantly better results ranging from (WiSMEA_{1}) to (WiSMEA_{600}). LIO does not perform well, while the SA provides competitive results for the time limit around .
For the lowdensity 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 WiSMEA 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 iWiSMEA in Fig. (b)b. Although the iWiSMEA provides outstanding performance in lowdensity instances, it performs similarly to the WiSMEA 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 iWiSMEA approach is beneficial for lowdensity 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 WiSMEA solvers. Note that WIRISEA 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 WiSMEA_{1} have the mean cost that is significantly lower than for the AA with , the SA_{4} with , but also the SA_{8} with with the mean computational time . The AA is the fastest approach with the mean and the SA_{4} needs . For large instances with , the proposed WiSMEA dominates the other methods. The WiSMEA_{10} provides the mean in and reaches in approximately ten minutes.
ViC 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 WiSMEA 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 WiSMEA gives the best results in units of seconds for small and mediumsized instances, while for the large instances with hundreds of locations, it needs only low tens of seconds. To improve the rate of convergence for lowdensity instances of the DTSP, we suggest using the ETSP initialization of the WiSMEA.
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 WiSMEA in medium and highdensity 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 highquality 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 curvatureconstrained 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 WiSMEA 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 WiSMEA outperforms the existing stateoftheart 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 highquality solutions of curvatureconstrained routing problems with Dubins vehicle.
References
 [1] (2007) The traveling salesman problem: a computational study. Princeton University Press. Cited by: §I.
 [2] (2012) Julia: A fast dynamic language for technical computing. CoRR abs/1209.5145. Cited by: §VI.
 [3] (2006) Pattern recognition and machine learning (information science and statistics). SpringerVerlag. External Links: ISBN 9780387310732 Cited by: §VB.

[4]
(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] (2003) Concorde TSP Solver. Note: [cited 22 Aug 2018] External Links: Link Cited by: §I, §VIB, §VI.

[6]
(1985)
Applying adaptive algorithms to epistatic domains.
In
International Joint Conference on Artificial Intelligence (IJCAI)
, pp. 162–164. Cited by: §IV.  [7] (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] (2019) Unsupervised learningbased flexible framework for surveillance planning with aerial vehicles. Journal of Field Robotics 36 (1), pp. 270–301. External Links: Document Cited by: §I, §VIC.
 [9] (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, §VA, §VA, §VI, §VI, §VI.
 [10] (2018) Surveillance planning with Bézier curves. IEEE Robotics and Automation Letters 3 (2), pp. 750–757. External Links: Document Cited by: §VIC.
 [11] (2011) Deep sparse rectifier neural networks. In International Conference on Artificial Intelligence and Statistics, pp. 315–323. Cited by: §VB.
 [12] (2013) Boundedcurvature shortest paths through a sequence of points using convex optimization. SIAM Journal on Computing 42 (2), pp. 662–684. Cited by: §III.

[13]
(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] (1997) Long shortterm memory. Neural computation 9 (8), pp. 1735–1780. Cited by: §V.
 [15] (1992) Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control and artificial intelligence. MIT Press. External Links: ISBN 9780262082136 Cited by: §IV.
 [16] (2013) Memetic algorithm for realtime combinatorial stochastic simulation optimization problems with performance analysis. IEEE Transactions on Cybernetics 43 (5), pp. 1495–1509. External Links: Document Cited by: §VB.
 [17] (200501) A comprehensive survey of fitness approximation in evolutionary computation. Soft Computing 9 (1), pp. 3–12. External Links: Document Cited by: §VB.
 [18] (2011) Surrogateassisted evolutionary computation: recent advances and future challenges. Swarm and Evolutionary Computation 1 (2), pp. 61–70. Cited by: §VB.
 [19] (2015) Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), Cited by: §VIA.
 [20] (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] (2012) On the Dubins Traveling Salesman Problem. IEEE Transactions on Automatic Control 57 (1), pp. 265–270. Cited by: §I.
 [22] (2018) LKH solver 2.0.9. Note: [cited 23 Dec 2018] External Links: Link Cited by: §I, §VI.
 [23] (2018) A survey on routing problems and robotic systems. Robotica 36 (12), pp. 1781–1803. External Links: Document Cited by: §I, §II.
 [24] (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] (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] (201802) 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] (201012) Today’s traveling salesman problem. Robotics & Automation Magazine, IEEE 17 (4), pp. 70–77. Cited by: §I, §I.
 [28] (2012) Samplingbased 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] (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] (2017) Dubins orienteering problem. IEEE Robotics and Automation Letters 2 (2), pp. 1210–1217. Cited by: §VIC.
 [31] (2005) On the pointtopoint and traveling salesperson problems for Dubins’ vehicle. In American Control Conference, pp. 786–791. Cited by: §I, §I, §VI.
 [32] (2014) Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research 15 (1), pp. 1929–1958. Cited by: §VIA.
 [33] (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] (2018) Optimal Solution of the Generalized Dubins Interval Problem. In Robotics: Science and Systems (RSS), Cited by: §VIC.
 [35] (2018) Global and Local SurrogateAssisted Differential Evolution for Expensive Constrained Optimization Problems With Inequality Constraints. IEEE Transactions on Cybernetics, pp. 1–15. External Links: Document Cited by: §I.
 [36] (2014) A memetic algorithm for path planning of curvatureconstrained UAVs performing surveillance of multiple ground targets. Chinese Journal of Aeronautics 27 (3), pp. 622–633. Cited by: §I, §IV, §VI.