DENG Ye,ZHU Wanhong,*,LI Hongwei,and ZHENG Yonghui

1.College of Field Engineering,Army Engineering University of PLA,Nanjing 210001,China;

2.Department of Fire Control,Armored Force Institute of PLA,Bengbu 214500,China

1.Introduction

Since Dantzig and Ramser[1],Clarke and Wright[2]proposed the vehicle routing problem(VRP),it has been a rich topic for researchers and practitioners.Many variants of the classic VRP are motivated by the complexity of supply chain management and the physical delivery of goods services in order to perform a realistic optimization.Taking into consideration the time windows of the customers and the depot(s)extends the problem into the vehicle routing problem with time windows(VRPTW).There are also many possible extensions of the VRPTW,such as pickup and delivery[3,4],VRPTW with backhauls[5],and period routing[6,7].No matter for VRP or VRPTW or their variants,most research focuses on static or deterministic vehicle routing in which all information is known at the time of planning.However,in the real-life world,travel speed or travel time may be in constant flux due to traffic congestion,weather conditions,moving targets,mobile obstacles,or endogenously whenever the decision maker can set the vehicles’speeds[8].Hence,the time dependent vehicle routing problem with time windows(TDVRPTW)model gets more and more attention.

Malandraki and Daskin[9] first addressed the time dependent vehicle routing problem(TDVRP)as a real strict time depending issue and proposed a mixed integer linear programming solution to TDVRPTW.After that,Ichoua et al.[10]proposed a step function for the speed and a piecewise linear function of time distribution that satisfied the first in first out(FIFO)principle[11]which was defined by Ahn et al.Most research after Ichoua adopted the speed model to study TDVRP or its variants and concentrated on proposing high performance algorithms.Fleischmann et al.[12]utilized a route construction savings method and an insertion method to solve uncapacitated TDVRP with and without time windows.Albiach et al.[13]proposed a genetic algorithm to solve the time dependent traveling salesman problem(TDTSP)with time windows(TDTSPTW).In the last decades,many algorithms for TDVRPTW were also explored.There are metahueristic algorithms[14–16],some hueristics[17–19]and hyper hueristics[20].More detailed information about TDVRPTW is specified in the review[8]by Michel Gendreau et al.

Recently,some research began to consider more complex situations.Behrouz and Alireza[21]studied the TDVRPTW with multi-depot and heterogeneous flee model.Setak et al.[22]considered a new TDVRPTW in multigraph where existence of more than one edge between the nodes was allowed.Besides,some other considerations like stochastic nature[23,24],congestion charge related to time periods[25]were added to TDVRPTW.All these work make TDVRPTW become one of the most attractive research topics in operational research and public logistics fields.

Ant colony optimization(ACO)has been widely used and developed to solve different types of NP-hard combinatorial optimization problems and performed well.It was first introduced by Dorigo et al.[26],inspired by the collaborative behavior of ant colonies in the real world.Since then,a number of improved ACO algorithms have been introduced such as ant colony system(ACS)by Gambardella et al.[27,28],max-min ant system by Stützle et al.[29],ASrank by Bullnheimer et al.[30],and a hyper-cube framework by Blum et al.[31].For VRPTW model,Gambardella et al.[32]presented a multiple ant colony system(MACS-VRPTW),which shows competitive results of most of Solomon benchmark instances with previousalgorithms.Then Donati et al.respectively applied the ACS algorithm[33]and multi ACS(MACS)algorithm[33–35]for the TDVRPTW in the real road network.After that,Balseiro et al.[16]found the ant colony algorithm for VRPTW had a weakness in that ants tend to create infeasible solutions with unserved clients,so he presented an aggressive insertion heuristics relying on the minimum delay metric which is hybridized with the MACS algorithm to fill the gap.Duan and Sun considered a stochastic TDVRPTW and used ACO to solve the model[24].Some other algorithms are also proven to be effective in TDVRPTW[36,37],but ACO seems more competitive and stable.However,just like other algorithms,the weakness of ACO always exists,mainly referred to easily trapping into local optimum,the slow convergence speed,the tendency to create infeasible solutions with unserved clients,etc.Therefore,it is necessary to combine ACO with some classic heuristics or other meta-heuristics[3].

Solomon et al. first proposed a successful sequential insertion heuristic which is called I1 for VRPTW[38],then Potvin et al.implemented it with a parallel version[39].Atkinson et al.proposed a greedy minimization function of look-ahead solution approach[40].Based on that,Ioannou et al.proposed new criteria for customer selection and insertion[41].Lau et al.proposed a Tabu search algorithm with a holding list containing the unserved clients for the m-VRPTW and defined an extended neighborhood space that allows relocations and exchanges with the holding list[42].Lim and Zhang presented three improvement approaches for their Tabu search algorithm[43].Based on the the three strategies,Balseiro et al.proposed an insertion hueristics that combines the local search’s algorithmic framework with an insertion procedure[16].

As for local search procedure,it can be classified into two categories which are edge-exchange algorithms and node-exchange algorithms.For edge-exchange algorithms,there are k-opt[44,45],Or-opt[46],2-opt∗[39],λ-interchange[47],CROSS-exchange[48],etc.For nodeexchange algorithms,there are relocate[45],exchange[49],GENI-exchange[50],etc.

To solve TDVRPTW in this paper,we propose a multitype ant system(MTAS)algorithm structure and an adaptive pheromone updating strategy to make the optimization procedure more efficient,which will be discussed in Section 3.Besides,three improvement heuristics are specified in Section 4,namely a nearest neighbor selection(NNS)mechanism,an insertion local search procedure and a local optimization procedure.

2.Problem model and formulation

