Efficient Constraint Handling in Electromagnetism-Like Algorithm for Traveling Salesman Problem with Time Windows

The traveling salesman problem with time windows (TSPTW) is a variant of the traveling salesman problem in which each customer should be visited within a given time window. In this paper, we propose an electromagnetism-like algorithm (EMA) that uses a new constraint handling technique to minimize the travel cost in TSPTW problems. The EMA utilizes the attraction-repulsion mechanism between charged particles in a multidimensional space for global optimization. This paper investigates the problem-specific constraint handling capability of the EMA framework using a new variable bounding strategy, in which real-coded particle's boundary constraints associated with the corresponding time windows of customers, is introduced and combined with the penalty approach to eliminate infeasibilities regarding time window violations. The performance of the proposed algorithm and the effectiveness of the constraint handling technique have been studied extensively, comparing it to that of state-of-the-art metaheuristics using several sets of benchmark problems reported in the literature. The results of the numerical experiments show that the EMA generates feasible and near-optimal results within shorter computational times compared to the test algorithms.


Introduction
The traveling salesman problem with time windows (TSPTW) is an important variant of the well-known traveling salesman problem (TSP) in which each customer has a service time (i.e., the time that should be spent during the visit to the customer), which should start within a given time window. In the TSPTW, the time window of each visit is bounded by an earliest arrival time and a latest arrival time. The TSPTW can be defined as the problem of finding a minimum cost tour that starts and ends at the same depot, where each node should be visited exactly once within its time window. A nonnegative cost is associated with each arc, and the total cost can be taken as the route completion time, total travel time, or total traveled distance [1]. In this study, the total traveled distance will be considered the cost value. The TSPTW can be considered a special case of the capacitated vehicle routing problem with time windows (VRPTW) to which the relaxation of capacity constraints is applied.
The TSPTW has various practical real-world applications such as package delivery, school bus routing, scheduling, automated guided machines, and routing problems in the context of lean manufacturing. Savelsbergh [2] has shown that the TSPTW is NP-hard, and even finding a feasible route is a NP-complete problem. Nevertheless, early study focused on exact optimization techniques to solve the TSPTW. Baker [1] proposed an exact algorithm using a branch and bound approach in which lower bounds are determined by dual relaxations of the model. Dumas et al. [3] developed a dynamic programming approach, and Langevin et al. [4] described a two-commodity flow model including two complementary flows. More recently, a branch and cut algorithm [5], linear time dynamic programming [6], and constraint programming [7,8] were proposed to solve the TSPTW.
Because it is difficult to solve the TSPTW within acceptable computation times using exact methods, Heuristic and metaHeuristic techniques have been analyzed extensively in the literature. Carlton and Wesley Barnes [9] proposed a static-penalty-based tabu search to solve the TSPTW. Gendreau et al. [10] presented an insertion Heuristic, and Calvo [11] used a construction Heuristic based on both greedy insertion and local search. Furthermore, Da Silva and 2 The Scientific World Journal Urrutia [12] proposed a general variable neighborhood search (VNS) Heuristic that consists of a constructive stage for finding a feasible solution and an optimization stage during which feasible solutions are improved using a general VNS Heuristic. Ohlmann and Thomas [13] proposed a compressed annealing approach, which is a variant of simulated annealing and utilizes a variable penalty method, and López-Ibáñez and Blum [14] proposed a Beam-ACO algorithm, which hybridizes Beam search and ant colony optimization to solve the TSPTW via extensive computational analysis.
In this study, an electromagnetism-like algorithm (EMA) using a new variable bounding strategy is proposed to solve the TSPTW efficiently. The EMA approach basically emulates electromagnetism theory in physics, in which charged particles exert attractive or repulsive forces on each other [15]. The basic idea behind the algorithm is to force particles to search for the optimum in a multidimensional space by applying a collective force on them. In recent work, an EMA was used to train neural networks [16] to solve fuzzy relations [17]. Debels et al. [18] solved resource constrained scheduling problems by combining an EMA with scatter search. Tsou and Kao [19] used an EMA to control and optimize multiobjective inventory models. Maenhout and Vanhoucke [20] used an EMA to solve nurse scheduling problems. Jhang and Lee [21] used an EMA for array pattern optimization in the field of electrical engineering. Chang et al. [22] combined a GA with an EMA to solve single machine earliness/tardiness scheduling problems and showed that the hybrid EMA performs better than the plain EMA. Moreover, the researchers used a random key procedure to decode particles into feasible schedules. Wu et al. [23] applied an EMA using a random key procedure to solve the traveling salesman problem. Yurtkuran and Emel [24] solved capacitated vehicle routing problems by using an EMA that hybridizes a local search procedure. Naji-Azimi et al. [25] combined preprocessing procedure, mutation, and local search with an EMA to solve the well-known unicost set covering problem. In [26] EMA framework was introduced for solving the response time variability problem. Guan et al. [27] used an EMA to solve flow path design problems for automated guided vehicles. Jamili et al. [28] proposed a simulated annealing and electromagnetism-like mechanism hybrid framework to solve the periodic job shop scheduling problem. In [29] EMA was used to detect circles on figures. Su and Lin [30] introduced an EMA mechanism for feature selection. Lee and Chang [31] used an improved EMA to optimize fractional-order PID controllers. Furthermore, EMAs have also been applied to address multimode project scheduling under uncertainty [32] and nonlinear system control [33].
In this study, new constraint handling techniques are introduced to cope with time-window constraints in the TSPTW. To the best of our knowledge, our proposed constraint handling technique is the first in the literature in which the feasibility is maintained without using any extra feasibility operators. The main goal of using VBS is to narrow the unbounded search space to a bounded search space to reach feasible solutions effortlessly. Moreover, to analyze the effectiveness of the proposed algorithm, first, the constraint handling performance of the proposed EMA framework is analyzed, comparing it to that of the traditional EMA. Then, the modified framework is compared to state-of-theart algorithms using benchmark problems.
The rest of this paper is organized as follows. Section 2 provides a brief introduction and a mathematical formulation of the TSPTW. Section 3 discusses the traditional EMA, and the proposed EMA for the TSPTW is presented in Section 4. The computational results are summarized in Section 5, and, finally, Section 6 concludes the paper.

