1 Introduction
In several applications such as taxiservices (Caramia et al., 2002; Fabri and Recht, 2006), emergencydispatches (Yang et al., 2005) and warehousepicking (SmolicRocak et al., 2009) the next location is unknown. In such cases the multiple Travelling Salesman Problem (mTSP) becomes dynamic. That means during the “building” of the route new nodes are revealed. Hence, there is either an implicit or explicit time dimension involved.
In the multiple Travelling Salesman Problem (mTSP) salesmen (vehicles) visit customers (nodes). Each customer will be visited once by one and only one salesman. The salesmen start and end their routes at a single depot (source node). The objective is to minimise the total distance travelled by the salesmen. The mTSP is called balanced when each vehicle visits “approximately” the same number of customers. This can be achieved by introducing an upper and lower limit of customers a salesman has to visit. Gouveia and SalazarGonzález (2010) discuss these capacity bounds in the context of the Vehicle Routing Problem. MartinezSykora and Bektaş (2015) introduce a general approach that transforms nodebalanced routing problems into generalised balanced TSPs. This allows the reduction of arcs in the underlying graph. Bektaş et al. (2019) solve balanced vehicle routing using a polyhedral analysis and a branchandcut algorithm. An interesting integer program formulation of the balancedmTSP is given in Kara and Bektas (2006). I adapted this program in Garn (2020) for the balancedmTSP and compared it to two heuristics, one of them is the precursor for this work. Furthermore, that work offers a continuous approximation formula for the balancedmTSP.
The aim of this work it to provide a continuous approximation model (CAM) for the balanceddynamic mTSP (BDmTSP). Hence, the rest of the paper is organised as follows. Section 2classifies types of dynamics. Sequential dynamics will be used to develop two BDmTSP heuristics (BDCVH and BDAVH in Section 3). They are compared against each other and an exact balancedstatic (BS) mTSP method (Section 3.3). Here, a number of classic TSP/mTSP testinstances are used. It must be stressed that the comparison is intended to derive the difference between dynamic and static mTSPs rather than whether one is better than the other one. This will give first insights of how much balancedstatic and balanceddynamic mTSPs differ from each other. Continuous approximation models for the static mTSP are reviewed in Section 4
. Most approaches focus on uniformly distributed nodes in the Euclidean plane, which will be adapted for the dynamic case. To the best of the author’s knowledge this is the first CAM model for the balanceddynamic case. Moreover, this appears to be that for the first time a structured Machine Learning approach was used to derive a CAM relationship. The here introduced BDmTSP continuous approximation models describe the total distance depending on the number of vehicles, customers and sequentialdynamics.
2 Dynamics
2.1 Evolution and Quality
Pillac et al. (2013) reviewed dynamic vehicle routing problems, which is also relevant to the mTSP. They suggested to differentiate between information evolution and quality. I will propose definitions to make their ideas more tangible. Information evolution refers to changes in the data available to the planner such as new customer requests. Let be the sequence of customers changing over an evolving time horizon . denotes the finite maximum number of customers for all . Assume that is a discrete time event where an information evolution occurred. Without restricting generality let be ordered over time, i.e. if then information evolution happened before or at the same time as .
Definition 2.1.
Node information evolution is the process of associating time to nodes . is a monotonic increasing sequence. is ordered such that each is mapped to .
Hence, the information evolution is the “visibility” of customers at a certain point in time. The above definition automatically ensures that related information such as the location or distances between customers becomes time dependent. That means, the customer locations can be deterministic but revealed over time. In the Euclidean plane with coordinates and . The above definition means that customers and locations are ordered according to a time dimension. A special case, which will be used later, is that all time steps are discrete and equidistant.
The sequence of vehicles (fleet) and related information can evolve over time as well. Hence, vehicle information evolution is defined similarly. Later sections explaining heuristics and CAMs will assume that the fleets existence is permanent.
Information quality is related to uncertainty in the input data.
Definition 2.2.
Node information quality
is the probability distribution associated to the positions of nodes or distances between nodes.
An example for requiring an exact (deterministic) location is the picking of an item. The eventual location of an object in a dropoff zone is an example for an uncertain (stochastic) position.
Table 1 shows the two dimensions and its related classes.
quality  
input  deterministic  stochastic  
evolution 
static  SD  SS 
dynamic  DD  DS 
The “classic” mTSP falls into the SD class, i.e. static (no) information evolution and deterministic information quality. The exact balanced mTSP introduced in Garn (2020) is in this class. Even mTSPs with timewindows are in this class. In the context of routing  stochastic means that information such as customer locations follows a known random distribution. Other typical stochastic VRP factors are demand, times and pickups. On a strategic level all the static continuous approximation models (see Section 4) fall into the SS class. Any VRP or mTSP implementation that is supplied with data based on a random distribution is called stochastic. The balanced dynamic closest vehicle heuristic (BDCVH) and balanced dynamic assignment vehicle heuristic (BDAVH) are able to operate in the dynamicdeterministic (DD) class where information evolution is dynamic and information quality is deterministic (all locations are known). The DD class means for the mTSP that some of the customers will not be known in advance, but at some time during the execution of the route. The dynamicstochastic (DS) class implies for the mTSP that new customers are revealed during the building, and their locations or distances are randomly distributed. Note that some of the customer locations could be known in advance. The knowledge of the underlying random distribution allows algorithms to anticipate (predict) the occurrence of new locations. The continuous approximation given later is a prime example for the DSclass. The BDCVH and DCAH operate in this class as well.
Pillac et al. (2013) mention dynamic VRP solution methods and divides them into continuous and periodic approaches. Here, continuous, implies as soon as new information becomes available reoptimisation is performed. An example for continuous reoptimisation is Gendreau et al. (1999). They use a parallel tabu search to find solutions for the dynamic VRP. Their premises are that new information is an event that occurs, and a reoptimisation is triggered, which allowed them to adapt a static mTSP implementation. The heuristics focus on accommodating minor changes. For comparison purposes they proposed a few algorithms. Their “insertion” algorithm adds a new customer to the planned routes such that the additional cost is minimised. A “rebuild” algorithm and a few more adaptive tabu search methods are briefly mentioned as options. Hence, if most information is known in advance and only minor adaptation are expected these algorithms are a good option. This type of dynamics will be called random node insertion.
Genetic algorithms (GA) are a popular approach, and so they find their presence in the continuous DmTSP world. Cheung et al. (2008) used a GA to update any changes to existing routes. Again, this method is based on the assumption of having most of the information available. Similarly, Haghani and Jung (2005)
developed a GA to solve a DVRP with timedependent travel times. DmTSP based on periodic reoptimisation fall back on static mTSP solution methods. Whenever, new information is available all data is assumed to be static and deterministic and a classic mTSP solution method is executed. Probably one of the first (if not the first) periodic reoptimisations is due to
Psaraftis (1980). He applied dynamic programming and periodic reoptimisations to solve a dial and ride problem.The dialaride problem (DARP) (Lois and Ziliaskopoulos, 2017; Kirchler and Calvo, 2013) is closely related to the mTSP. It differs from the DmTSP by having a pickup and dropoff location. However, those locations could be aggregated into a single abstract node; which allows to transform the dialaride problem into an asynchronous mTSP. Typical DARP formulations emphasis time windows. Surprisingly, almost all DARP formulations are in the static category. Cordeau and Laporte (2003) reviewed the DARP and confirmed this view. Madsen et al. (1995) is one of the few works which analyse the DARP in a dynamic environment intended for online scheduling. The practical requirement was to handle up to 300 requests having 24 vehicles available to be scheduled. Furthermore, a solution had to be returned within 2 seconds. Multiple objectives criteria and constraints had to be considered. They started their work based on Jaw et al. (1986) algorithm. Again, their approach falls into the class of random insertion heuristics, where a known schedule is improved. The heuristic looks for all feasible insertions in the existing routes and adds the new request such that a minimal change to the objective occurs. If a feasible insertion cannot be found the request (job) is not served.
Given the gap of a systematic dynamics’s classification with the exception of information evolution and quality several dynamic scopes will be proposed.
2.2 Scopes
Measures for the dynamic scopes will be introduced, which include absolute, relative and dependent values.
Absolute dynamics is defined as the fixed number of customer requests known at time with the exception of the final period. To accomodate the number of vehicles the term absolute dynamics is proposed: .
Example 2.1 (Absolute dynamics).
Customers ordered by “reveal” time and let the absolute dynamics be . That means at five customers and related deterministic or stochastic quality information such as their locations are known. So, at “any” point of time the customers are known and the previously revealed customers . Obviously, when reaching the end () of the scenario only the remaining customer information is revealed.
Assume there are 3 vehicles and the absolute dynamics is customers. This is equivalent to absolute dynamics .
An interesting special case occurs when and the next customer locations are randomly distributed in the “vicinity” of the last visited customers, because this is equivalent to random walks in higher dimensions. It would be interesting to investigate the relationship between dimensional Gaussian random walks and the DmTSP with the closest vehicle heuristic. This could lead to interesting insight and relationships to the BlackScholes derived approaches.
If the total number of customer requests is known then relative node dynamics can be defined:
(1) 
where controls the fraction of revealed data. The use of these rational numbers has the advantage that one can convert between relative and absolute dynamics without rounding. Alternatively, a given percentage value is converted to the relative node dynamics via . Generally, it is more convenient to allow percentages as relative node dynamics.
Example 2.2 (Relative dynamics).
Customers ordered by “reveal” time and let the relative dynamics be . That means 5% of the customer data is revealed at each time step .
To consider the number of vehicles needs to be incorporated. Equally to before, relative dynamics is defined:
(2) 
Example 2.3 (relative dynamics).
The number of vehicles is , customers are ordered by “reveal” time and let the relative dynamics be . That means 6% of the customer data is revealed at each time step, giving 2% location choice to each vehicle.
In most cases it is necessary to convert to absolute node dynamics: , and .
The next level of dynamics occurs, when the number of customers varies at each time step . I will call this variable dynamics . is a sequence of customers visible at . Note, that at at . Generally, due to the finite nature of the formulation the following set of inequalities hold: .
Example 2.4 (variable dynamics).
The number of vehicles is , customers are ordered by “reveal” time and let the variable dynamics be . That means at customers are visible. At information about customers is available. At information about customers is available. Now the visibility reduces to customers at . This is interesting because it offers considerations such as allowing memory of . Alternatively, once visibility is lost it could the location of could change.
As mentioned before the most frequently used type of dynamics is random insertion. In essence customers are known and is added. Hence, this is a special case of the absolute dynamic scope with . It should be noted that in the previously reviewed literature time windows for customer tasks were considered. This add another timeline to the problem, or adds constraints to the existing timeline . To distinguish these two the term eventknowledgetimeline and scheduledtimeline are proposed. Future work shall look at the dynamics of the mTSP with timewindows.
On top of these fundamental dynamics the stochastic elements of customers’ quality of the location (e.g. position improves over time) must be taken into account. There are two aspects  where and when. A trivial case is that the customer location is static over the entire period. Each customer’s location can change over time (“moving target”). A potential approach is allowing on top of inserting customers the removal of customers. Additionally, the accuracy of the static or changing location of each customer needs to be considered. Now, these measures could apply to classes of customers or the entire customer set.
Vehicle dynamics can be considered using a varying number of active vehicles over time. Additionally, uncertainty in the position of the vehicles position or state can occur as well.
Figure 1 gives an overview of the abovementioned dynamics that can occur in the mTSP.
3 BDmTSP Heuristics
We will consider two greedy balanceddynamic mTSP (BDmTSP) heuristics. The first heuristic allocates the vehicle closest to “visible” customers. The second heuristic assigns vehicles minimising their distance.
The input for the algorithms is a distance matrix , the number of vehicle and the absolute dynamics . The balancing threshold (capacity limit) is determined automatically. Please note that the algorithms can be easily amended to contain the capacity as input and as variable. The output of the algorithms is a list of routes.
The dynamic algorithms are examined using and introducing several testinstances. The difference of the dynamic results to the static exact mTSP testinstances are explained. The heuristics will be used to derive a CAM.
3.1 Closest Vehicle
This algorithm is the dynamic version of the CVH. I introduced the CVH in Garn (2020). This algorithm uses absolute customer dynamics, i.e. is the number of customers visible at . As shown in the previous section relative dynamics and dynamics measures can be converted to absolute dynamics. It should be noted that the algorithm can be easily adapted to variable customer dynamics. vehicles are available. The customer nodes are balanced between those vehicles. We will assume that customers and their locations are revealed sequentially. This happens when all vehicles have completed their assigned tasks (customer visits). It is assumed that the total number of customers is known in advance. can be deduced from the vehicle’s capacity limit and the balancing constraint. Algorithm 1 shows the details for this implementation.
It begins with requiring some input such as . In realworld implementations would not be known in advance. For road networks shortest path algorithms can be used to derive . In most of the subsequent testinstances customer locations were given and the Euclidean distance between customer location and was derived . Sometimes it may be more practical to compute during execution as represented in line 7.
The algorithm ensures to return routes for all vehicles. Again in an online version it is recommended to use the route output from line 16.
The first five lines of the algorithm are initialisation and an overall loop start. Line 6 controls the nodes visibility (information evolution). This represents the limitations imposed by sequentialdynamics. Consequently, the distance matrix under consideration is reduced to (line 7). As mentioned before, for realworld implementations, this is the point where the new customer locations are known (rather than revealed) and is computed instead of being extracted from . Although preprocessed distance matrix constitutes computational savings. However, for the study of dynamics this implementation detail does not make a difference. Line 812 assign the closest vehicle to the visible nodes. Line 9 returns the row and column of the minimum matrix element. This is equivalent of identifying the closest vehicle. In line 10 that vehicle is assigned to the customer () using the assignment matrix . Note it can happen that . In order to prevent the choice of the same vehicle and customer the respective row and column are set to . Line 13 returns the original node numbers by making use of the matrix multiplication. Line 1422 add the vehicle to the route and enforce the balancing constraint. Line 23 removes the visited nodes from all not visited nodes. Syntactically, the case may need to be considered (depending on the programming language). Line 24 updates the fromnodes with the tonodes used in the last step, leaving the unused vehicles (some from nodes) as they are.
3.2 Dynamic Assignment Vehicle Heuristic (BDAVH)
This algorithm builds on the conceptual framework of the BDCVH. Line 812 were responsible for assigning the closest vehicle, which is a greedy heuristic. A better solution is obtained by assigning available vehicles such that the distance is minimised. This is achieved with the following binary program (if ):
min  (3)  
subject to  
In the case of the constraints change slightly and the following program is solved instead
min  (4)  
subject to  
Consequently, the algorithm for the BDAVH is the same apart from line 812 being replaced with the above binary programs, i.e. if then program (3) else program (4).
The subsequent section gives insights about the performance of the heuristics.
3.3 TestInstances
The balanced dynamic closest vehicle heuristic (BDCVH) and balanced dynamic assignment vehicle heuristic (BDAVH) are examined using several testinstances. Some of them are new (available on http://wiki.smartana.org/index.php/MTSP) and others were adapted from TSPLIB (http://elib.zib.de/pub/mptestdata/tsp/tsplib/tsp/). For the first set of mostly small testinstances (, except for and ) the number of vehicles was chosen arbitrarily. The second set contains medium testinstances and the number of vehicles is varied between 2 and 5. Furthermore, absolute sequential dynamics and absolute sequential dynamics were considered. The absolute sequential dynamics were . The absolute dynamics were obtained by setting the node visibility for all instances to . Table 2 shows the results for the first set of testinstances.
garn9m2  

BDAVH  68.9  66.1  64.8  44.8  44.8  44.8  44.8 
BDCVH  68.9  66.1  66.5  44.8  44.8  44.8  44.8 
garn13m3L4  
BDAVH  90.7  89.9  92.1  71.0  79.4  79.4  79.4 
BDCVH  91.4  91.0  93.2  71.9  87.8  87.8  87.8 
garn13m3L6  
BDAVH  90.7  89.9  92.1  71.0  79.4  79.4  79.4 
BDCVH  91.4  91.0  93.2  71.9  87.8  87.8  87.8 
garn20m3  
BDAVH  122.0  117.1  119.3  107.2  107.2  113.0  137.2 
BDCVH  123.3  120.0  120.9  106.0  106.0  106.7  140.3 
bays29m4  
BDAVH  4.93k  4.13k  3.98k  4.12k  3.41k  3.36k  3.91k 
BDCVH  5.31k  4.26k  4.14k  3.75k  3.53k  3.48k  3.93k 
berlin52m5  
BDAVH  25.7k  18.9k  17.1k  16.3k  15.2k  15.8k  17.4k 
BDCVH  26.8k  22.0k  18.2k  15.8k  15.9k  15.9k  17.6k 
eucln100m7  
BDAVH  42.8k  37.7k  31.6k  29.0k  29.5k  29.6k  29.6k 
BDCVH  44.7k  37.4k  30.7k  29.7k  31.0k  30.9k  31.4k 
lin318m20  
BDAVH  343.3k  322.5k  315.9k  253.3k  243.1k  207.8k  203.6k 
BDCVH  349.2k  332.0k  324.0k  253.9k  252.2k  212.7k  209.7k 
The instance garn9m2 shows that both BDCVH and BDAVH are identical. It is also interesting to compare them to their static versions, which I proposed in Garn (2020). I used the same testinstances and the results are shown in (Garn, 2020, p4, Table 1). Note that when dynamics the static CVH must give the same results as the BDCVH. However, we will focus on the exact static mTSP for the following comparisons. The exact static mTSP solution value is 44.8, which the dynamic version achieves with . The other garn instances in a dynamic setting differ substantially from the the static mTSP solution values. The static exact solution values for garn13m3L4, garn13m3L6, garn20m3 are 57.7, 48.9, 90.5. For garn13m3 dynamics an added distance of roughly 50% is observed. The total route distance differences between static and dynamic TSP are in the region of 30% to 50% for instances bays29m4 and berlin52m5 where static solution values are 2603 and 9845.2 respectively.
Interesting is that with some testinstances increased knowledge of future events does not always lead to better results. In almost all testinstance the absolute dynamics of
In the appendix Table 5 shows the corresponding absolute sequential dynamics. More medium testinstances are given in the appendix, with Table 6 and Table 7 showing results for absolute and absolute sequential dynamics respectively. In general, it can be observed that the BDAVH delivers shorter routes (see Table 3 for average absolute percentage differences).
mrelative  absolute  

small  2.74%  3.59% 
medium  1.43%  1.59% 
However, occasionally the BDCVH returns better solutions, e.g. instance bayes29m4 with . This is justifiable by observing that although the assignment step is optimal the overall algorithm is a heuristic. Hence, the assignment can send the vehicles to an unfavourable location for subsequent customer destinations.
In summary this section introduced two heuristics, which can be used for onlinerouting. The BDAVH returns slightly better solution values than the BDCVH. Dynamic solution values for a selection of testinstances are about 30% to 50% larger than their static exact counterparts.
4 Continuous Approximation Model
In this section we will give the functional relationship that relates the solution distance to the number of salesmen (vehicles), customers (nodes) and sequential dynamics.
4.1 Related Work
One of the first researches addressing the above for the static TSP was the work by Beardwood et al. (1959) called a shortest path through many points. This led to the famous BeardwoodHaltonHammersley theorem, which gives an asymptotic formula for the length of a TSP route.
Theorem 1 (The Beardwood–Halton–Hammersley (BHH) Theorem).
Let
be a set of random variables in
independently and identically distributed with bounded support. Then the length of a shortest TSP tour through the points satisfies(5) 
where is the absolutely continuous part of the distribution of the and depends on but not on this distribution.
Its usage influenced probability theory, physics, computer science and operational research. The case when
is of interest to us, i.e. a TSP length is ‘almost always’ asymptotically proportional to , where represent the number of points in a bounded plane region with area . The exact value of is still unknown (Applegate et al., 2006, p23), but was approximated as 0.7313 for . Applegate et al. (2006, p497ff)cite a few more estimates of
from various authors and provided a table with between 100 and 2,500 nodes and between 0.7765 and 0.7241 respectively. Most approaches to find a function for the TSP use heuristics or regression models. Building on the above work Çavdar and Sokol (2015) gave a distributionfree TSP tour length estimation model for random graphs. Their approach used a regression model based on sampling probability distributions.An interesting empirical formula for the CVRP was developed by Christofides (1971) , where the depot is in the centre and points are uniformly distributed in a square. Here, represents the general location to reach a point and the detour distance. , where is the area of interest, and is the number of items. This was the starting “stone” for Daganzo (1984) , who consider a special case of the CVRP, which agrees with the mTSP definitions given earlier. They derived a formula for the CVRP by extending the above TSP findings. This was achieved by clusterfirst and routesecond, i.e. a heuristic model. They emphasise the importance of choice of shape (e.g. slenderness) for the clusters. Their length formula was expressed as:
(6) 
where is the average Euclidean distance of the node locations to the depot, the number of points in an sector, the total area, and is the density of the area. In (Garn, 2020, p5,eq14) I proposed a continuous approximation model for the balanced static mTSP:
(7) 
where is a an discrete 100 by 100 area and customers are uniformly distributed on the grid.
Franceschetti et al. (2017) review literature on continuous approximation models in freight distribution management. They identified Beardwood et al. (1959) and Daganzo (1984) as seminal works (see above for details). Their concise review contains several more formulations and applications of continuous approximation models containing several valuable ideas and formulations. However, their work did not identify an approximation formula for the dynamicmTSP, which indicates a potential gap in the bodyofknowledge. Similarly, Ansari et al. (2018) discuss advancements in continuous approximation models for logistics and transportation systems by focusing their review on literature between 1996 and 2016. Their work covers a fast range of applications and problem instances. It appears that approximation models for the dynamic mTSP are sparse. Nevertheless, there is one report by Erera and Daganzo (2003)
, which introduces a “threshold global sharing” scheme utilising a realtime reoptimisation control for the VRP. At first it appears to have similarities with the CVH. However, it differs because of an initial partitioning approach, which is typical for Daganzo’s studies. Roughly speaking, their model is based on a spherical area. Customers’ density, expected demand and standard deviation are used. Vehicles’ capacity is set. Additionally, a buffer factor
is given which defines the number of standard deviations for the total customer demand. Overall, their approximation model has a high level of complexity, which makes it difficult to use for verification purposes.4.2 Continuous BDmTSP approximation
The idea is to run the BDAVH on uniformly distributed customers in the Euclidean plane. Regressing on the route lengths and feature matrices will provide the approximation formula. The analysis will be restricted to the unit square area. The advantage of this is the scalability to other area sizes.
Figure 2 and 3 show the average distances found by the BDCVH and BDAVH respectively. These heuristics used sequential dynamics.
(a)  (b)  (c) 
(a)  (b)  (c) 
It can be seen that the BDAVH finds better results in general. Taking averages over all () percentage deviations show that BDAVH is 1.61% better than BDCVH. It is interesting to note the upper bound , if sequential dynamic visibility is one and customers are always in opposite corners. The instance has the upper bound 708.5 and BDAVH distance 135.7. The lower bound is trivial. The previous two figures indicate a distinction between the cases and may be necessary. In general, the distance increases with the number of vehicles if , but in the case instance has lower distance values than . The above considerations and analyses identified the BDAVH as algorithm of choice. This algorithm is used to solve 30 testinstances for each configuration. A configuration is an element from:
Example 4.1 (Configuration 57).
is the 57 configuration having 3 vehicles, 100 customers and sequential dynamics of 15. Here, 100 customer locations are randomly created 30 times. Leading to 30 distance values, where is the average distance.
The total number of configurations is . The configurations are captured in . The total number of testinstances is . Since, the average travel distance for each configuration is taken .
The realworld model (systematic information) is explained by:
(8) 
where is the output (response),
is the input (known also as predictors, variables, features) and
is the random error. All inputs and outputs are real with and . Here, is the number of observations and is the number of features.The prediction model (estimate for ) generates predictions using existing (real) input (rows , columns ):
(9) 
is influenced by the degree of flexibility, e.g used features.
The quality of the fitted can be described using the rootmeansquarederror (RMSE):
(10) 
where and is the th row. The meanabsolutepercentageerror (MAPE) focuses on relative deviations and is commonly found in forecasting:
(11) 
where . Note: some authors define the MAPE as . The main issue of the MAPE is its sensitivity with small values, and that it is not defined when any . However, the advantage is its explainability.
We will derive an approximation formula (model) using multivariate regression, i.e.
(12) 
Here, is the total travelled distance and are features. represents “observed” values using the BDAVH algorithm. The features will vary depending on our approach. The prediction model is given by:
(13) 
A naive approach is to use the features asis. That means the number of vehicles , customers and sequential dynamics . This defines . I will refer to these as the basefeatures (or feature basis), i.e.
consists out of three basefeature vectors. Figure
4 (a) compares the observed and modelled distances with each other.(a)  (b) 
Distances modelled with lineare regression compared to observed ones from the BDAVH using (a) features asis, (b) feature map.
In a perfect model the points (observed, modelled) lay on the black line, i.e. the modelled distances are the same as the observed ones. However, it can be seen that this is not the case. Moreover, the RMSE is 12.64 and the MAPE is 26.41%, which is improvable.
We know from our previous and related work that the distance is proportional to the square root when dynamics are disregarded. Additionally, from Figure 2 there is polynomial behaviour visible. Hence, I build the features systematically using combinations of powers. The powers are ^{2}^{2}2Powers perform equally well: RMSE 1.67, MAPE 1.69%. . The resulting feature map consists out of 64 features:
This increases the flexibility of the model. Figure 4 (b) compares the modelled and observed distances. The RMSE is 1.55, which is low considering that is 51.9. The MAPE is 1.30% which means that the model is very well fitted. The computation of inverse of returns a matrix that is close to singular or badly scaled (absolute value magnitudes range between and ). Using the pseudoinverse instead leads to absolute magnitude range between and . However, the RMSE and MAPE increase to 4.55 and 7.41% respectively.
In order to improve the interpretability of the model and avoid overfitting  regularisation can be applied. Three common approaches are subset selection, shrinkage and dimension reduction (James et al., 2013). The bestsubset selection approach requires the fitting of regression models. This is computationally infeasible. There are two alternatives: forward or backward stepwise selection. These are computationally efficient and lead to similar results. I have used the backward stepwise selection algorithm. Choosing the optimal model can be achieved using Colin Lingwood Mallows’s crossvalidation prediction criteria , Akaike information criteria (AIC), Bayesian information criteria (BIC) or adjusted. Figure 5 shows , BIC and adjusted for the best model for each feature step.
Models with 0, 1 and 2 features were omitted in the display due to their large criteria values and large RMSEs ().
The crossvalidation prediction errors suggest a model with 62 features,
which appears to be a poor choice considering that the RMSE changes only slight
between feature 20 (rmse: 1.56, mape: 1.39%) and 62 (rmse: 1.55 , mape: 1.27%).
The case with three selected features is of interest,
because this is the same number of features as in the naive approach but with a RMSE of 5.44 instead of 12.64 (MAPE: 9.83% instead of 26.41%).
Recall the three basefeatures are vehicles , customers and dynamics .
The selected features (as polynomials of the basefeatures) for our first model are:
polynomial  b 

0.39051  
0.05477  
0.00023 
An advantage is the ease of using this model via:
(14) 
Example 4.2 (Predicted distance with CAM).
Let be an input configuration. The observed average distance is , which was computed using the BDAVH. The predicted average distance is obtained from equation (14). The corresponding percentage error is 4.8%.
being in all terms emphasises the importance of the total number of nodes (customers). However, unlike to the static mTSP the length is not proportional to the square root of customers. A reason could be that the visibility of the customers is limited by sequential dynamics. However, the square root is applied to the sequential dynamics scope basefeature. This makes sense, since this represents the visible nodes, which are used when assigning the vehicles to the customers. The number of vehicles is found in the term.
The recommended model according to BIC derived from the backward selection algorithm is the model with 16 features and a RMSE of 1.68 (MAPE: 2.11%) is shown in Table 4.
This (second) model contains the previously mentioned model including the three features. can be found in five terms confirming the importance of sequential dynamics. The total number of nodes (customers) is present in almost all terms with the exception of term 7.
A better compromise between model quality and interpretability
is model three with 9 selected features:
polynomial  b 

0.52829  
0.29958  
0.17818  
0.08168  
0.03651  
0.02354  
0.01102  
0.00927  
0.00126 
This model has a RMSE of of 2.35 and a MAPE of 2.82%.
The CAM models can be used to approximate the expected distance when nodes are uniformly distributed in the unit square area. The three proposed models vary in flexibility and accuracy. The first model is simple to use (three features) and has a MAPE of 9.83%. The second model has an accuracy of 2.11% (MAPE) but requires 16 features. The third model is a compromise using 9 features and having an accuracy of 2.82%.
5 Conclusion
Dynamic routing is essential in many real world scenarios. This work focused on dynamics for the balanced mTSP. However, this work can be easily applied to the capacitated vehicle routing problem (CVRP). It is closely related to dialaride problem (DARP) as explained in Section 2.2. Several types and scopes of dynamics were proposed. This work focused on the sequential dynamic scope. Future work can investigate variable dynamics. Two algorithms: balanceddynamic closest vehicle heuristic (BDCVH) and balanceddynamic assignment vehicle heuristic (BDAVH) were developed. It would be interesting to compare the BDAVH to insertionheuristics. Several testinstances (including derivates from the TSPLIB) gave insights about the behaviour of dynamic routing. The proposed testinstances and solutions can be used as benchmark reference. The distances indicated that in general the BDAVH algorithm leads to slightly better solutions (about 3%). A comparison to static routing on these testinstances gave an idea about the expected additionally travelled distance (up to about 50% to an exact static mTSP solution).
A continuous approximation model (CAM) for the BDAVH allows to predict the expected travel distance, when the number of vehicles, number of customers and sequential dynamics scope is given (or can be estimated). Three models were derived using a Machine Learning approach. Each offering a different degree in flexibility (ease of use) and accuracy (prediction quality). Currently, the method was used using uniformly distributed customers in the Euclidean plane. However, it is possible to use any other stochastic customer location distribution. Currently, the models were derived to operate for medium sized testinstances  in a space limited by the number of vehicles, number of customers and scope of dynamics. These limitation were chosen arbitrarily, i.e. setting different limits will return models that can be applied to larger testinstances.
This work’s algorithms were designed to work with sequential time consideration. In particular, discrete time events and a finite time horizon were suggested. However, it would be interesting to consider stochastic time distributions for customers and vehicles. Focus to the dynamics of vehicles needs to be addressed in future work. Once these fundamental dynamics are well understood, Finite State Machines (FSM) and Discrete Event Simulations (DES) may proof to be useful tools in developing general dynamicrouting solutions. Other future investigations could include timewindows for the mTSP and make used of eventknowledge and schedule timelines.
6 Acknowledgement
My thanks goes to the inspiring and creative HysterYale UK team, who’s work on the next generation of intelligent autonomous forklifttrucks motivated me to look for mTSP solutions. Mark’s vision of forklifttrucks bargaining for the next job inspired the closest vehicle mTSP heuristic, which is suitable for the onlinemTSP, i.e. a mTSP without knowing the pallet tasks (or customer locations) in advance. Chi (Ezeh) encouraged me to keep a wholistic view that reached from strategic mTSP formulations to vehicles’ Artificial Intelligence. Chris drew my attention to concurrency issues, that vehicles encounter when using the same path. There were many more inspirations, which demonstrated the importance of the onlinemTSP and its solution approaches in practical applications.
I am grateful for the anonymous reviewer of my work (Garn, 2020); who provided many interesting pointers to valuable research articles, which was the motivation for this work.
7 Appendix
instance  algorithm  

garn9m2  BDAVH  44.8  44.8  44.8  44.8  44.8  44.8 
BDCVH  44.8  44.8  44.8  44.8  44.8  44.8  
garn13m3L4  BDAVH  89.9  79.4  79.4  79.4  79.4  79.4 
BDCVH  91.0  87.8  87.8  87.8  87.8  87.8  
garn13m3L6  BDAVH  89.9  79.4  79.4  79.4  79.4  79.4 
BDCVH  91.0  87.8  87.8  87.8  87.8  87.8  
garn20m3  BDAVH  117.1  107.2  130.9  116.3  116.3  116.3 
BDCVH  120.0  106.0  134.0  120.2  120.2  120.2  
bays29m4  BDAVH  4.49k  4.12k  4.24k  3.73k  3.62k  3.46k 
BDCVH  4.79k  3.75k  4.05k  3.64k  3.99k  3.64k  
berlin52m5  BDAVH  25.7k  17.1k  15.2k  17.4k  16.2k  15.5k 
BDCVH  26.8k  18.2k  15.9k  17.6k  16.8k  16.9k  
eucln100m7  BDAVH  37.1k  39.3k  33.1k  30.6k  29.6k  29.9k 
BDCVH  39.5k  38.4k  30.5k  29.8k  30.9k  30.9k  
lin318m20  BDAVH  211.7k  261.4k  298.6k  343.3k  333.1k  322.5k 
BDCVH  215.9k  264.9k  300.8k  349.2k  342.2k  332.0k 
References
 Advancements in continuous approximation models for logistics and transportation systems: 1996–2016. Transportation Research Part B: Methodological 107, pp. 229–252. Cited by: §4.1.
 The traveling salesman problem: a computational study. Princeton university press. Cited by: §4.1.
 The shortest path through many points. In Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 55, pp. 299–327. Cited by: §4.1, §4.1.
 Balanced vehicle routing: polyhedral analysis and branchandcut algorithm. European Journal of Operational Research 273 (2), pp. 452–463. Cited by: §1.
 Routing a fleet of vehicles for dynamic combined pickup and deliveries services. In Operations Research Proceedings 2001, pp. 3–8. Cited by: §1.
 A distributionfree tsp tour length estimation model for random graphs. European Journal of Operational Research 243 (2), pp. 588–598. Cited by: §4.1.
 Dynamic routing model and solution methods for fleet management with mobile technologies. International Journal of Production Economics 113 (2), pp. 694–705. Cited by: §2.1.
 Distribution management: mathematical modelling and practical analysis. New York, Hafner. Cited by: §4.1.
 The dialaride problem (darp): variants, modeling issues and algorithms. Quarterly Journal of the Belgian, French and Italian Operations Research Societies 1 (2), pp. 89–101. Cited by: §2.1.
 The distance traveled to visit n points with a maximum of c stops per vehicle: an analytic model and an application. Transportation science 18 (4), pp. 331–350. Cited by: §4.1, §4.1.
 A dynamic scheme for stochastic vehicle routing. Report, Georgia Institute of Technology. Cited by: §4.1.
 On dynamic pickup and delivery vehicle routing with several time windows and waiting times. Transportation Research Part B: Methodological 40 (4), pp. 335–350. Cited by: §1.
 Continuous approximation models in freight distribution management. Top 25 (3), pp. 413–433. Cited by: §4.1.
 Closed form distance formula for the balanced multiple travelling salesmen. External Links: 2001.07749 Cited by: §1, §2.1, §3.1, §3.3, §4.1, §6.
 Parallel tabu search for realtime vehicle routing and dispatching. Transportation science 33 (4), pp. 381–390. Cited by: §2.1.
 On the vehicle routing problem with lower bound capacities. Electronic Notes in Discrete Mathematics 36, pp. 1001–1008. Cited by: §1.
 A dynamic vehicle routing problem with timedependent travel times. Computers & operations research 32 (11), pp. 2959–2986. Cited by: §2.1.
 An introduction to statistical learning. Vol. 112, Springer. Cited by: §4.2.
 A heuristic algorithm for the multivehicle advance request dialaride problem with time windows. Transportation Research Part B: Methodological 20 (3), pp. 243–257. Cited by: §2.1.

Integer linear programming formulations of multiple salesman problems and its variations
. European Journal of Operational Research 174 (3), pp. 1449–1458. Cited by: §1.  A granular tabu search algorithm for the dialaride problem. Transportation Research Part B: Methodological 56, pp. 120–135. Cited by: §2.1.
 Online algorithm for dynamic dial a ride problem and its metrics. Transportation research procedia 24, pp. 377–384. Cited by: §2.1.
 A heuristic algorithm for a dialaride problem with time windows, multiple capacities, and multiple objectives. Annals of operations Research 60 (1), pp. 193–208. Cited by: §2.1.
 Transformations of nodebalanced routing problems. Naval Research Logistics (NRL) 62 (5), pp. 370–387. Cited by: §1.
 A review of dynamic vehicle routing problems. European Journal of Operational Research 225 (1), pp. 1–11. Cited by: §2.1, §2.1.
 A dynamic programming solution to the single vehicle manytomany immediate request dialaride problem. Transportation Science 14 (2), pp. 130–154. Cited by: §2.1.
 Time windows based dynamic routing in multiagv systems. IEEE Transactions on Automation Science and Engineering 7 (1), pp. 151–155. Cited by: §1.
 Online dispatching and routing model for emergency vehicles with area coverage constraints. Transportation Research Record 1923 (1), pp. 1–8. Cited by: §1.
Comments
There are no comments yet.