Research Article  Open Access
Efficient Constraint Handling in ElectromagnetismLike Algorithm for Traveling Salesman Problem with Time Windows
Abstract
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 electromagnetismlike algorithm (EMA) that uses a new constraint handling technique to minimize the travel cost in TSPTW problems. The EMA utilizes the attractionrepulsion mechanism between charged particles in a multidimensional space for global optimization. This paper investigates the problemspecific constraint handling capability of the EMA framework using a new variable bounding strategy, in which realcoded 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 stateoftheart metaheuristics using several sets of benchmark problems reported in the literature. The results of the numerical experiments show that the EMA generates feasible and nearoptimal results within shorter computational times compared to the test algorithms.
1. Introduction
The traveling salesman problem with time windows (TSPTW) is an important variant of the wellknown 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 realworld 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 NPhard, and even finding a feasible route is a NPcomplete 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 twocommodity 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 staticpenaltybased 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 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ópezIbáñez and Blum [14] proposed a BeamACO algorithm, which hybridizes Beam search and ant colony optimization to solve the TSPTW via extensive computational analysis.
In this study, an electromagnetismlike 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. NajiAzimi et al. [25] combined preprocessing procedure, mutation, and local search with an EMA to solve the wellknown 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 electromagnetismlike 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 fractionalorder 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 timewindow 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 stateoftheart 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.
2. 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 is the node (customer) set and is the arc set. Node 0 denotes the depot, and denotes the number of customers. There is a cost associated with every . 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 , ;
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 timewindow 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).
3. ElectromagnetismLike Algorithm
The EMA is a populationbased 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 .
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 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 . Index in each particle ’s position vector identifies a dimensional element , where and . 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 is determined by calculating the vector sum of the forces exerted on a particle . Then, is updated by 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., ), whereas if, particle will repel particle.
The traditional EMA has four phases: initialization, calculation of particle charges and force vectors, movement according to the resultant force, and 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 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 absolutevaluebased; 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.
4. EMA for TSPTW
This section provides the details of the proposed EMA.
4.1. 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 and are the worst and best OFVs in the population, respectively. For a minimization problem, if , then will be negative, and the reverse is true if . If , then will be zero.
After calculating , the force vector exerted on particle by particle is calculated as follows:
For a minimization problem, if particle is a better solution than particle , that is, , then particle will repel particle because will be negative.
In the modified EMA, first, and are calculated for every combination of particle pairs; then, the resultant cumulative force on particle is determined by . Once the cumulative force exerted on particle is determined, the particle will move in the direction of to yield improved solutions. A uniformly distributed random step 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 , , is used to control the cooling effect. The motion of particle along the 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, , is preserved.
4.2. 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, realcoded 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 realcoded values such that several metaHeuristic operators can be implemented without concern over feasibility issues. Because all position vectors are realcoded, 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 realcoded 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.
4.3. Handling Time Window Constraints
In the proposed EMA, two approaches are combined to prevent infeasible routes: (i) new variable bounding strategy (VBS) and (ii) penalty approach. VBS is used to eliminate infeasible candidate solutions when the problem consists of nonoverlapping time windows for nodes, whereas a penalty strategy is used to cope with infeasibilities resulting from overlapping time windows.
4.3.1. VBS with Nonoverlapping Time Windows
Nonoverlapping time window infeasibilities can be described as follows. Consider two customers and with time windows and , respectively. These time windows are said to be nonoverlapping if and only if either or , where by definition (Figure 2(a)). Because a waiting time up to is applied for early visits, the earliest time () in which a customer can be left is the early time of that customer (). Therefore, any customer sequence that ensures the following constraint will always be infeasible for a nonoverlapping customer pair: 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.
(a)
(b)
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 , where , , , and 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.
4.3.2. Penalty Strategy with Overlapping Time Windows
However, this variable bounding strategy has drawbacks in 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 , . Because we bound index 4 of the position vector with customer 4’s normalized time window , index 4 will always be larger than the other variables during the search process. As a result, when the realcoded 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 . 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.

4.4. Initialization and Boundary Control
The EMA framework begins with the initialization mechanism. The Initialize() procedure generates PopSize() solutions as a starting population using the normalized time windows. The procedure is shown in Algorithm 2, where and denote the population size and length of the position vector (i.e., number of variables), respectively, and draws samples from a uniform distribution. Each position vector is initialized randomly within the time window of each corresponding variable, .

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 boundhandling 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.

4.5. 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.