Traditional flow-arc formulation[51]can be used to describe TDVRPTW as follows.Let G=(V,D)be a complete directed graph with nodes set V=(v0,...,vn+1)and edge set D={(vi,vj)|vi,vj∈V,i/=j}.Let K be a vehicle set and each vehicle has the same load capacity Q.In nodes set,nodes v0and vn+1denote the depots while others represent clients and consist of the client set C.Each node has a demand qi≥0,a service time si≥0,and a service time window[ei,li].In edge set,each edge(vi,vj)has a distance dij≥ 0.If we denote the vehicle’s arrival time at client i as Ati,the departure time as Dti,the travel time tij(dij,f vij(Dtj))or tij(Dti)from node i to node j can be a continuous piecewise linear function of the distance dijand the speed distribution f vij(Dti)which is a discrete step function of departure time Dtifrom node i.

To construct a formulation for TDVRPTW,some assumptions and constraints must hold:

(i)Any individual client is possible to be served from the depot and visited exactly once by one vehicle.

(ii)There are enough vehicles for service to ensure that every client can be served.

(iii)Quantity constraint:the demand quantity qiof client i must be completely satisfied once client i is visited and the load capacity of vehicle k cannot be violated.

(iv)Time constraints:the arrival time at client i Atimust be earlier than the time window upper bound li,however,if Atiis earlier than ei,the vehicle must wait until eistarts to serve.Client i will be left immediately once it has been served with time si.When vehicle k leaves node i respectively at time Dti1 and Dti2(Dti1＜D ti2),the arrival time at node j is Atj1=Dti1+tij(Dti1)and Atj2=D ti2+tij(Dti2)mustsatisfy Atj1≤Atj2(FIFO property),which means any differentiable segment’s slopeof the piecewise function tij(Dti)satisfiesIn fact,using discrete step function for speed always meets this constraint.An example of a speed function f vij(D ti)and travel time function tij(Dti)between nodes i and j is presented in Fig.1.The result of the travel time function is calculated based on the speed function in Fig.1(a)and an assumption of the edge distance dij=50 units.The detail of the calculation procedure should refer to Ichoua et al.[10].

Fig.1 Example of speed function and travel time function

In this paper,we establish the hierarchical structure of multi-objectives for the TDVRPTW.It involves the minimization of the following objectives in order of priority:i)total number of vehicles used;ii)total travel distances and total arriving time including travel time,service time and waiting time.

The objectives are arranged in decreasing order of i)and ii)because the vehicle number cost is much more important than the time cost or distance cost.In the real world,time cost may be comparatively more important than distance cost in TDVRPTW.This approach is first used by Lau et al.[42]and its obvious superiority to the“big M”approach is that if a solution is better than another one on a higher hierarchy level,the lower level objective costs need not to compute whereas“big M”approach computes all the time.

(i)A binary decision variable:

(ii)A positive number decision variableindicates start time of service for customer i served by vehicle k.

Finally,the TDVRPTW can be formulated as a multiobjective mixed integer linear programming(MMILP)model.

Objectives:

Subject to

The hierarchical structures of multi-objectives are defined by(1)and(2).Constraints are defined by follows:(3)capacity of vehicle k cannot be violated;(4)every client must be visited and only once;(5)a vehicle arriving at a client must also depart from the same one;(6)each vehicle must start and end at the depot;(7)if a vehicle starts to travel,it must leave and return to the depot only once;(8)service start time should fall within the range in the client’s time window;(9)service start time must allow for travel time between each paire of clients.

3.MTAS for TDVRPTW

The famous MACS-VRPTW[52]consists of two artificial ant colonies,each one dealing with one of the objectives of the optimization:the first colony is named ACS-VEI and deals with number of active vehicles minimization while the second colony is named ACS-TIME and deals with distance or time minimization.The two colonies co-operate by sharing the best solution through pheromone updating.In this paper,based on the classic MACS-VRPTW,we propose a new MACS which consists of two types of ant colony for the TDVRPTW,i.e.,ACS-VEI and MMASTIME.Specifically,we use the new MMAS-TIME procedure instead of ACS-TIME procedure in MACS-VRPTW to optimize solutions,and apply an adaptive pheromone updating strategy to achieve a balance between good scouting performance and fast convergence rate(details in Section 3.3).Besides,we introduce an NNS mechanism to improve the efficiency of insertion hueristics based on Balseiro et al.[16]and Lim et al.[43](details in Section 4.1).

3.1 Design of MTAS-TDVRPTW

The main procedure of MTAS-TDVRPTW is shown in the following steps:

Step 1Initialize the number of vehicles v with function#vehicles(ψ0),the global best solution ψgband a solution ψVEIwith v−1 vehicles.

Step 2Activate ACS-VEI with v −1 vehicles,if ψVEIis feasible,then#new solution(v−1).

Step 3Activate MMAS-TIME with v vehicles,if ψTIMEis better than ψgb,then

Step 4Repeat from Step 2 until the termination criterion is met and output ψgb.

At the beginning,we use a traditional constructive heuristics Solomon I1[38],a nearest neighbor heuristic algorithm called#new solution to find an initial feasible solution ψ0and an infeasible solution ψVEIwith v − 1 vehicles.Then the function#vehicles returns the v active vehicles in the solution ψ0which is taken as a temporarily global best solution ψgb.After that,two ant colonies ACS-VEI and MMAS-TIME are activated simultaneously to search better solutions.Once a feasible solution ψVEIis found,the global best solution ψgbis updated by it and both colonies are restarted with reduced number of vehicles,and if a better solution ψTIMEis found,though ψgbis replaced by it,the number of vehicles will remain unchanged.Considering the hierarchical structure of multiobjectives,a feasible solution ψVEIwith smaller vehicles number always has priority over a better solution ψTIMEwith more vehicles.The algorithm keeps repeating until it meets a termination criterion.The criterion can be a maximum number of iterations,a total computing time,an acceptable solution value,a confirmation of the algorithm showing stagnation behavior or any other given ones.