Traveling Salesman Problem with Time Windows
The traveling salesman problem with time windows can be briefly defined as follows. Let = ( , ) be an undirected complete graph, where = {0, 1, . . . } is the node (customer) set and = {( ), = 1, . . . , , = 1, . . . , , ̸ = } is the arc set. Node 0 denotes the depot, and denotes the number of customers. There is a cost associated with every arc ∈ . The cost generally represents the distance or time between nodes and , plus a service time at customer . In addition, each node ∈ has a time window [ , ], where denotes the earliest time and the latest time in which the service can begin. In most of the formulations, waiting times are permitted; that is, arrival to node before is feasible, but a waiting time till is applied. On the other hand, arrival to node after is not permitted. The TSPTW can be mathematically formulated as follows: decision variables: : position of node , where ∈ within the tour; : arrival time at node , ∈ ; : waiting time at node , ∈ ; parameters: : set of nodes; : travel time (distance) from node to , where , ∈ ; : service time at node , ∈ ; : earliest arrival time at node , ∈ ; : latest arrival time at node , ∈ ; subject to ∑ ∈ = 1, ∀ ∈ , ̸ = . (3) − + ≤ − 1, ∀ , ∈ , ̸ = 1, ̸ = . (5) The Scientific World Journal 3 ≤ + ≤ , ∀ ∈ .
Formula (2) denotes the objective function of the problem. The objective is to minimize the total time required to travel from each node to node and the service time at each node . Constraints (3) and (4) ensure that every node is visited only once. Constraint (5) aims to maintain the sequence of the route. Constraints (6) and (7) specify the service time windows, where denotes a large real number. Constraint (8) bounds the sequence for each node . Depending on the tight bounds in time-window constraints (7), TSPTW results in a tighter search space of feasible solutions. If an efficient neighborhood search strategy which is part of a metaHeuristic solver takes advantage of this feasible solution space, it is expected that the performance of such a metaHeuristic algorithm can be improved. Therefore, in the following sections, EMA will be used to host our proposed strategy for handling constraint (7).