5. 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 stateoftheart metaHeuristics using an extensive set of benchmark instances reported in the literature.
5.1. 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 (EMAVP) was compared with the EMA with only a penalty strategy (EMAP). 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 () 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.
The categorized VOLs are presented in Table 2. Furthermore, Table 3 summarizes the selected problem parameters and the corresponding VOLs. Problems with a similar number of customers from different classes (RC 201.3, RC 202.1, RC 203.2, and RC 204.2), a small problem with an average VOL (RC 205.1), and a relatively more complex problem (RC 208.1) are selected for the experiments. In Table 3, denotes the total number of customers and the depot. Figures 5, 6,7, 8, 9, and 10 show the ratio of the number of feasible solutions to the total population as the algorithm proceeds from each benchmark problem. For those experiments, the population size is set to 25 and the penalty coefficients and are set to 1.5 and 0.35, respectively. A single run represents 1000 iterations, and the ratio of feasible solutions is recorded at the end of 25 iterations. Box plots show 10 independent data, each of which represents an average of 25 consecutive runs.


The experiments show that the EMAVP clearly outperforms EMAP (Figures 5–10). The following conclusions can be drawn from the figures. In all six experiments, the EMAVP demonstrates better performance than the EMAP as expected. In all six experiments, EMAVP finds feasible solutions, whereas the EMAP finds a limited number of feasible solutions only for the smallest problem (RC 205.1). Except for the benchmarks with high VOL (RC 204.2 and RC 208.1), the VBS generates a feasible initial population with a minimum ratio of 0.1, whereas the EMAP always starts searching with an infeasible initial population. The EMAVP quickly converges to a significant ratio of feasible solutions in the first 150–250 iterations and finds a maximum of 100% (RC 205.1) and a minimum of 3% (RC 204.2) feasible solutions in the final populations. The negligible reduction (clearly observed for benchmark problem RC 203.2, Figure 7) in the ratio of feasible solutions in the first 50 iterations indicates that particles are generally stagnant at the bounds under high magnitudes of resultant forces. Particles are more likely to leave the search space because relatively high resultant force values are applied in the early phases of the searching process as force magnitudes are reduced iteratively (see (13), Section 4.1).
5.2. 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 stateoftheart studies by LópezIbáñez and Blum [14] and Ohlmann and Thomas [13] for an accurate and objective comparison. The three benchmark sets used to test the proposed EMAVP 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 EMAVP mechanism. After these preliminary computational studies, the parameters were set as follows: , and . The results are presented in terms of relative percentage deviation (RPD), which is calculated as . 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 and standard deviation of relative percentage deviation results and the mean CPU times (second) are reported here as well.
The results of the comparison between the EMAVP 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 boldtyped 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 EMAVP is compared to the BeamACO [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 EMAVP finds the optimal or the best known values for 19 out of 30 instances without any solution value variability. The BeamACO outperforms the EMAVP 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 EMAVP and CA yield very similar results in all of the instances, and the results that are obtained by the EMAVP are better than those obtained by DP and Heuristic. Furthermore, the EMAVP is able to find a feasible solution for all instances.
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 EMAVP is compared to compressed annealing (CA) [13], BeamACO [13], and the best known values (BKV) [36]. The EMAVP yields promising results and the optimal values for the first 4 instances (i.e., and ) and achieved the best known values for the large instances (i.e., ). Moreover, no variance in the solution quality is reported. As a result, the EMAVP is compatible with CA and BeamACO 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 EMAVP shows promising convergence rates.
Table 6 reports the results of the EMAVP 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 EMAVP with the exact method developed by Dumas et al. [3], BeamACO [14], and CA [13], and the last column is titled Heuristic, which represents the best value obtained by the algorithms developed by [9–11]. The EMAVP yields the optimal values in 19 out of 27 instances. The EMAVP yields results that are better than those of CA and Heuristic and similar to those of BeamACO. EMAVP outperforms all other algorithms on instances ( and ) and ( and ). Moreover, the EMAVP surpasses BeamACO in instances ( and and ) with respect to both solution quality and consistency. On the other hand, BeamACO is slightly better than the EMAVP in instances ( and ), ( and ), and ( and ).
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], BeamACO 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 EMAVP is faster than BeamACO and Heuristic for all of the benchmarks. However, the performances of Heuristic and the EMAVP are very similar on the benchmark set 1. Moreover, CA is faster than the EMAVP on benchmark set 3, and EMAVP outperforms CA on the other benchmark sets. Table 7 clearly shows that, in general, the proposed EMAVP is an effective algorithm and outperforms other novel algorithms in 9 out of 11 test cases in terms of computational time.