The proposed ACS approach uses two pheromone trails to manage the knowledge of ants,weak and strong pheromone trails.The former is updated only locally,while the latter is updated only globally.For searching a feasible solution with fewer vehicles,this approach performs well with some insertion operators,however,for searching a better solution in the long term,MMAS is proven to be more competitive[53].In practice,we can also use MMAS to replace the former ACS,the better solution is recognized as one with more clients served in routs than the best-so-far solution,but comparison results show that it performs worse than ACS as it easily attempts to getlocal optimum.More detailed comparison results of using ACS and MMAS respectively for searching fewer vehicles solutions and better solutions with constant vehicles number are shown in Section 5.

3.2 Design of ACS-VEI

The procedure of ACS-VEI(s)is shown in the following steps:

Step 1Initialize pheromone and data structres using s vehicles.

Step 2Construct solutions ψafor all ants with function#ant construct solution(s)and updating local pheromone.

Step 3If unserved clients in ψaexist,then perform Insertion local search to insert these clients and I Ni←I Ni+1.

Step 4If unserved(ψa) ＜ unserved(ψVEI)or unserved(ψa)= unserved(ψVEI)&M D(ψa) ＜M D(ψVEI),then,perform global pheromone double updating with ψgband ψVEI.

Step 5Repeat from Step 2 until the termination criterion is met and output ψVEI.

ACS-VEI procedure uses a subprogram#ant construct solution to construct unfeasible or feasible solution ψa.If the number of unserved clients turns to 0,ACS-VEI has successfully found a feasible solution with s vehicles and update global best solution ψgbin main procedure,otherwise,another subprogram#Insertion-Local Search is activated to try to insert unserved clients directly or after local searching.Here,not all the edges are tested to insert a chosen unserved client,we use an NNS mechanism to quickly preclude some edges with low insertion probability and sort the left edges in descend insertion probability at the same time.The colony maintains counters{INi}i∈Cis used to record the times of client i that has not been served in the entire searching procedure[52].These counters can make clients unserved more likely to be served in next iteration.A better solution is found in two conditions:i)a solution with more served clients,ii)a solution with the number of served clients is the same with most-sofar solution and with lower minimum delay metric(MD)[54].Finally the pheromone matrix will be updated with solution ψgband ψVEIas in(10)and(11).When a better solution ψVEIis found.The double global pheromone updating can prevent the algorithm from searching stagnation.If we also adopt MMAS to search feasible solutions,i.e.,MMAS-VEI,the algorithm will drop into local optimum more easily than ACS-VEI as shown in experiments in Section 5.

3.3 Design of MMAS-TIME

MMAS was first proposed by Stützle et al.[53].The procedure of MMAS-TIME(s)is shown in the following steps:

Step 1Initialize pheromone and data structres using s vehicles.

Step 2Construct solutions ψafor all ants with function#ant construct solution(s)but without updating local pheromone.

Step 3If unserved clients in ψaexist,then perform Insertion local search to insert these clients.

Step 4If unserved(ψa)=0 and ψais worse than ψTIME,then perform Local optimization procedure,else if ψais better than ψTIME,then

Step 5Perform global pheromone double updating with ψgbor ψTIMEby adaptive strategy,limit pheromone trail between τminand τmaxand smooth pheromone trails.

Step 6Repeat from Step 2 until the termination criterion is met and output ψTIME.

The program#ant construct solution is also used here but without local pheromone updating in comparison with that in AVS-VEI.If ants find solutions not feasible,the#ant construct solution is activated again like in ACS-VEI to insert clients.After all ants have constructed their solutions,the global updating is done as

where ρ is the pheromone evaporation rate which is a decimal greater than 0,ψbestdenotes either the global-best solution ψgbor the time-local-best solution ψTIME,corresponding to its total times value ψbest.

According to the research of Stützle[55],when pheromone updates are always performed by the bestso-far ant,the search focuses very quickly around ψgb,whereas by the iteration-best ant it becomes less directed and needs a larger search space.Hence,we use an adaptive strategy by gradually increasing the frequency of the global-best tour ψgbto be chosen for the pheromone updating as

where Piterdenotes the probability of using the globalbest tour for pheromone updating and iter denotes the number of iterations.In each iteration,if Piteris greater than a random number which is between 0 and 1,the global-best ant will perform pheromone updating,otherwise,the time-local-best ant will do this job.The exponent is set to an adaptive value λ(itermax− iter)to increase Piterwith nonlinear scaling that can make the chosen probability of global-best tour smaller at the very beginning and bigger near the end.λ is a proportional control factor.When λ is set to 0,Piterwill be only related to objective values,and λ is set to 1,Piterwill be strongly impacted by iteration simultaneously.

The differences between Piterwithout exponent and with adaptive exponent in three convergence rates are presented in Fig.2.

Fig.2 Diagram of different pheromone updating strategies

It shows that,when the convergence speed becomes faster,the change of Piteralso becomes faster,which has a strong connection with the objective value.This self related adaptive pheromone updating strategy is comparatively more reasonable than those static strategies which only update with the iteration-best solution or the global best solution,which lies in that it finds a balance between good scouting performance and fast convergence rate.Experiments results can be seen in Section 5.

In Fig.2,the dotted lines denote Piterwithout the exponent(or exponent equals 1)and the solid lines denote Piterwith an adaptive exponent value of itermax−iter(λ =1).

In order to avoid search stagnation,lower and upper limits τminand τmaxwith the possible pheromone values on any edge are imposed.The value of τmaxis set to 1/(ρTgb)and

where av g is the average number of different choicesavailable to an ant at each step when it constructs a solution[29].In addition to trail limits,another useful mechanism called smoothing of trails[56]is adopted in case of the stagnation of search as(15)where c is a constant value set by manual and reflects the reinitialized to what extent when the algorithm meets a stagnation.If c is set toward 1,the pheromone will be reinitialized completely.