Electromagnetism-Like Algorithm
The EMA is a population-based Heuristic method introduced by Birbil and Fang [15] for solving bound constraint optimization problems. The algorithm was inspired by the theory of electromagnetism in physics, in which there are attractive and repulsive forces between charged particles. The EMA can be used easily and effectively to solve optimization problems with bounded variables of the following form: where [L,U] = {x ∈ R | ≤ ≤ , = 1, . . . , }. In the problem formulation, x is a vector that represents a solution point position in -dimensional space in a population of position vectors.
denotes the variable of the th dimension (i.e., axis). Each variable has an upper and lower bound, and , respectively, and (x) indicates the objective function value (OFV) of the candidate solution x.
In the EMA, a candidate solution is associated with a charged particle in a multidimensional space using a realcoded position vector x. Index in each particle 's position vector x identifies a dimensional element , where = 1, . . . , and = 1, . . . , . In other words, and are the population size and the total number of variables, respectively. The OFV of the th candidate solution is calculated by using its position vector. The charge of particle , , depends on the quality of the OFV. The better the OFV of the particle is, the greater amount of charge the particle has. Each particle exerts a repulsive or attractive force on other population members according to the charges they carry. The resultant force F (1) Initialize (PopSize) (2) Set Iter ← 1 (3) While Iter < MaxIter do (4) LocalSearch(LsIter) (5) CalculateCharges() (6) CalculateForces() (7) Movement() (8) Set Iter ← Iter + 1 (9) End While is determined by calculating the vector sum of the forces exerted on a particle . Then, x is updated by F at each iteration. The key idea of the EMA is that, for a minimization problem, a candidate particle will attract particle if particle has a better OFV than particle (i.e., (x ) < (x )), whereas if (x ) < (x ), particle will repel particle .
The traditional EMA has four phases: (1) initialization, (2) calculation of particle charges and force vectors, (3) movement according to the resultant force, and (4) local search to exploit the local minima [15]. The general scheme of the algorithm is presented in Algorithm 1 (for more details, readers are referred to [15]). In Algorithm 1, represents the population size, and and are the maximum iteration number for the algorithm and the local search procedure, respectively. An () procedure is used to generate number of points randomly from the search space. A ℎ() procedure is applied to the particles to improve the solution quality and to force the algorithm to search for unvisited regions. Then, ℎ (), (), and V () procedures are applied to the particles within the population at each iteration.
In this study, the charge and force calculations and the subsequent particle movement procedures are implemented using the modified EMA proposed by [18]. Here, the charge and force calculations are not absolute-value-based; instead, relative charge and force calculations are used for each particle pair in the population. The details of the proposed algorithm will be described in the next section.

EMA for TSPTW
This section provides the details of the proposed EMA.

Charge and Force Calculations and Movement of Particles.
In our proposed EMA, the charge of particle is defined as , which is relative to that of particle , whereas it has been defined as in previous works [20,21,25,26,29,32,33]. The value of can be obtained by calculating the relative difference between the OFVs of particles and [18]: where (x worst ) and (x best ) are the worst and best OFVs in the population, respectively. For a minimization problem,

4
The Scientific World Journal if (x ) < (x ), then will be negative, and the reverse is true if (x ) > (x ). If (x ) = (x ), then will be zero.
After calculating , the force vector F exerted on particle by particle is calculated as follows: For a minimization problem, if particle is a better solution than particle , that is, (x ) < (x ), then particle will repel particle because will be negative.
In the modified EMA, first, and F are calculated for every combination of particle pairs; then, the resultant cumulative force on particle is determined by F = ∑ =1, ̸ = F . Once the cumulative force F exerted on particle is determined, the particle will move in the direction of F to yield improved solutions. A uniformly distributed random step 0 < < 1 is used to force the algorithm to explore unvisited regions. Similarly to the cooling effect in the simulated annealing algorithm, the current iteration number, , is used to decrease the step length as the algorithm proceeds. Additionally, a preset parameter , 0 ≤ ≤ 1, is used to control the cooling effect. The motion of particle along the F direction is calculated as follows: It is important to note that the proposed EMA uses an elitist strategy. In other words, the best solution's position vector in the population, x best , is preserved.