6. 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 easytocode, straightforward metaHeuristic algorithm that emulates the attractionrepulsion interactions of charged particles in analogy to Coulomb’s law in electromagnetic theory. The proposed algorithm uses two important approaches to handle timewindow 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 wellknown benchmark instances shows that the EMAVP converges to feasible regions in a search space and finds the best known or nearoptimal results. Furthermore, the EMAVP outperforms other novel metaHeuristics in terms of computational time. Future work may involve combining the VBS technique with other metaHeuristics using realcoded 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is fully supported by Uludag University Scientific Research Council under the Grant no. KUAP(M)2012/32 and is a part of Ph.D. thesis of the first author.
References
 E. K. Baker, “Technical note—an exact algorithm for the timeconstrained traveling salesman problem,” Operations Research, vol. 31, no. 5, pp. 938–945, 1983. View at: Google Scholar
 M. W. P. Savelsbergh, “Local search in routing problems with time windows,” Annals of Operations Research, vol. 4, no. 1, pp. 285–305, 1985. View at: Publisher Site  Google Scholar
 Y. Dumas, J. Desrosiers, E. Gelinas, and M. M. Solomon, “An optimal algorithm for the traveling salesman problem with time windows,” Operations Research, vol. 43, no. 2, pp. 367–371, 1995. View at: Publisher Site  Google Scholar
 A. Langevin, M. Desrochers, J. Desrosiers, S. Gélinas, and F. Soumis, “A twocommodity flow formulation for the traveling salesman and the makespan problems with time windows,” Networks, vol. 23, no. 7, pp. 631–640, 1993. View at: Google Scholar
 N. Ascheuer, M. Fischetti, and M. Grötschel, “Solving the asymmetric travelling salesman problem with time windows by branchandcut,” Mathematical Programming B, vol. 90, no. 3, pp. 475–506, 2001. View at: Google Scholar
 E. Balas and N. Simonetti, “Linear time dynamicprogramming algorithms for new classes of restricted TSPs: a computational study,” INFORMS Journal on Computing, vol. 13, no. 1, pp. 56–75, 2001. View at: Google Scholar
 F. Focacci, A. Lodi, and M. Milano, “A hybrid exact algorithm for the TSPTW,” INFORMS Journal on Computing, vol. 14, no. 4, pp. 403–417, 2002. View at: Publisher Site  Google Scholar
 G. Pesant, M. Gendreau, J.Y. Potvin, and J.M. Rousseau, “An exact constraint logic programming algorithm for the traveling salesman problem with time windows,” Transportation Science, vol. 32, no. 1, pp. 12–29, 1998. View at: Google Scholar
 W. B. Carlton and J. Wesley Barnes, “Solving the travelingsalesman problem with time windows using tabu search,” IIE Transactions, vol. 28, no. 8, pp. 617–629, 1996. View at: Google Scholar
 M. Gendreau, A. Hertz, G. Laporte, and M. Stan, “A generalized insertion heuristic for the traveling salesman problem with time windows,” Operations Research, vol. 46, no. 3, pp. 330–335, 1998. View at: Google Scholar
 R. W. Calvo, “New heuristic for the traveling salesman problem with time windows,” Transportation Science, vol. 34, no. 1, pp. 113–124, 2000. View at: Google Scholar
 R. F. Da Silva and S. Urrutia, “A General VNS heuristic for the traveling salesman problem with time windows,” Discrete Optimization, vol. 7, no. 4, pp. 203–211, 2010. View at: Publisher Site  Google Scholar
 J. W. Ohlmann and B. W. Thomas, “A compressedannealing heuristic for the traveling salesman problem with time windows,” INFORMS Journal on Computing, vol. 19, no. 1, pp. 80–90, 2007. View at: Publisher Site  Google Scholar
 M. LópezIbáñez and C. Blum, “BeamACO for the travelling salesman problem with time windows,” Computers and Operations Research, vol. 37, no. 9, pp. 1570–1583, 2010. View at: Publisher Site  Google Scholar
 Ş. İ. Birbil and S.C. Fang, “An electromagnetismlike mechanism for global optimization,” Journal of Global Optimization, vol. 25, no. 3, pp. 263–282, 2003. View at: Publisher Site  Google Scholar
 P. Wu, W.H. Yang, and N.C. Wei, “An electromagnetism algorithm of neural network analysis: an application to textile retail operation,” Journal of the Chinese Institute of Industrial Engineers, vol. 21, no. 1, pp. 59–67, 2004. View at: Google Scholar
 Ş. İ. Birbil and O. Feyzioğlu, “A global optimization method for solving fuzzy relation equations,” in Fuzzy Sets and SystemsIFSA, pp. 718–724, Springer, Berlin, Germany, 2003. View at: Google Scholar
 D. Debels, B. De Reyck, R. Leus, and M. Vanhoucke, “A hybrid scatter search/electromagnetism metaheuristic for project scheduling,” European Journal of Operational Research, vol. 169, no. 2, pp. 638–653, 2006. View at: Publisher Site  Google Scholar
 C.S. Tsou and C.H. Kao, “Multiobjective inventory control using electromagnetismlike metaheuristic,” International Journal of Production Research, vol. 46, no. 14, pp. 3859–3874, 2008. View at: Publisher Site  Google Scholar
 B. Maenhout and M. Vanhoucke, “An electromagnetic metaheuristic for the nurse scheduling problem,” Journal of Heuristics, vol. 13, no. 4, pp. 359–385, 2007. View at: Publisher Site  Google Scholar
 J.Y. Jhang and K.C. Lee, “Array pattern optimization using electromagnetismlike algorithm,” AEUInternational Journal of Electronics and Communications, vol. 63, no. 6, pp. 491–496, 2009. View at: Publisher Site  Google Scholar
 P.C. Chang, S.H. Chen, and C.Y. Fan, “A hybrid electromagnetismlike algorithm for single machine scheduling problem,” Expert Systems with Applications, vol. 36, no. 2, pp. 1259–1267, 2009. View at: Publisher Site  Google Scholar
 P. Wu, K.J. Yang, and B.Y. Huang, “A revised EMlike mechanism for solving the vehicle routing problems,” in Proceedings of the IEEE 2nd International Conference on Innovative Computing, Information and Control (ICICIC '07), Kumamoto, Japan, September 2007. View at: Publisher Site  Google Scholar
 A. Yurtkuran and E. Emel, “A new hybrid electromagnetismlike algorithm for capacitated vehicle routing problems,” Expert Systems with Applications, vol. 37, no. 4, pp. 3427–3433, 2010. View at: Publisher Site  Google Scholar
 Z. NajiAzimi, P. Toth, and L. Galli, “An electromagnetism metaheuristic for the unicost set covering problem,” European Journal of Operational Research, vol. 205, no. 2, pp. 290–300, 2010. View at: Publisher Site  Google Scholar
 A. GarciaVilloria and R. P. Moreno, “Solving the response time variability problem by means of the electromagnetismlike mechanism,” International Journal of Production Research, vol. 48, no. 22, pp. 6701–6714, 2010. View at: Publisher Site  Google Scholar
 X. Guan, X. Dai, and J. Li, “Revised electromagnetismlike mechanism for flow path design of unidirectional AGV systems,” International Journal of Production Research, vol. 49, no. 2, pp. 401–429, 2011. View at: Publisher Site  Google Scholar
 A. Jamili, M. A. Shafia, and R. TavakkoliMoghaddam, “A hybridization of simulated annealing and electromagnetismlike mechanism for a periodic job shop scheduling problem,” Expert Systems with Applications, vol. 38, no. 5, pp. 5895–5901, 2011. View at: Publisher Site  Google Scholar
 E. Cuevas, D. Oliva, D. Zaldivar, M. PérezCisneros, and H. Sossa, “Circle detection using electromagnetism optimization,” Information Sciences, vol. 182, no. 1, pp. 40–55, 2012. View at: Publisher Site  Google Scholar
 C.T. Su and H.C. Lin, “Applying electromagnetismlike mechanism for feature selection,” Information Sciences, vol. 181, no. 5, pp. 972–986, 2011. View at: Publisher Site  Google Scholar
 C.H. Lee and F.K. Chang, “Fractionalorder PID controller optimization via improved electromagnetismlike algorithm,” Expert Systems with Applications, vol. 37, no. 12, pp. 8871–8878, 2010. View at: Publisher Site  Google Scholar
 P. Godinho and F. G. Branco, “Adaptive policies for multimode project scheduling under uncertainty,” European Journal of Operational Research, vol. 216, no. 3, pp. 553–562, 2012. View at: Publisher Site  Google Scholar
 C.H. Lee and Y.C. Lee, “Nonlinear systems design by a novel fuzzy neural system via hybridization of electromagnetismlike mechanism and particle swarm optimisation algorithms,” Information Sciences, vol. 186, no. 1, pp. 59–72, 2012. View at: Publisher Site  Google Scholar
 J. C. Bean, “Genetic algorithms and random keys for sequencing and optimization,” ORSA Journal on Computing, vol. 6, no. 2, pp. 154–160, 1994. View at: Publisher Site  Google Scholar
 J.Y. Potvin and S. Bengio, “The vehicle routing problem with time windows part II: genetic search,” INFORMS Journal on Computing, vol. 8, no. 2, pp. 165–172, 1996. View at: Google Scholar
 M. M. Solomon, “Algorithms for the vehicle routing and scheduling problems with time window constraints,” Operations Research, vol. 35, no. 2, pp. 254–265, 1987. View at: Google Scholar
 J. J. Dongarra, Performance of Various Computers Using Standard Linear Equations Software, Department of Computer Science, University of Tennessee, 1989.
Copyright
Copyright © 2014 Alkın Yurtkuran and Erdal Emel. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.