3.4 Solution construct procedure

The#ant construct solution subprogram is used in both ACS-VEI and MMAS-TIME.It describes how an ant constructs a solution in detail.The ant chooses unserved nodes to travel one by one according to the result of constraints checking by C ons C heck and nodes chosen rules by Decision Rule.The C ons C heck procedure includes quantity constraint checking and time constraints checking referring to Section 2.All nodes which have not been visited but can be visited from node i are stored in nodes set Ni.The Decision Rule is the procedure employed by ants to choose the next node to visit as(16).In our Decision Rule,artificial ants make decisions autonomously depending on a component probability of three parts:i)the pheromone trail τij,which is used to record ants’pheromone of movement information;ii)the attractiveness ηij,which is calculated as(17)is used to attract ants to explore possible new solutions;iii)the saving hueristics value uij,which is calculated as(18),is used to strengthen the exploitation efficiency in time dependent situation.The exponents α,β and γ reflect the relative importance of τij, ηijand μij.Equation(17)successfully puts three factors into a combination which is the distances and travel times between two nodes,time windows width and nodes unvisited times denoted by I Nj.The I Njwill be set to 0 in MMAS-TIME searching procedure.Here,we also adopt the probabilistic rule[52]to increase solutions diversity.After each movement,ants always leave information with pheromone and the local updating procedure is performed.However,the local updating procedure is performed only in ACS but not in MMAS because ants of MMAS only do global updating.The parameterτ0is the pheromone increment which can be set as τ0=1/(nT0),where n is the sum of nodes and T0is the second objective value of initial solution ψ0which has been mentioned in Section 3.1.

4.Improvement heuristics

4.1 Nearest neighbor selection mechanism

The NNS mechanism we propose is a heuristic strategy to improve the efficiency of insertion procedure.It is like the Nearest-Neighbor List in traveling salesman problem(TSP).In TSP,the nearest clients are only selected by a constant radius in distance.However,the NNS mechanism concentrates on sorting the high possible insertion places around an unserved client by considering both the distance and time factors.It consists of two procedures,one is selecting insertion places with descending insertion priority order,called P re-selection,and the other is selecting unserved clients with ascending potential insertion probabilities order to insert,called Post-selection.The potential insertion places are those infeasible but have potential to be feasible,especially after the local search operations.These two procedures together make the Insertion local search procedure become more efficient.

(i)P re-selection

An unserved client v must be inserted to a rational place,supposing the insertion place is after client i,that as near client i as possible not only in distance but also in time factors such as waiting time,travel time and time windows.Solomon[38]once proposed three metrics and a total value with different weights to measure the closeness of customers as

where δ1+ δ2+ δ3=1,Tiv=[ev− (ei+si)]and uiv=lv−[ei+si+tiv(ei+si)].

However,the exact value of weights cannot be determined always and the waiting time is not considered either.Hence,we propose a multiple formative function with consideration of waiting time as

where the minimum waiting time denotes wiv,wiv=max{ev−[li+si+tiv(li+si)]}.

If uiv＜0,the client v is impossible to insert after client i and then results in civ＜0 which means the insertion place after client i is unfeasible and excludes it in advance.However,this metric cannot make sure the insertion places are feasible because the real departure time from a client i is not taken into calculation and whether an insertion is feasible must depend on the result of feasiblity checking under quantity and time constraints in the later proceedings.Hence,the nearest neighbor insertion places,a subset of all clients for client v,i.e.,N earv,can be built before the algorithm’s iterations begin.

(ii)Post-selection

If an insertion place is feasible under constraints,an unserved client can be inserted successfully.However,there are always some clients which remain unserved and cannot be inserted anyway.The Post-selection is then activated to sort the unserved clients in ascending insertion possibilities by the minimum delay metric M D(v).The minimum delay metric M D(v)was first proposed by Homberger and Gehring[54]and further extended by Bent[57],Lim[58]and Balseiro[16],in order to measure how hard it is to insert an unserved client v into a sub-route r,i.e.,r=(v0,v1,...,vm)where v0and vmdenote the depot.Considering the insertion of client v in the sub-route r between client viand client vi+1,M D(v)is calculated by adding the three time shifts md1r(vi,v),md2r(vi,v)and md3r(vi,v)as(21),where md1r(vi,v)indicates the violation of time window’s upper bound lvof client v,md2r(vi,v)indicates the time of arriving later than permitted latest arrival time of client vi+1and md3r(vi,v)indicates the time window’s lower bound evlater than the permitted latest arrival time of client v.

where

Lavidenotes the permitted latest arrival time of client vi,which can be calculated as(22).The latest arrival time of the last client vmis its time window’s upper bound,i.e.,La vm=l vm.

gives the travel time from client i to client j based on the arrival time Atjat client j.Ichoua[10]also proposed a procedure to calculate the departure time which is equal to calculating the travel time.If we only want to measure how hard it is to insert client v to a sub-route r,we can use the M D(v,r)as(23).And if we want to measure the potential insertion ability of a solution ψ,we can use M D(ψ)as(24).

Besides,if we have obtained the nearest neighbor clients set N earvfor a chosen unserved client v by using P reselection,we only need to calculate the M D within the clients set N earvas

After calculating the minimum delay values of each subroute r,i.e.,M D(v,r),the best potential insertion place for an unserved client will be found and the corresponding minimum delay M D(v)will also be calculated.We sort these potential insertion places by the value of M D(v)in a descending order which logically makes the probabilities of successful insertion in an ascending order.

4.2 Insertion local search procedure

The Insertion local search procedure is activated both in ACS-VEI and MMAS-TIME procedures and it means to improve the algorithm’s ability in constructing feasible solutions.It applies the NNS mechanism to improve the efficiency of the searching procedure.Besides,it also uses a data structure,called ejection pool(EP),proposed by Lim[43],to hold unserved clients temporarily and exchange some served clients with unserved ones to improve insertion probabilities.