Solution Representation.
As mentioned above, the EMA was originally designed to cope with continuous optimization problems. In order to solve combinatorial optimization problems such as vehicle routing problems with EMA, real-coded position vectors (candidate solutions) have to be decoded into permutations of customers. To the best of our knowledge, most researchers have introduced a random key (RK) procedure into the EMA to facilitate solving combinatorial optimization problems [18,22,25,26]. The RK representation was proposed by [34] to maintain feasible solutions after crossover operations in genetic algorithms. In [34] a random number encoding structure was proposed for the chromosomal representation of solutions. The main advantage of using the RK procedure is that each candidate solution can be represented by real-coded values such that several metaHeuristic operators can be implemented without concern over feasibility issues. Because all position vectors are real-coded, integrating the random key procedure into the EMA is a very straightforward and easy process, thus making the EMA an efficient search algorithm for combinatorial optimization problems. A sample random key procedure is shown in Figure 1. In the random key procedure, when the real-coded coordinate values of the position vector are sorted in a nondecreasing order, the new permutation of the indexes of this position vector represents a route for the TSPTW as a sorted index. In Figure 1, because the smallest coordinate value of the position vector is 0.04 for index = 2, customer 2 will be visited first. The other customers are visited following the sorted index in a similar manner, and the resulting route will be 2 → 1 → 5 → 4 → 3.  (Figure 2(a)). Because a waiting time up to is applied for early visits, the earliest time ( earliest ) in which a customer can be left is the early time of that customer ( earliest ≥ ). Therefore, any customer sequence that ensures the following constraint will always be infeasible for a nonoverlapping customer pair:

Handling Time
where represents the set of customer pairs that customer precedes . In other words, if < , then customer should be visited before customer ; otherwise, if < , then customer should be visited before customer and any tour that contains a sequence in which is visited before customer is infeasible. VBS relies on the ability of the EMA to operate with bounded variables. In VBS, the time windows of customers are normalized between predetermined lower and upper global bounds [ , ], and the variables are then bounded within their corresponding normalized time windows [ nor , nor ], where nor = ( / min ), nor = ( / max ), min = min { }, and max = max { } and represents the customer number. Combining the VBS and the nondecreasing sorting step in RK, any solution point that has infeasible customer pairs with nonoverlapping time windows is thus eliminated from the candidate solution population. In other words, EMA is forced to search in feasible regions using VBS.

Penalty Strategy with Overlapping Time Windows.
However, this variable bounding strategy has drawbacks in The Scientific World Journal the case of highly overlapping time windows. The effect of VBS will decrease in going from nonoverlapping to overlapping windows and will be ineffective if the time windows of a customer pair ( , ) are fully overlapping ( = and = ) (see Figure 2(b) for an overlapping time window example).
To overcome the infeasibility problems for highly overlapping time windows, a penalty strategy is introduced. A penalty cost that is calculated from a linear penalty function is added to the OFV if the solution violates the time window of any customer. This penalty is assumed to be a linear function of the amount of time that is violated. The penalty cost is calculated as follows: where denotes the penalty cost at customer , represents the vehicle arrival time to customer , and is the penalty coefficient, which will be determined experimentally. By using the penalty approach, infeasible solutions will have high OFVs and will exert repulsive forces on better solutions. The effect of VBS and RK in eliminating time window infeasibilities is illustrated by considering a sample problem. Consider a TSP with 4 customers having time windows, as indicated in Table 1. For ease of analysis, we ignore the normalization step and time windows are directly used as the bounds. As shown in Figure 3, customers (1,2), (1,4), (2,4), and (3,4) have nonoverlapping time windows. By definition, customer 4 should be visited last because 4 ≥ , = 1, 2, 3. Because we bound index 4 of the position vector with customer 4's normalized time window [ nor , nor ], index 4 will always be larger than the other variables during the search process. As a result, when the real-coded coordinate values of the position vector are sorted in a nondecreasing order in RK, it will be ensured that customer 4 will be the last. Similarly, by combining VBS and RK, customer 1 will always precede customer 2 because 2 ≥ 1 . To summarize, by combining VBS and RK as a solution representation strategy, infeasibilities associated with nonoverlapping time windows are eliminated, and those infeasible regions are abandoned. Figure 4 shows the variable boundaries and possible positions  of customers after using the VBS and RK. Furthermore, a penalty approach will help to discern feasible solutions from the possible feasible solution set because infeasible solutions will have higher OFVs.  At the end of each iteration, boundary control is applied to the position vector of the particles to determine whether any boundary violations occur. The proposed algorithm adopts an absorbing bound-handling scheme, where, in the case of boundary violation, the corresponding variable is relocated to the bound. Algorithm 3 is the boundary control mechanism; as indicated, if the particle flies outside the boundary, the corresponding variable is set to its normalized early or late arrival time.