The Insertion local search procedure combines with two main stages of insertion.In the first stage,it uses the NNS mechanism of P re-selection to reduce the searching scope of an unserved client v as the set of N earvand use Insertion Hueristics proposed by Solomon[38]to insert clients directly as

whereβ1+β2=1,β1≥ 0,β2≥ 0.The original start time of service for client j isand the start time of service for client j after insertion of client v is S tj=max{Dti+tij(Dti),ej}.

In the second stage,it uses the NNS mechanism of Post-selection to sort those still unserved clients in the set EP with the ascending insertion probabilities order and use local search operators to explore potential new insertion places.While an insertion of client v makes a sub-route infeasible,some served clients,except client v,must be continuously taken off the sub-route and added into the EP until the target sub-route becomes feasible.The number of clients removed from route to EP will be kept as less as possible until all clients have been inserted during limited iterations.The metric of deciding which the served client is removed best is based on the“Kick Saving”proposed by Lim[43].The“Kick Saving”is presented as

where θ1+ θ2=1,θ1≥ 0,θ2≥ 0 and

4.3 Local optimization procedure

In Section 4.2,the local search operators are applied with the insertion procedure to search feasible solutions.After the feasible solutions have been found,the MMAS-TIME is activated to optimize the solution quality with the Local optimization procedure.It can be obtained by iteratively moving to a neighbor solution by using some local search operators which has been already applied in the Insertion local search procedure.However,the objective of these local search operators in the Insertion local search procedure is to decrease the minimum delay or increase the maximal free time.In the Local optimization procedure,it turns to decrease the objective of total times or the total distances as(2).In this paper,we apply Exchange,Relocate,2-opt,Or-opt,2-opt∗and Cross to improve the solution quality.In order to decrease the computational time,a moving rule is proposed that the Exchange and Relocate operators are applied in the early iterations,and other operators are applied when no more improvement are obtained.All operators are presented in the Figs.3–6.

Fig.3 Diagram of Exchange operator

Fig.4 Diagram of Relocate operator

Fig.5 Diagram of 2-opt and Or-opt operator

Fig.6 Diagram of 2-opt∗and cross operator

5.Experiments and results

In this section,we will first introduce the classical benchmark instances of Solomon’s 100 clients Euclidean problems for VRPTW and the extension of Solomon’s instances proposed by Figliozzi[18]for TDVRPTW.Then,we test the efficiency of the two types of ant colony in the MTAS algorithm based on the VRPTW instances.Finally,we adopt MTAS-TDVRPTW to solve the instances of TDVRPTW and compare the results with that of Figliozzi[18].

5.1 Benchmark instances

As in the constant speed case,the well-known 56 Solomon benchmark instances for VRPTW are based on six groups of problem instances with 100 clients.The six problem classes are named C1,C2,R1,R2,RC1,and RC2.Clients locations are clustered in groups in sets C1 and C2,randomly generated by a uniform distribution in sets R1 and R2,and mixed with randomly generated and clustered customer locations in sets RC1 and RC2.Instances sets R1,C1,and RC1 have a shorter scheduling horizon,tighter time windows,and lower vehicle capacity,thus permitting fewer clients to be served per route than instances sets R2,C2,and RC2 respectively.In these instances the primary objective is to minimize the number of vehicles and the secondary objective is to minimize traveled distance.

As in the varying speed case,Ichoua et al. first provided a time-dependent extension of Solomon’s instances[10].However,because the instances were solved by considering the time windows as “soft”,results showed that when hard time windows constraint was imposed strictly,some instances became infeasible and made comparisons impossible.Figliozzi said the only approach to leverage the existing Solomon instances to evaluate time dependent problems was to introduce speed variability with higher travel speeds[18].And the new test instances he proposed were readily replicable for time-dependent vehicle routing problems with hard time windows.The instances capture the typical speed variations of congested urban settings which generally include slower travel speeds during morning and evening rush hours.The maximum ratio in travel speeds is set from 2.5 to 1 which is realistic to assume that in an urban freeway the free- flow or uncongested travel speed is around 60 miles per hour and with heavy congestion the speed drops to 24 miles per hour.Another difference from Ichoua’s instance is that the speed distribution is assigned to all edges to facilitate reproducibility[10].

In Figliozzi’s instance,the entire time horizon,i.e.,the depot working time[e0,l0],is divided into five time periods of equal durations:[0,0.2l0),[0.2l0,0.4l0),[0.4l0,0.6l0),[0.6l0,0.8l0),[0.8l0,l0].

Four distinct cases are also proposed corresponding with different time periods.In each case,the vehicle would travel an extra 25%,50%,and 75%more for speeds TD1,TD2,and TD3 respectively than in the original Solomon instances.Here,we only adopt the case of type(a)as an example,other cases can be seen in Figliozzi et al.[18].The corresponding travel speeds are as follows:

In case of this travel speed type,the vehicles leave the depot very early in the day and travel without congestion during the first 1/5 of the workday before morning congestion taking place.The following 1/5 of the workday coincides with morning congestion followed by the somewhat less congested midday period.Another 1/5 of the workday coincides with evening congestion and the last 1/5 of the workday takes place after congestion has subsided.

All the test work is completed on a laptop Lenovo G550,with Windows 7 Ultimate editions,an Intel Core(TM)2 Duo CPU 2.10 GHz and 5.00 GB RAM.The algorithm is coded and compiled in Matlab R2014b.For the MTAS algorithm structure,based on an initial solution by Solomon I1 heuristics,we test the ACS and MMAS algorithms for solutions construction and optimization with Solomon’s VRPTW instances.For the test of TDVRPTW benchmark instances,we adopt MTAS with all improvement strategies and make a general comparison between the results of our algorithm and Figliozzi[18].

5.2 Experiments with two types of ants colony in MTAS

In order to test the efficiency of two types of ants colony in the MTAS.Two experiments are conducted with VRPTW instances.The first experiment tests the algorithm’s ability of searching feasible solutions,using ACS and MMAS respectively and the MMAS with two pheromone updating strategies,one uses only the best-so-far solution and another uses an adaptive pheromone updating strategy as(14).The objective is to minimize the number of vehicles as(1).Six instances are randomly picked out to test.For each instance,the algorithm is run for 10 times.Results are presented in Table 1,where MNV is the minimum number of vehicles.Among all solutions the algorithm has searched during 10 times running tests,NGap is the average number gap more than MNV.The second experiment tests the algorithm’s ability of optimizing candidate solutions which are initialized by Solomon I1 hueristics,using ACS and MMAS algorithms respectively as well.The objective is to minimize the total distances as(2)with C d=1,C t=0.Results are presented in Table 2.BTD is the best objective value of total distance during 10 test periods and NV is the number of vehicles corresponding to the solution with BTD,Gap is the average gap in distance larger than BTD.

Table 1 Comparing the results of different strategies for feasible solutions construction

Table 2 Comparing the results of different strategies for candidate solutions optimization

However, we only use the classic Insertion Hueristics as(26)to make solutions feasible during the two experiments,the NNS mechanism and the local search procedure are all unused here for the reason of avoiding interference.Refer to Gambardella[52]and Stützle[53],the essential parameters in this paper are set as follows:itermax=100,m=20 ants,α =1,β =2,γ =1,q0=0.8.In the ACS-VEI algorithm,the pheromone evaporation rate ρ=0.1 and in the MMAS-VEI algorithm it turns out ρ=0.01,and if the adaptive strategy mentioned in Section 3.3 is used in MMAS-VEI,the specific parameters are set as follows:λ=0.5,c=0.5.

Results in Table 1 show that,the best algorithm for searching feasible solutions with fewer vehicles is ACSVEI,though it does not get the currently best known vehicles number(Column 2 in Table 1)for some instances,like R108,RC104 and RC106,it shows acceptable results considering it performs well with a small number of iterations and without any local search procedure.The MMASVEI with an adaptive strategy is slightly better than the MMAS-VEI in the MNV,however,it still performs worse than ACS-VEI.

Results in Table 2 show that,the best algorithm for optimizing initial solutions is the MMAS-TIME with the adaptive strategy,for its shortest TD and relatively lower Gap.On the contrary with the first experiment,this test does not adopt any special way to decrease the number of vehicles in candidate solutions.However,the optimization algorithms can still decrease the vehicles number to some extent as shown in value of NV and the three algorithms all get the same value of NV for each instance.

Based on the results of the two experiments,we make the decision to use the ACS-VEI algorithm to decrease the total NV of solutions and use MMAS-TIME algorithm to decrease the total distances or total times of solutions,and these two algorithms together comprise the main part of the MTAS algorithm structure.We also apply this algorithm structure to test the TDVRPTW instances in the last part of this section.

5.3 Results of TDVRPTW benchmark instances

In this experiment,all improvement strategies are applied in the MTAS algorithm for TDVRPTW benchmark instances proposed by Figliozzi[18].Each instance is run five times for 600s and the best result is kept and shown in Tables 3–5.Figliozzi also proposed an algorithm called iterative route construction and improvement(IRCI)to solve all 56 TDVRPTW instances and the results are also presented in Tables3–5.However,the difference from former experiments is that the algorithms’primary objective is to minimize the NV as(1),the secondary objective is to minimize the total time as(2)which sets C d=0,C t=1.The former experiments in[18]made the secondary objective as a combination of total distance and total time that always result in different solutions obviously influenced by parameter C d and C t.

Table 3 Comparison of the results in travel speed model TD1

Table 4 Comparison of the results in travel speed model TD2

Table 5 Comparison of the results in travel speed model TD3

In Tables 3–5,for all instances of each class,ANV indicates the average NV,ATT indicates the average total travel time,ATDindicates the average total distance and AV indicates the average travel speed.The MTAS algorithm shows more competitive results than IRCI where the ANV values are matched with or outperform the latter and the specific NV for each instance is listed in Table 6.

Table 6 Comparing the minimal NV in the constant speed model(VRPTW)with the varying speed model(TDVRPTW)for each instance

We notice a fact that when the general speed rate becomes higher(TD1 to TD3),the whole NV will become smaller.The ATT and ATD values of most instance classes obtained by our algorithm are lower than IRCI though some of the ATD values are slightly higher in a few in-stance classes,considering the time cost weights much more than the distance cost in the real world,the results shown here should be more reasonable.In general,the average AV values by our algorithm are 1.22,1.43 and 1.62,which have improved by 1.21%,2.35%and 4.27%separately for each speed model in comparison with the IRCI algorithm.

In Table 6,the NV for each instance in the second column is the best ever known result in the constant speed model(VRPTW)and in the third to the fifth column are the bestresults of three varying speed model sever obtained during five run times for 600s by our algorithm.The results of C class instances between two models have no difference(NV of C1 class is 10 and NV of C2 class is 3),hence,we only present the results of R and RC class instances for short.

6.Conclusions

In this paper,a variation of the VRPTW,the time dependent vehicle routing problem has been studied.This model is common to see in the real world especially in the urban logistics field where the times factor is deemed to be very important and traffic flow fluctuation takes place all the time.To solve this problem,we propose an MTAS which is combined with the ACS and the MMAS,absorbing the merits of both,in solutions construction and optimization procedures.Besides,an efficiency improving strategy is proposed called NNS mechanism for insertion procedure,combining with some local search operators,making the MTAS algorithm structure searching feasible solutions with fewer vehicles more efficient.And an adaptive pheromone updating strategy is adopted by MMAS to avoid searching stagnation in solutions optimization which makes the MTAS algorithm structure search better solutions in total distances or total times effectively.