General Scheme of EMA.
The general scheme of the proposed EMA for solving the TSPTW is summarized in Algorithm 4. Each step is executed according to the explanations provided above.

Computational Results
The proposed algorithm, described in the previous sections, was implemented in Visual Basic. Net on a PC with an Intel Core 2 Duo CPU running at 2.0 GHz with 2 GB RAM for computational experiments. Two types of experiments were carried out to assess the effectiveness of the proposed EMA. First, to evaluate the performance of the VBS, the EMA with VBS and that without VBS are compared with respect to selected benchmark instances. Second, the proposed EMA is compared to state-of-the-art metaHeuristics using an extensive set of benchmark instances reported in the literature.

The Effectiveness of Variable Bounding Strategy for the TSPTW.
To understand the role of VBS in finding feasible solutions, the proposed EMA with VBS and a penalty strategy (EMA-VP) was compared with the EMA with only a penalty strategy (EMA-P). Because the level of time window overlapping is the key criterion in analyzing the effect of VBS, six different problem instances selected from the benchmark set provided by Potvin and Bengio [35] are categorized as exhibiting a low, average, or high level of time window overlapping. An explicit indicator value of the overlap level (VOL) is calculated by adding two percentages calculated from problem instances, that is, (i) the percentage of time windows of two or more customers that intersect over the total time line (min and max ) and (ii) the ratio of customers with an overlapping time window of at least one unit length. In other words, the length of the overlapped time and number of customers having overlapping time windows are calculated as basic indicators. Therefore, a VOL of 200 corresponds to the full overlap of time windows, whereas 0 denotes a nonoverlapping problem.