Experiments on Solomon VRPTW instances set and TDVRPTW instances set are performed with the same parameters.Results confirm the benefits of using all mentioned strategies in the MTAS algorithm for all experiments.In the constant speed case,the MTAS algorithm outperforms the classic MACS algorithm during the same computational time.In the varying speed case,the MTAS algorithm outperforms the IRCI by Figliozzi in the total travel time values and vehicle number values.The experiments for TDVRPTW have also shown a fact that,the higher the average speed is,the more vehicles will be saved.The future efforts may develop the ACO to solve more complicated variants of the VRP model by using the logistics data or traffic information from the intelligent transportation system.

[1]DANTZIG G B,RAMSER J H.The truck dispatching problem.Management Science,1959,6(1):80–91.

[2]CLARKE G,WRIGHT J W.Scheduling of vehicles from a central depot to a number of delivery points.Operations Research,1964,12(4):568–581.

[3]AVCI M,TOPALOGLU S.An adaptive local search algorithm for vehicle routing problem with simultaneous and mixed pickups and deliveries.Computers&Industrial Engineering,2015,83:15–29.

[4]ZHANG T,CHAOVALITWONGSE W A,ZHANG Y.Integrated ant colony and tabu search approach for time dependent vehicle routing problems with simultaneous pickup and delivery.Journal of Combinatorial Optimization,2014,28(1):288–309.

[5]BELLOSO J,JUAN A A,MARTINEZ E,et al.A biasedrandomized metaheuristic for the vehicle routing problem with clustered and mixed backhauls.Networks,2017,69(3):241–255.

[6]LUO Z,QIN H,CHE C H,et al.On service consistency in multi-period vehicle routing.European Journal of Operational Research,2015,243(3):731–744.

[7]CHENG C,QI M,WANG X,et al.Multi-period inventory routing problem under carbon emission regulations.International Journal of Production Economics,2016,182:263–275.

[8]GENDREAU M,GHIANI G,GUERRIERO E.Timedependent routing problems:a review.Computers&Operations Research,2015,64(1):189–197.

[9]MALANDRAKI C,DASKIN M S.Time dependent vehicle routing problems:formulations,properties and heuristic algorithms.Transportation Science,1992,26(3):185–200.

[10]ICHOUA S,GENDREAU M,POTVIN J Y.Vehicle dispatching with time-dependent travel times.European Journal of Operational Research,2003,144(2):379–396.

[11]AHN B H,SHIN J Y.Vehicle-routeing with time windows and time-varying congestion.Journal of the Operational Research Society,1991,42(5):393–400.

[12]FLEISCHMANN B,GIETZ M,GNUTZMANN S.Timevarying travel timesin vehicle routing.Transportation Science,2004,38(2):160–173.

[13]ALBIACH J,SANCHIS J,SOLER D.An asymmetric TSP with time windows and with time-dependent travel times and costs:an exact solution through a graph transformation.European Journal of Operational Research,2008,189(3):789–802.

[14]OSVALD A,STIRN L Z.A vehicle routing algorithm for the distribution of fresh vegetables and similar perishable food.Journal of Food Engineering,2008,85(2):285–295.

[15]KUOY.Using simulated annealing to minimize fuel consumption for the time-dependent vehicle routing problem.Computers&Industrial Engineering,2010,59(1):157–165.

[16]BALSEIRO S R,LOISEAU I,RAMONET J.An ant colony algorithm hybridized with insertion heuristics for the time dependent vehicle routing problem with time windows.Computers&Operations Research,2011,38(6):954–966.

[17]HASHIMOTO H,YAGIURA M,IBARAKI T.An iterated local search algorithm for the time-dependent vehicle routing problem with time windows.Discrete Optimization,2008,5(2):434–456.

[18]FIGLIOZZI M A.The time dependent vehicle routing problem with time windows:benchmark problems,an efficient solution algorithm,and solution characteristics.Transportation Research Part E:Logistics and Transportation Review,2012,48(3):616–636.

[19]KOK A L,HANS E,SCHUTTEN J.Vehicle routing under time-dependent travel times:the impact of congestion avoidance.Computers&Operations Research,2012,39(5):910–918.

[20]MINOCHA B,TRIPATHI S.Two phase algorithm for solving VRPTW problem.International Journal of Artificial Intelligence and Expert Systems,2013,4(1):1–16.

[21]BEHROUZ A N,ALIREZA A N.A constructive heuristic for time-dependent multi-depot vehicle routing problem with time-windows and heterogeneous fleet.Journal of King Saud University-Engineering Sciences,2014,29(1):29–34.

[22]SETAK M,HABIBI M,KARIMI H,et al.A time-dependent vehicle routing problem in multigraph with FIFO property.Journal of Manufacturing Systems,2015,35:37–45.

[23]LECLUYSE C,VAN WOENSEL T,PEREMANS H.Vehicle routing with stochastic time-dependent travel times.4OR,2009,7(4):363–377.

[24]DUAN Z,SUN S,SUN S,et al.Stochastic time-dependent vehicle routing problem:mathematical models and ant colony algorithm.Advances in Mechanical Engineering,2015,7(11):1–16.

[25]WEN L,EGLESE R.Minimum cost VRP with timedependent speed data and congestion charge.Computers&Operations Research,2015,56:41–50.

[26]DORIGO M,MANIEZZO V,COLORNI A.Ant system:optimization by a colony of cooperating agents.IEEE Trans.on Systems,Man,and Cybernetics,Part B:Cybernetics,1996,26(1):29–41.

[27]GAMBARDELLA L M,DORIGO M.Ant-Q:a reinforcement learning approach to the traveling salesman problem.Proc.of the 12th International Conference on Machine Learning,1995:252–260.

[28]DORIGO M,GAMBARDELLA L M.Ant colony system:a cooperative learning approach to the traveling salesman problem.IEEE Trans.on Evolutionary Computation,1997,1(1):53–66.

[29]STÜTZLE T,HOOS H H.MAX-MIN ant system.Future Generation Computer Systems,2000,16(8):889–914.

[30]BULLNHEIMER B,HARTL R F,STRAUSS C.A new rank based version of the ant system-a computational study.Central European Journal of Operations Research,1997,7(1):25–38.

[31]BLUM C,DORIGO M.The hyper-cube framework for ant colony optimization.IEEE Trans.on Systems,Man,and Cybernetics,Part B:Cybernetics,2004,34(2):1161–1172.

[32]GAMBARDELLA L M,TAILLARD E,AGAZZI G,et al.MACS-VRPTW:a multiple ant colony system for vehicle routing problems with time windows.Proc.of the New Ideas in Optimization,1999:63–76.

[33]DONATI A V,GAMBARDELLAL M,CASAGRANDE N,et al.Time dependent vehicle routing problem with an ant colony system.IDSIA Report,2003,185(3):1174–1191.

[34]DONATI A V,MONTEMANNI R,CASAGRANDE N,et al.Time dependent vehicle routing problem with a multi ant colony system.European Journal of Operational Research,2008,185(3):1174–1191.

[35]DONATI A V,MONTEMANNIR,GAMBARDELLA L M,et al.Integration of a robust shortest path algorithm with a time dependent vehicle routing model and applications.Computational Intelligence for Measurement Systems&Applications,2003,48(5):26–31.

[36]KUMAR S N,PANNEERSELVAM R.A comparative study of proposed genetic algorithm-based solution with other algorithms for time-dependent vehicle routing problem with time windows for e-commerce supply chain.Journal of Service Science&Management,2015,8(6):844–859.

[37]HUART V,PERRON S,CAPOROSSI G,et al.A Heuristic for the time-dependent vehicle routing problem with time windows.Cham,Switzerland:Springer International Publishing,2016.

[38]SOLOMON M M.Algorithms for the vehicle routing and scheduling problems with time window constraints.Operations Research,1987,35(2):254–265.

[39]POTVIN J Y,ROUSSEAU J M.A parallel route building algorithm for the vehicle routing and scheduling problem with time windows.European Journal of Operational Research,1993,66(3):331–340.

[40]ATKINSON J B.A greedy look-ahead heuristic for combinatorial optimization:an application to vehicle scheduling with time windows.Journal of the Operational Research Society,1994,45(6):673–684.

[41]IOANNOU G,KRITIKOS M,PRASTACOS G.A greedy look-ahead heuristic for the vehicle routing problem with time windows.Journal of the Operational Research Society,2001,52(5):523–537.

[42]LAU H C,SIM M,TEO K M.Vehicle routing problem with time windows and a limited number of vehicles.European Journal of Operational Research,2003,148(3):559–569.

[43]LIM A,ZHANG X.A two-stage heuristic with ejection pools and generalized ejection chains for the vehicle routing problem with time windows.Informs Journal on Computing,2007,19(3):443–457.

[44]LIN S.Computer solutions of the traveling salesman problem.Bell System Technical Journal,1965,44(10):2245–2269.

[45]FOSIN J,CARIC T,IVANJKO E.Vehicle routing optimization using multiple local search improvements.Automatika,2014,55(2):124–132.

[46]ILHAN O R.Traveling salesman-type combinatorial problems and their relation to the logistics of regional blood banking.Evanston,United States:Northwestern University,1976.

[47]OSMAN I H.Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem.Annals of Operations Research,1993,41(4):421–451.

[48]TAILLARD E,BADEAU P,GENDREAU M,et al.A tabu search heuristic for the vehicle routing problem with soft time windows.Transportation Science,1997,31(2):170–186.

[49]DOHN A,RASMUSSEN M S,LARSEN J.The vehicle routing problem with time windows and temporal dependencies.Networks,2011,58(4):273–289.

[50]GENDREAU M,HERTZ A,LAPORTE G.New insertion and postoptimization procedures for the traveling salesman problem.Operations Research,1991,40(6):1086–1094.

[51]DESROCHERS M,LENSTRA J K,SAVELSBERGH M W P,et al.Vehicle routing with time windows:optimization and approximation.Proc.of the Vehicle Routing:Methods&Studies,1987:65–84.

[52]GAMBARDELLA L M,TAILLARD E,AGAZZI G.MACSVRPTW:a multiple ant colony system for vehicle routing problems with time windows,IDSIA-06-99.Lugano,Switzerland:Technical Report Idsia,1999.

[53]STÜTZLE T,HOOS H H.Improving the ant system:a detailed report on the MAX-MIN ant system.FG Intellektik,FB Informatik,TU Darmstadt,1996,12(1):1–22.

[54]HOMBERGER J,GEHRING H.Two evolutionary metaheuristics for the vehicle routing problem with time windows.Information,1999,37(3):297–318.

[55]STÜTZLE T.Local search algorithms for combinatorial problems:analysis,improvements and new applications.Donauworth:STM Publishing House,1999.

[56]STUTZLE T,HOOS H.Improvements on the ant-system:introducing the MAX-MIN ant system.Proc.of the Artificial Neural Nets and Genetic Algorithms,1998:245–249.

[57]BENT R,VAN HENTENRYCK P.A two-stage hybrid local search for the vehicle routing problem with time windows.Transportation Science,2004,38(4):515–530.

[58]LIM A,ZHANG X.A two-stage heuristic with ejection pools and generalized ejection chains for the vehicle routing problem with time windows.INFORMS Journal on Computing,2007,19(3):443–457.