Comparison Using Benchmark Instances.
To verify the effectiveness of the proposed algorithm, we performed an extensive analysis using several benchmark instances. The results of this study are presented in a manner similar to that of state-of-the-art studies by López-Ibáñez and Blum [14] and Ohlmann and Thomas [13] for an accurate and objective comparison. The three benchmark sets used to test the proposed EMA-VP mechanism are as follows.
(1) Benchmark Set 1. The benchmark set was introduced by Potvin and Bengio [35], which was originally developed for VRPTW problems by Solomon [36]. This set includes 29 problems and is the most widely used benchmark set for the TSPTW. Originally, the set included 30 problems; however, there was a conflict among the researchers about the node number and best solution value for problem RC 204.1 [13,14].
(Thus, it is not considered here.) (See Table 4.) (2) Benchmark Set 2. The set of 70 problems in seven instance classes was introduced by Langevin et al. [4]. These instances range from 20 to 60 customers ( Table 5).  (3) Benchmark Set 3. The set of benchmark instances was introduced by Dumas et al. [3] This set consists of 135 instances, and the customer numbers range from 20 to 200 (Table 6).
Because parameter calibration is the key task in metaHeuristic applications, we performed a set of pilot studies to determine a good set of parameters for the EMA-VP mechanism. After these preliminary computational studies, the parameters were set as follows: = 500, = 1.5, = 0.35, and = 50. The results are presented in terms of relative percentage deviation (RPD), which is calculated as 100×(EMA value−the Best known value)/(the Best known value). Because the EMA (particularly the movement procedure) is stochastic in nature, each result is reported as the average of 10 runs. As presented in the studies by [13,14], both the mean RPD and standard deviation RPD of relative percentage deviation results and the mean CPU times (second) are reported here as well.
The results of the comparison between the EMA-VP and novel metaHeuristics for the benchmark set 1 are summarized  in Table 4. The first column of Table 4 represents the problem name; denotes the number of customers in the problem; VOL indicates the value of the time window overlap level and BKV is the best known solution value of the problems presented in the literature. Furthermore, whereas bold-typed values of BKV are the optimal values reported by others, the values indicated by an asterisk are the optimal values determined using CPLEX 12.1 in this study. The EMA-VP is compared to the Beam-ACO [14], compressed annealing (CA) proposed by [13], dynamic programming (DP) [6], and the best values reported in the studies by Gendreau et al. [10] and Calvo [11] as Heuristic. As shown in Table 4, the EMA-VP finds the optimal or the best known values for 19 out of 30 instances without any solution value variability. The Beam-ACO outperforms the EMA-VP in some of the high VOL instances; nevertheless, the differences between the mean RPDs (i.e., RC 204.2, RC 204.3, and RC 208.3) are quite small. Moreover, the EMA-VP and CA yield very similar results in all of the instances, and the results that are obtained by the EMA-VP are better than those obtained by DP and Heuristic. Furthermore, the EMA-VP is able to find a feasible solution for all instances.
The Scientific World Journal 9  The results of the benchmark set 2 are presented in Table 5. These results are the averages of 10 instances of 10 runs, as in the other studies performed by [13,14]. The EMA-VP is compared to compressed annealing (CA) [13], Beam-ACO [13], and the best known values (BKV) [36].
The EMA-VP yields promising results and the optimal values for the first 4 instances (i.e., = 20 and = 40) and achieved the best known values for the large instances (i.e., = 60). Moreover, no variance in the solution quality is reported. As a result, the EMA-VP is compatible with CA and Beam-ACO  for the benchmark set 2. It is important to note that all of the instances in the problem set are low and average VOLs; thus, the EMA-VP shows promising convergence rates. Table 6 reports the results of the EMA-VP for the benchmark set 3. The results are presented as averages of 10 instances in each class and 10 runs for each instance. Table 6 compares the EMA-VP with the exact method developed by Dumas et al. [3], Beam-ACO [14], and CA [13], and the last column is titled Heuristic, which represents the best value obtained by the algorithms developed by [9][10][11] For a better evaluation of metaHeuristic algorithms, not only the solution quality but also the computation times should be investigated. However, it is not very easy to make an objective comparison between metaHeuristics because both the programming languages and the machine configurations are generally not comparable, and, in most studies, the complexities of the algorithms are not reported. Nevertheless, an approximate comparison can be made based on the MFlop (million floating point operations per second) values of the processors on which the algorithms were coded and run [37]. The CA algorithm was coded in C++ and implemented on an Intel Pentium 4 CPU operating at 2.66 GHz [13], Beam-ACO was implemented in C++ on an Intel Xeon X3350 processor with a 2.66 GHz CPU [14], and Heuristic [11] was coded in FORTRAN and executed on an Intel 486 CPU operating at 66 MHz. The MFlop values of the processor speeds based on the benchmark values obtained from the site http://www.netlib.org/benchmark/linpackjava/timings list .html and the reported and normalized mean CPU times on the benchmark sets for the algorithms are summarized in Table 7. DP [6] is not included in the comparison because the CPU was not reported in the study. As shown in Table 7, the EMA-VP is faster than Beam-ACO and Heuristic for all of the benchmarks. However, the performances of Heuristic and the EMA-VP are very similar on the benchmark set 1. Moreover, CA is faster than the EMA-VP on benchmark set 3, and EMA-VP outperforms CA on the other benchmark sets. Table 7 clearly shows that, in general, the proposed EMA-VP is an effective algorithm and outperforms other novel algorithms in 9 out of 11 test cases in terms of computational time.

Conclusion
This paper has presented an EMA with a variable bounding strategy (VBS) as a novel constraint handling technique for solving the traveling salesman problem with time windows. The EMA is an easy-to-code, straightforward metaHeuristic algorithm that emulates the attraction-repulsion interactions of charged particles in analogy to Coulomb's law in electromagnetic theory. The proposed algorithm uses two important approaches to handle time-window constraints, the penalty approach and VBS. VBS is one of the main contributions of this study, in which the upper and lower bounds of a variable are set using the corresponding time window for serving a customer. The main goal of using VBS is to narrow the unbounded search space to a bounded search space to reach feasible solutions effortlessly. We clearly show that our approach competes other approaches reported in the literature. An extensive computational analysis using well-known benchmark instances shows that the EMA-VP converges to feasible regions in a search space and finds the best known or near-optimal results. Furthermore, the EMA-VP outperforms other novel metaHeuristics in terms of computational time. Future work may involve combining the VBS technique with other metaHeuristics using real-coded particles as in particle swarm optimization, differential evolution, or artificial bee colony algorithms to solve combinatorial optimization problems that have constraints similar to time windows, such as scheduling with precedence constraints or resource constraint project management.