Time-dependent vehicle routing problem with path flexibility

Time-dependent vehicle routing problem with path flexibility

Citation for published version (APA):

Huang, Y., Zhao, L., van Woensel, T., & Gross, J-P. (2017). Time-dependent vehicle routing problem with path

flexibility. Transportation Research. Part B: Methodological, 95, 169-195.

https://doi.org/10.1016/j.trb.2016.10.013

Document license:

TAVERNE

DOI:

10.1016/j.trb.2016.10.013

Document status and date:

Published: 01/01/2017

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be

important differences between the submitted version and the official published version of record. People

interested in the research are advised to contact the author for the final version of the publication, or visit the

DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page

numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

ContentslistsavailableatScienceDirect

Transportation

Research

Part

B

journalhomepage:www.elsevier.com/locate/trb

Time-dependent

vehicle

routing

problem

with

path

flexibility

Yixiao Huang

a

_{, Lei Zhao}

a,∗

_{, Tom Van Woensel}

b

_{, Jean-Philippe Gross}

a a Department of Industrial Engineering, Tsinghua University, Beijing, 10 0 084, China

b School of Industrial Engineering, Eindhoven University of Technology, Eindhoven, Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 7 July 2016 Revised 23 October 2016 Accepted 24 October 2016 Available online 11 November 2016

Keywords:

Time-dependent vehicle routing problem Path flexibility

Geographical graph Stochastic travel time City logistics

a

b

s

t

r

a

c

t

Conventionally, vehicle routing problems are defined on a network in which the customer locations and arcs are given. Typically, these arcs somehow represent the distances or ex- pected travel time derived from the underlying road network. When executed, the quality of the solutions obtained from the vehicle routing problem depends largely on the quality of the road network representation. This paper explicitly considers path selection in the road network as an integrated decision in the time-dependent vehicle routing problem, denoted as path flexibility (PF). This means that any arc between two customer nodes has multiple corresponding paths in the road network (geographical graph). Hence, the de- cisions to make are involving not only the routing decision but also the path selection decision depending upon the departure time at the customers and the congestion levels in the relevant road network. The corresponding routing problem is a time-dependent vehicle routing problem with path flexibility (TDVRP–PF). We formulate the TDVRP–PF models un- der deterministic and stochastic traffic conditions. We derive important insights, relation- ships, and solution structures. Based on a representative testbed of instances (inspired on the road network of Beijing), significant savings are obtained in terms of cost and fuel con- sumption, by explicitly considering path flexibility. Having both path flexibility and time- dependent travel time seems to be a good representation of a wide range of stochasticity and dynamics in the travel time, and path flexibility serves as a natural recourse under stochastic conditions. Exploiting this observation, we employ a Route-Path approximation method generating near-optimal solutions for the TDVRP–PF under stochastic traffic condi- tions.

© 2016 Elsevier Ltd. All rights reserved.

1. Introduction

Urbanizationhasbecomeaworldwidephenomenonduringthepastdecades.Around54%oftheworld’spopulationlive incitiesin2014andthisisexpectedtobe66%by2050.Althoughinthemosturbanizedregion,i.e.,NorthernAmerica,82% ofthepopulationliveinurbanarea,only40%and48%ofthepopulationliveinurbanareasinAfricaandAsia,respectively (UnitedNations,2015).However,theratesofurbanizationinAfricaandAsiaarefasterthanthoseofother regions,which resultsinmoreandmoreemergingmega-citiesinthesetworegions,suchasChinaandIndia.

Handinhandwithurbanization(SavelsberghandVanWoensel,2016),trafficcongestionbecomesoneofthemajor chal-lengesnot only forcommutersbutalso forlogisticscompanies. Forexample, Beijingisranked 15 amongover 200 large

∗ _{Corresponding author.}

E-mail address: lzhao@tsinghua.edu.cn (L. Zhao). http://dx.doi.org/10.1016/j.trb.2016.10.013

(a) Path off the ring road is faster at 07:30 am on a workday

(a) Path off the ring road is fasteff r at 07:30 am on a workday

(b) Path on the ring road is faster at 10:00 am on a workday

Fig. 1. Two popular morning commute paths from a residential area to Tsinghua University in Beijing, China.

cities(population> 800K)intheworld,basedontheoverallcongestionlevelin2014(TomTom,2015).Congestion signifi-cantlyincreasesthefuelconsumptionofvehicles,andthereforeincreasesthecostoflogisticscompanies.In2011,thecostof delaysonroadsinthesixlargestcapitalcitiesinAustraliaisestimatedtobearound $13.7billion(InfrastructureAustralia, 2015).In China,fuel costaccountsforover46% ofthe totaloperationalcost foralogistics company(China Federationof Logistics&Purchasing,2014).

Firstmilepickupsandlast miledeliveries,oftenmodeledasavehiclerouting problem,aremoreandmoresituatedin thesecongestedurbanizedregions.Vehicleroutingproblemsaremostlydefinedonacustomer-basedgraph.Arcsbetween thesecustomers areassumedto adequately representthedistances orexpectedtravel time.As such, thepath usedby a vehicle(intheroadnetwork)betweenapairofcustomersisaggregatedintothearccost.Whenexecuted,thequalityofthe solutionsobtaineddependslargelyonthequalityofthisroadnetworkrepresentation.Typically,intheliterature (seee.g.,

Fleischmannetal., 2004;Gendreauetal.,2015;Ichouaetal., 2003;Jabalietal.,2012),thistrafficdependencyismodeled implicitly,resultingintime-dependentand/orstochastictraveltimevariantsofthevehicleroutingproblem.

This papertakesan alternative approach andexplicitlyconsiders path selection inthe roadnetwork asan integrated decision in the vehicle routing problem that minimizes the total cost, including fuel cost and vehicle depreciation cost. Thismeans thatanyarcbetweentwo customernodeshas multiplecorrespondingpaths inthe roadnetwork (geographi-cal graph). Hence, thedecisions tomake are involvingnot onlythe routing decisionbutalso thepath selection decision dependingupon thedeparture time atthecustomersandthe congestionlevelsinthe relevantroadnetwork. Wedenote thepathselectiondecisionaspathflexibility(PF).The correspondingroutingproblemisatime-dependentvehiclerouting problemwith path flexibility (TDVRP–PF). We formulate the TDVRP–PF models underdeterministic and stochastic traffic conditions.

Consider the examplegivenin Fig.1, thereare two typical path options,one via the 4thRingRoad (the lower path) andtheotheroff theringroad(theupperpath).Betweenthesetwopaths,thefasterpathdependsontheactualdeparture time andon the (spread ofthe) traffic congestionin theroad network. Notethat the congestionin theroad networkis alsoregion-dependent.ConsidertheexampleinFig.2,wherethecolorsoftheroadsrepresentthetrafficconditions(green forfreeflow, yellowformediumcongestion, andred forheavy congestion).Observethat thetrafficconditionsofthering roads changesignificantlybetweenpeak andoff-peak hours,whilethose oftheroads withinthe2nd ringremain almost unchanged. Jointly considering routing andpath selection decisionsmakes it possible toexplicitly considerthe temporal andspatialdifferencesofcongestionintheroadnetwork.

2nd Ring Road 3rd Ring Road 4th Ring Road (a) 08:00 am 2nd Ring Road 3rd Ring Road 4th Ring Road (b) 10:00 am 2nd Ring Road 3rd Ring Road 4th Ring Road (c) 12:00 pm

Fig. 2. Estimated traffic conditions of urban area in Beijing on Mondays (Source: Google Map).

There have been few studies explicitly considering path flexibility in the vehicle routing problem. Garaix et al. (2010) workedon a multi-attribute vehicle routing problem andproposed a multigraph to representthe network. Setak etal.(2015) solved thetime-minimizing time-dependentvehicle routing problemin a multigraphwithFirst-In-First-Out (FIFO)propertyviatabusearch.Ehmkeetal.(2016b)andQianandEglese(2016)bothdealtwiththe time-dependent vehi-cleroutingproblemtominimizeemissions,whichalignwithourworkwell.Ourpaperfitsintothisstreamofliterature.We focusontheanalysisofoptimizationmodelsandtheiroptimalsolutions,whileEhmkeetal.(2016b)andQianandEglese (2016)proposed efficient andeffectivetabusearch methods tosolve largerinstances. We flexiblyandexplicitlykeep the precomputedpaths ascandidatepaths andthepath selection isembeddedasdecisionvariables inthe TDVRP–PF,under deterministicandstochastictrafficconditions.

Ourpapermakesanumberofcontributionstotheliterature.

1.Weexplicitlyconsiderpathflexibilityinthetime-dependentvehicleroutingproblem(TDVRP–PF).Wemodelthe TDVRP-PF under deterministic and stochastic traffic conditions, the latter ofwhich leadsto a two-stagestochastic program. Ratherthandevelopingadvancedsolutionalgorithmsforourproblemon-hand,weuseanindustrialsolver(CPLEX)and derive importantinsights,relationships, andstructures. Weobservethat,whentheTDVRPandTDVRP–PFsolutionsare structurally different, the cost saving of path flexibility is usually more significant. Thisis important in practice. For instance,home appliance deliverycompanies loadtheappliances (normally bigandheavy)accordingtotheir delivery sequence.Itisdifficulttochangethesequenceenroute.Theseresultsserveasagoodbasisforfurtheralgorithmicwork, whichisleftforfutureresearch.

2.We demonstrate the benefits of path flexibility based on a representative testbed of instances inspired on the road networkofBeijing(i.e.,similartoQianandEglese,2016,whousedtheLondonroadnetworkfortheirinstances).Forthis richsetofinstances,significantsavingsareobtainedintermsofcostandfuelconsumptionbyexplicitlyconsideringpath flexibility into therouting andscheduling. Itis alsodemonstrated that theseimprovementsare more significant than timeflexibility(i.e.,flexibledeparturetimeatthedepotand/orpost-servicewaitingtimeatcustomersites).Additionally, theseobservationsarerobustbasedonourextensivesensitivityanalysis,withregardtothecustomernetworktopology andthecorrespondingchoicesonthenumberofalternativepathsinthegeographicalgraph.

3.We observethattheroutesgeneratedwiththeTDVRP–PFunderdeterministictraffic conditionsareverygood approx-imations for theTDVRP–PF understochastic traffic conditions. Havingboth path flexibility andtime-dependenttravel time seems to be a good representation of a wide range of stochasticity and dynamics inthe travel time, and path flexibility servesasanaturalrecourseunderstochastic conditions.Exploitingthisobservation,weemploya Route-Path approximationmethodgeneratingnear-optimalsolutionsfortheTDVRP–PFunderstochastictrafficconditions.

Intheremainderofthispaper,wefirstreviewtherelevantliteratureinSection2,andthenformallydescribetheTDVRP– PFin Section 3. In Section 4, we provide some preliminarieson fuel consumption and time-dependentmodels, andwe formulatetheTDVRP–PFasmixed-integer programsinSection 5,includingbothdeterministic andstochastic versions.We describe thedesign of experiments anddiscussthe numerical resultsin Section 6. Finally,we summarize ourwork and providethesuggestionsforfutureresearchinSection7.

2. Literaturereview

Our work is closely related to the green vehicle routing problem (green VRP), which considers carbon emission or fuel consumption ofvehicles, andthe time-dependentvehicle routing problem(TDVRP),which reflects the varying traf-ficconditions.Interestedreaderscan referto Linetal.(2014)andDemir etal.(2014b)forreviewsofthe greenVRPand

Gendreau et al.(2015) fora recentreview ofthe TDVRP.As shortlyoutlined in the introduction,our paperextends and buildsupontheliteratureinanumberofdirections.

Theconventionaltime-minimizingTDVRPprovidesthefoundationofourresearchonthetime-dependencymodelling.It isstandardtousetime-dependenttimematricestostorethefastestpaths ofdifferentdeparturetime.Duetothe First-In-First-Out(FIFO)network,itisoptimaltoalwayschoosethefastestpathcorrespondingtoacertaindeparturetime.However, thismethodcannot beextendedto acost-minimizingTDVRP,wherethe challengeis,asbrieflymentioned inSavelsbergh andVanWoensel(2016),thecostsaredependentonspeed,distance,andvehicleload.Duetothecomplexstructure,itmay notbegloballyoptimalfortheroutingsolutionstochoosethelocallycost-minimizingpathbetweeneachpairofcustomers inaroute.Alternatively,wekeepmultiplepathsasthepoolofcandidatepathsandproposeacorrespondingmathematical modelallowingthesimultaneousdecisionsofbothroutingandpathselection.Intheirpaper,Koç etal.(2016)analyzedthe impactofspatialdifferencesofcongestion(speedzones).Wefocusonpathflexibilityconsideringbothspatialandtemporal differences.Ehmkeetal.(2016b)andQianandEglese(2016)proposedmetaheuristicmethods,whilewefocusonobtaining insightsandtheanalysisofoptimalsolutions.

Althoughourworkissimilarto thePollution-RoutingProblem(PRP, Bekta¸s andLaporte, 2011),thereareseveralmajor differencesintheproblemsettings.First,wedonotconsiderthedrivers’salary.Second,wedonotcontrolthevehiclespeed aswefocusonthedecisionsofroutingandpathselection,whichisusuallydecidedbytheplanningdepartment,ratherthan thedetaileddrivingbehavior.Further,astheroadsareusuallycongested,weassumethatthevehiclesjustfollowtraffic.Last, weconsidertheroutingprobleminatime-dependenttrafficconditionandexplicitlyincorporatedecisionsofpathselection intheroutingproblem.

2.1. Time-dependentrouting

The TDVRPhas beenintroduced tocapture congestionina traffic network.Malandraki andDaskin(1992)proposed a mixed integer formulation of the time-dependent vehicle routing problem where the travel time on each arc is a step function of the departure time. Malandraki andDial (1996) applied restricted dynamic programmingto solve the time-dependent travelingsalesmanproblem. HillandBenton(1992) introduceda modelfortime-dependenttravel speedsand methodsforparameterestimation,whichcanbeappliedforthetime-dependentvehicleroutingproblems.However,allthe studiesdisrespecttheFIFOpropertyinthetime-dependentnetwork,thatis,vehiclesleavinglatermayarriveearlieronthe samearc.

Twomethods (Fleischmannetal., 2004 andIchouaetal., 2003) appeared tomodel thetime-dependentnetwork that respectstheFIFOproperty.Fleischmannetal.(2004)modeledthetraveltimeviathesmoothedtraveltimefunctionbased onthestepfunctionoftraveltime.Ichouaetal.(2003)presentedamodelwheretravelspeedsarestepfunctionsandtravel timebecomepiecewiselinearfunctions.Thesemethodsarecommonlyappliedinresearchafterwards.

Theunderlyingtime-dependentshortestpathproblem(TDSPP)isessential totheTDVRP.Fleischmannetal.(2004) an-ticipatedthepotentialvalue ofconsideringpath changingintheTDVRPandsuggestedit astheresearch direction. Eglese etal.(2006)constructedatimetableoftime-dependentshortestpathsbasedonthereal-worldroadnetworkandthe cor-respondingfloatingvehicle datainEngland.Koketal.(2012)solved theTDVRPwithtime-dependentshortestpaths viaa restricteddynamicprogrammingheuristicandshowedthesignificant reductionin traveltime tousethe time-dependent shortestpathsintheTDVRP.Similarly,Ehmkeetal.(2012)alsoshowedthesuperiorityoftime-dependentvehiclerouting overstaticvehicleroutinginacity logisticscasestudy.Theyalsoprovideddetailedanalysesofspatio-temporal structures ofthetime-dependentroutingsolutions.

Madenetal.(2010)andFigliozzi(2011)derived thedailyspeedlevelsbasedonthespeedobservedin15-minutetime binsaveragedover3monthperiodofrealtrafficdata,inordertoprovideagoodapproximation.Alternatively,studiessuch asKuo(2010),Jabalietal.(2012),andFranceschettietal.(2013)usefewer(from2to5)speedlevelstodepicttypicaltraffic conditionsin urban areas.Followingthistrend,we studythesituation of3speed levels inthemorning shift,wherethe secondspeedlevelisofthepeakhour.AlthoughEhmkeetal.(2016b)andQianandEglese(2016)usefinerspeedlevelsto depictthe24hourtrafficconditions,thedeliverytimehorizonintheirinstancesusuallycontainonlyonepeakhourperiod.

2.2. Greenrouting

There have been increasing and various studies concentrating on the green vehicle routing problemin recent years.

ErdoˇganandMiller-Hooks(2012)consideredafleetofalternativefuel-poweredvehicleswiththerestrictionoflimited driv-ingrangeandlimitedrefuelinginfrastructure,andthusthenecessitytoincludethedecisionofvisitingfuelstationinthe vehicleroutingproblem.KopferandKopfer (2013)andKwonetal. (2013)consideredtheeffectofaheterogeneousfleetof vehiclesinthegreenvehicleroutingproblem.ResearchonthecasestudyofrealworldpracticeincludesUbedaetal.(2011),

Figliozzi(2011),etc.

Pollution-RoutingProblemisanimportantvariantofthegreenVRP.Bekta¸sandLaporte(2011)firstcoinedthe Pollution-RoutingProblem(PRP),whichalsofirstusedtheComprehensiveModalEmissionModel(CMEM)tocalculatefuel consump-tion.Demir etal.(2014a) extendedthePRPtoabi-objective PRPwheretheobjectivefunctionsaretominimize bothfuel consumption anddriving time.Koç et al.(2014) considered thePRP withaheterogeneous fleet ofvehicles, todetermine thenumber ofvehiclesofeach vehicletype andtherouting ofthesevehicles, inordertominimize therouting cost and

thefixed cost.Tajik etal. (2014) dealtwiththe PRP withpickupsand deliveries underuncertainty andused the robust optimizationapproachtosolvetheproblem.Koç etal.(2016)consideredemissionminimizationwithinthelocation-routing problemwithaheterogeneousfleet,wherethecityisdividedintodifferentzones,eachwithaconstant(time-independent) speed.Asthevehicletravelsamongthezoneswithdifferentload,itmaychoosedifferentpathsbetweenzonestominimize emissions.

Ontheperspectiveofalgorithmicdesign,Demiretal.(2012)proposedanadaptivelargeneighborhoodsearchalgorithm tosolvethePRP.Krameretal.(2015b)developedahybriditeratedlocalsearchapproachwithasetpartitioningprocedure andaspeedoptimizationalgorithm.Krameretal.(2015a)furtherproposed anexactmethodonspeedanddeparturetime optimization,which improves the performance of the iterated local search method developed by Kramer et al. (2015b).

Dabiaetal.(2016)presentedanexactmethodbasedonabranch-and-pricealgorithm,includingnewdominancecriteriato discardunpromisinglabels.

2.3.Greentime-dependentrouting

As congestionsignificantly affects fuel consumption,there are a stream ofstudies taking fuelconsumption or carbon emissionintoconsiderationinthetime-dependentvehicleroutingproblems,whichwerefertoastheGreenTDVRP.

Madenetal.(2010)proposedaLANTIMEheuristicalgorithmtominimizethetotaltraveltime.Theycalculatedthe emis-sionsofthetime-minimizingsolutionsandshoweda7% emissionreduction inthecasestudy.Kuo (2010)introduced the TDVRPtominimize fuelconsumptionandsolved theproblemusingsimulatedannealing.Theresultshowedaclear trade-off betweentime-minimizingandfuel-minimizingsolutions.Figliozzi(2011)minimizedfuelconsumption intheTDVRPin urbanfreightdistributionbasedontherealnetwork ofPortland.Jabalietal.(2012) studiedatime-dependentgreen vehi-cleroutingproblem,whichusestheMethodologyforCalculatingTransportationEmissionsandEnergyConsumption(MEET,

Hickmanetal.,1999)emissionmodeltoevaluatetheCO2emissions.However,thesestudiesignoretheloadwhen calculat-ingtheemissions.

Franceschettietal.(2013)introducedthetime-dependentpollution-routingproblem,usingamixedintegerprogramming formulation.They obtained the analytical results of speed and departure time optimization based on a two-step speed functionwithcongestionfirstandfreeflowsecond,wheretravelspeedsinthefreeflowaredecisions.However,theirresults aredifficult togeneralize tospeed functionswithmorethantwo steps. Soysaletal.(2015) considered atime-dependent two-echelonvehiclerouting problem,whichminimizesthetotalcostincludingfuelcostforthetwoechelons,drivers’cost forthe two echelons, andhandling cost inthe satellites(second-echelon depots). The numerical experimentalso uses a two-step speed function with congestionfirst andfree flow second. Bothstudies depict time-dependency asa two-step speedfunction oneach arcwithout consideringthe underlyingpath(s). Therefore,thecorresponding travel timefunction andfuelconsumptionfunctiononeacharcarepredeterminedbeforerouting.Whentheunderlyingpathsareconsidered,as inourwork,thetraveltimefunctionandfuelconsumptionfunctionassociatedwitheacharcdependontheselectedpath andvehicleload,wherethelatterisdeterminedbytheroutingdecisions.

Thetime-dependentemission-minimizingpathproblem,theunderlyingproblemofthegreenTDVRP,issomehow under-studied.UnliketheTDSPP,itisgenerallyvery difficultorimpossibletocalculatethetime-dependentemission-minimizing path,asitrequirestheinformationofbothdeparturetimeandvehicleloadoneacharc,whichdependontherouting de-cision.Itmightbepossibletoprecomputethepathsbasedoneverycombinationofdiscretizeddeparturetime andvehicle load, whichis by itself computationallydemanding and thusmakes the routing problemovercomplicated. Therefore, the existing studies(e.g., QianandEglese,2014 andEhmke etal., 2016a) useemission models thatneglect the vehicleload, suchasMEETandNationalAtmosphericEmissionInventory(NAEI,2013).

Ehmkeetal.(2016b) dealtwiththetime-dependentvehiclerouting problemtominimizeemissions. Theyinvestigated howtosimplifytheprecomputationofemission-minimizingpathsvia identifyingtheoptimalpaths thatare independent ofthevehicleloadanduseda tabusearch methodto solvetheroutingproblem.Theirtestinstancesare basedonalarge realroadnetworkandextensiveexperimentsaretestedforanalysis.QianandEglese(2016)proposedacolumngeneration basedtabusearchmethodtosolvethetime-dependentemission-minimizingrouting problem,whichimplicitlyconsidered the path selection in the pricing subproblems. Besides, Qian andEglese (2016) use the NAEI emission model, which is independentofthevehicleload.

3. Problemdescription

Fortheconvenienceofthereader,weprovideTable1listingthenotationusedinthispaper.

A logistics company operates a homogeneous fleet of vehicles to deliver goods to customers that are geographically scatteredinanurbanarea.Thedeliveriesaresubjecttovehiclecapacityanddeliverytimehorizonrestrictions.Theobjective istominimize theoperationalcost,includingfuelconsumptioncostandvehicledepreciationcost.Inthelogisticspractice inChina,alogisticscompany usuallypaysdrivers witha fixedbasesalaryplusthe bonusbasedonthe deliveryquantity. Drivers’costisthereforeindependentoftherouting decisionandisaconstant.Thus,we donotconsiderthedrivers’cost intheobjectivefunction.

Table 1 Notation.

Notation Description First appeared section

Gc _{Customer graph 3}

Nc _{Set of customer nodes 3}

Ac _{Set of arcs connecting pairs of customer nodes 3}

Gg _{Geographical graph 3}

Ng _{Set of intersections 3}

Ag _{Set of road segments 3}

[ L, U ] Delivery time horizon 3

Pij Set of paths connecting arc ( i, j ) ∈ A c 3

P p

i j A path in the set P ij 3

qi Demand of customer node i ∈ N c ࢨ{0} 3

si Service time of customer node i ∈ N c ࢨ{0} 3

K Number of vehicles 3

Q Capacity of a vehicle 3

R Number of time periods with constant speed 3

ψe , ψs , ψw Coefficients in CMEM model to calculate fuel consumption 4.1

μ Vehicle curb-weight 4.1

F ( t, v, f ) Fuel consumption function of constant speed v , duration t , and vehicle load f 4.1 vpr

i j Speed in time period r on path P p

i j 4.2

w Departure time from a node 4.2

τp

i j(w) Piecewise linear travel time function of departure time w from node i on path P p i j 4.2 Mp

i j Intervals of the piecewise linear travel time function τ p

i j(w) on path P p

i j 4.2

[ b pm

i j , b p,mi j ] +1 m th interval of travel time function τ p

i j(w) on path P p

i j 4.2

tpmr

i j Travel time at speed v pr

i j if the vehicles departs from node i at time b pm

i j on path P p i j 4.2 τpr

i j(w) Piecewise linear travel time function of departure time w at speed v pr

i j on path P p i j 4.2 θpmr

i j Slope of piecewise linear function τ pr i j(w) within interval [ b pm i j , b p,m+1 i j ] on path P p i j 4.2 ηpmr

i j Intercept of piecewise linear function τ pr i j(w) within interval [ b pm i j , b p,mi j ] on path P +1 p i j 4.2 di jp Distance of path P p i j 4.2

f_{i j}p Vehicle load carried on path P p

i j 4.2

F_{i j}p(w) Fuel consumption of departure time w from node i on path P p

i j 4.2

S Number of discretized time points to calculate the time-dependent shortest path 4.3 xij Equal to 1 if arc ( i, j ) is on the optimal route 5.1 xp

i j Equal to 1 if the vehicle travels on path P p

i j 5.1

xpm

i j Equal to 1 if the departure time at node i toward node j is in [ b pm

i j , b p,mi j ] +1 5.1

fij Load carried on arc ( i, j ) 5.1

ϑpm

i j Earliest possible departure time at node i on path P p

i j 5.1

ζpm

i j Time flexibility at node i on path P p

i j 5.1

wpmi j Actual departure time at node i (in [ b pm

i j , b p,mi j ] ) on path P +1 p

i j 5.1

cf _{Unit cost of fuel consumption 5.1}

cd _{Unit cost of vehicle depreciation 5.1}

Thetrafficconditionoftheurbanroadnetworkvariesbothtemporallyandspatially.Therefore,whentravellingbetween twolocationsinthecityatdifferenttimeoftheday,avehiclemaychooseadifferentpath,dependingonthetraffic condi-tion.

The customer locationsandthe roadnetwork canbe represented bytwo layers ofa graph: customer graphand geo-graphicalgraph(Fig.3).Followingtraditionalsettingsofvehicleroutingproblems,inthecustomer graphGc₌

₍

_Nc_,Ac

₎

_,the

nodesetNc_{representsthesetofcustomerlocationsandthedepotlocation,whilethedirectedarcsetA}c_{connectsthenode}

pairs. Specifically,node 0 ∈Nc _{represents thedepot. Each customer node i ∈N}c_{ࢨ{0}is associatedwitha demand q} i and

servicetimes_i.ThehomogenousfleethasKvehicleseachwithacapacityQ.Vehiclesareroutedtosatisfyallthecustomer demandswithinthedeliverytimehorizon[L,U].

The geographicalgraphGg₌

₍

_Ng_,Ag

₎

_{representstheactual roadnetworkunderlying thecustomergraph G}c_.Ng_denotes

thesetofintersectionsandAg _{denotesthesetofroadsegments intheroadnetwork.Eacharc(i,j)∈A}c _{inthecustomer}

graphcancorrespondtoa setofmultiplepaths inthegeographicalgraph,denotedasPi j=

{

Pi j1,Pi j2,...,

}

.ApathP p i j, p=

1_,...,

|

P_{i j}

|

_, connectsnodes i and j in Nc_{,and P}p

i j consists of a finitesequence of road segments in A

g_{, denoted asP}p i j=

{

a₁p,a₂p,...

}

.

Fig. 4shows the speed patternsof the roadnetwork in Beijingfrom October14th to 18th(Monday to Friday), 2013. Whilethespeedmayvaryinregionordistrict,thespeedpatternsoftheentireroadnetworkremainsimilar.Therefore,we partitionthetimeofthedayintoRperiods,andassumethat,withineachtime periodr=1,...,R,theaveragespeedofa roadsegmentaremainsstable,denotedas v r

a.

Besidespathflexibility,wealsostudytheeffectoftimeflexibility,whichrefers totheflexibledeparture timeafterLat thedepot(i=0)and/orthepost-service waitingtime atacustomer nodei ∈Nc_{ࢨ{0},toavoidcongestionontheroad.The}

Fig. 3. Customer graph and geographical graph. 0 10 20 30 40 50 60 0:00 4:00 8:00 12:00 16:00 20:00 0:00 A v erage speed (km /h) Time

Entire road network Inside 2nd ring 2nd ring to 3rd ring 3rd ring to 4th ring 4th ring to 5th ring Dongcheng District Xicheng District Haidian District Chaoyang District Fengtai District Shijingshan District

Fig. 4. Speed patterns of the road network in Beijing (Source: Beijing Municipal Commission of Transport).

restrictionoftheloading/unloadingzone.Notealsothattime flexibilityisdifferentfromthewaitingtimebeforethe start ofserviceatacustomernodeinthecaseofearlyarrivalinthevehicleroutingproblemwithtimewindows(VRPTW).

Furthermore,thetrafficconditionoftheroadnetworkisalsostochastic.Thatis,whenthelogisticscompanymakesthe routingdecisionsatthebeginningofthetimehorizon,theactual roadcondition(orcongestionlevel)ofaparticularroad segmentataparticulartimeperiodremainsuncertain.

4. Preliminaries

Inthissection,wedescribethefundamentalsofthefuelconsumptionmodel,thestructureofthetime-dependenttravel timemodel,andamethodtoobtainthetime-dependentshortestpaths.Thissectionbasicallypreparesforourmodel for-mulationinSection5.

0 5 10 15 20 25 30 35 40 45 50 10 20 30 40 50 60 70 80

Fuel consumption (liter/100km)

Speed (km/h)

Total fuel consumption Engine module Speed module Weight module

Fig. 5. Fuel consumption function of speed (with 20 0 0 kg load).

0 1 2 3 0 2 4 6 T ravel speed Departure time

(a) Step speed function

0 1 2 3 0 1 2 3 4 5 6 T ravel tim e Departure time

(b) Piecewise linear travel time function

Fig. 6. An example of converting step speed function to piecewise linear travel time function on a road segment, with 3 units of distance.

4.1. Fuelconsumptioncalculation

Same asrelevant papers in the literature (e.g.,Bekta¸s and Laporte, 2011; Franceschettiet al., 2013), we alsouse the Comprehensive Modal EmissionModel(CMEM, Barthet al., 2005 andBarthand Boriboonsomsin, 2009) to calculatefuel consumption.

Foravehicletravelingadistancedataconstantspeedvindurationt(t₌d_/

v

)andcarryingaloadf,thefuel consump-tionF(t,v,f)canbecalculatedas

F

(

t,

v

,f

)

=

ψ

et+

ψ

s

v

3t+

ψ

wd

(

μ

+f

)

, (1)

where

μ

isthe vehiclecurb-weight,and

ψ

e,

ψ

s,and

ψ

w are the coefficientsoftheengine, speed,and weightmodules,

respectively.Fig.5illustrates thefuelconsumptionfunctionofatypicallight-dutyvehicleinChina (parametersandvalues areinTableA.16).ThefunctionhasaU-shapewiththeminimumfuelconsumptionatthespeedof52.27km/h.

DetaileddescriptionoftheCMEMisinAppendixA.

4.2. Time-dependenttraveltimemodel

As discussed in Section 3,we assume that the delivery time horizon[L, U] consists ofR time periods andthe travel speed ofa particular road segment a is constant within a period r. Thus the travel speed of road segment a is a step function ofthetime oftheday. Thecorresponding travel timefunction ofroadsegment a isapiecewise linearfunction, whereweusethesamemethodasdescribedinIchouaetal.(2003)fortheconversionfromthestepspeedfunctiontothe piecewiselineartraveltime functiononaroadsegment,illustratedinFig.6.ForapathPp

i j withafinitesequenceofroad

segments,eachofwhichhasacontinuouspiecewiselineartraveltimefunction,weusethemethoddescribedinGhianiand Guerriero (2014)to generate

τ

_{i j}p

(

·

)

,thetravel time function onthe entirepathPp

i j,which isalsoa continuouspiecewise

linearfunction.

Consideringthepiecewiselineartraveltime function

τ

_{i j}p

(

·

)

ofpathPp

i j ofarc(i,j),itconsistsofM p

i j intervalsdefinedby

thatbp1_{i j} =Landbp,M

p i j+1

i j =U,asthedomain(setofinputs)ofthetraveltimefunction

τ

p

i j

(

·

)

isthedeliverytimehorizon[L,

U].

Lett_{i j}pmr denotethetraveltimethevehicletravelsatspeed v _{i j}pr alongthepathPp

i j,ifitdepartsfromnodeiattimeb pm i j .

Ifthevehicledepartsfromnodeiattimew∈[bpm_{i j} ,b_{i j}p,m+1],thetraveltimeittravelsatspeed v pr_{i j} is

τ

pr i j

(

w

)

= t_{i j}p,m+1,r− tpmr i j b_{i j}p,m+1− bpm i j

(

w− bpm i j

)

+t pmr i j =

θ

pmr i j w+

η

pmr i j , (2) where

θ

pmr i j = t_{i j}p,m+1,r− tpmr i j bp_{i j},m+1− bpm i j , and (3)

η

pmr i j =t pmr i j − b pm i j t_{i j}p,m+1,r− tpmr i j b_{i j}p,m+1− bpm i j . (4) Clearly,t_{i j}pmr=

τ

pr i j

(

b pm

i j

)

.ThetotaltraveltimeonpathP p i j isthen

τ

p i j

(

w

)

= R  r=1

τ

pr i j

(

w

)

= R  r=1

θ

pmr i j w+ R  r=1

η

pmr i j . (5)

Therefore,thefuel consumption ofa vehicle travelingonpath _{Pi j}p ofarc (i, j) anddeparting fromnode i at time w_∈

[bpm_{i j ,}bp,m_{i j} +1]can be calculated as F_{i j}p

(

w

)

=

ψ

e



 r∈R

θ

pmr i j w+  r∈R

η

pmr i j



+

ψ

s



 r∈R

(

v

pr i j

)

3

_θ

pmr i j w+  r∈R

(

v

pr i j

)

3

_η

pmr i j



+

ψ

wd_{i j}p

(

μ

+f_{i j}p

)

, (6)

whered_{i j}p denotesthedistanceofpath_{Pi j}pand f_{i j}pdenotestheloadcarriedbythevehicleonpath_{Pi j}pofarc(i,j).

4.3.Time-dependentshortestpathmethod

In a time-dependent network where the travel time function of each arc is piecewise linear, the goal of the time-dependent shortest path problem is to find the shortestpaths from i to j (measured by travel time),for all possible de-parturetime atnode i.This problemisNP-hard inits generalform(Dehneetal., 2012). However, fora particulargiven departuretimewatnodei,we canuseamodifiedDijkstra’salgorithm(Algorithm1)tofindthetime-dependentshortest pathfromnodei toj efficiently(BrodalandJacob,2004).Notethat,inStep2.2(a),we calculatethevalueTnbasedonthe

time-dependenttraveltimefunction

τ

_nn

(

Vn

)

.

Algorithm1 ModifiedDijkstra’salgorithmforthetime-dependentshortestpathonarc(i,j). Step1.Initialization.

1.1.Labeltheorigincustomernodeias[0,−],where[Vi,nipre]representsthevalueViand

thepredecessornoden_ipreofnodei. 1.2.LabelallothernodesinNg_as[∞,−].

1.3.Setcustomernodeias“scanned.” Step2.Labelupdate.

2.1.Denotethelastscannednoden . 2.2.

∀

(

n ,n

)

∈Ag_andn∈Ng_notscanned,

a)ComputeTn=Vn+

τ

nn

(

Vn

)

,where

τ

nn

(

Vn

)

representsthetraveltimeon

theroadsegment

(

n ,n

)

whenthevehicledepartsfromnoden attimeV_n. b)IfTn<Vn,relabelnodenwith[Tn,n ].

Step3.Scananode.

FindtheunscannednodewiththesmallestVnandscanthenode.

Step4.Terminationcheck.

Ifthecustomernode jisscanned,stop. Else,gotostep2.

In this paper, we solve the time-dependent shortest path problem by Algorithm 1 with S discretized time points

w=w1,...,wS asthe givendeparturetime atnode i, tocollecta setofdistinctive time-dependentshortestpaths

(time-minimizingpaths),Pi j=

{

Pi j1,Pi j2,...,

}

.Further,besidesthetime-minimizingpaths,wealsoincludethedistance-minimizing

pathonarc(i,j)inP_ij.Thesetofpathsserveascandidatepathstobeselectedwhenthevehicletravelsonarc(i,j).

5. Mathematicalmodels

Inthissection,weformulatethedeterministictime-dependentvehicleroutingproblemwithpathflexibility(TDVRP–PF) asamixedintegerprogramandthestochasticTDVRP–PFasatwo-stagestochasticmixedintegerprogram.

5.1. DeterministicTDVRP–PF

Underdeterministictrafficconditions,thevehicleroutingandpathselectiondecisionsaremadesimultaneouslyto mini-mizetheoperationalcost,thatis,thefuelconsumptionandvehicledepreciationcost.ThesetofcandidatepathsPijofeach

arc(i,j)∈Ac_{arepreprocessedasdiscussedinSection4.3.}

Wedefinethedecisionvariablesas xi j =



1, ifarc

(

i,j

)

isontheoptimalroute, 0, otherwise.

xp_{i j} =



1, ifthevehicletravelsonpathPp i j,

0, otherwise.

x_{i j}pm =



1, ifthedeparturetimeatnodeitowardnode jisin[bpm_{i j} ,b_{i j}p,m+1], 0, otherwise.

fi j =loadcarriedonarc

(

i,j

)

.

f_{i j}p=loadcarriedonpathPp i j.

ϑ

pm

i j =earliestpossibledeparturetimeatnodeionpathP p i j.

ζ

pm

i j =timeflexibilityatnodeionpathP p i j.

w_{i j}pm =

ϑ

pm i j +

ζ

pm

i j =actualdeparturetimeatnodei(in[b pm i j ,b p,m+1 i j ])onpathP p i j.

Weintroducexp_{i j} and f_{i j}ptospecifythepathselectiononarc(i,j),andintroducexpm_{i j ,}w_{i j}pm_,

ϑ

_{i j}pm_,and

ζ

_{i j}pm tospecifythe departuretimeatnodeiwithinatimeinterval[b_{i j}pm_,b_{i j}p,m+1]ofthepiecewiselineartraveltimefunctionalongpath_{Pi j}p.

Note that when we preprocess the set of candidate paths Pij for each arc (i, j) ∈ Ac, the information related to the

underlyingroadnetworkGg₌

₍

_Ng_,Ag

₎

_{isall embeddedinP}

ij.Therefore,intheTDVRP–PF,wefocusonthecustomer graph

Gc₌

₍

_Nc_,Ac

₎

_{andthe explicitcandidatesetof pathsP}

ij oneach arc(i,j) ∈Ac,obtainedfromthegeographical graph.For

model simplicity,we drop the superscript“c” and useG=

(

N,A

)

to representthe customer graphGc₌

₍

_Nc_,Ac

₎

_.Further,

foreacharc(i,j)∈A,thereare |Pij|candidatepaths, andforeach pathPi jp∈Pi j,thereare

|

Mi jp

|

intervalstorepresentthe

piecewiselineartraveltime function.Inthefollowing,inordertosimplifythenotation,wedropthesuper/subscriptsand usePforP_ijandMforM_{i j}p.

TheTDVRP–PFisformulatedasamixedintegerprogramasfollows.

Min

(

TDVRP−PF

)

cf



 (i, j)∈A  p∈P  m∈M

ψ

e

(

 r∈R

θ

pmr i j w pm i j +  r∈R

η

pmr i j x pm i j

)

(7) +  (i, j)∈A  p∈P  m∈M

ψ

s



 r∈R

(

v

pr i j

)

3

_θ

pmr i j w pm i j +  r∈R

(

v

pr i j

)

3

_η

pmr i j x pm i j



(8) +  (i, j)∈A  p∈P

ψ

wd_{i j}p

(

μ

x_{i j}p+f_{i j}p

)



(9) +cd  (i, j)∈A  p∈P d_{i j}px_{i j}p (10)

s.t.  j∈N x0j≤ K, (11)  i∈N xi j=1,

∀

j∈N

\

{

0

}

, (12)  j∈N xi j=1,

∀

i∈N

\

{

0

}

, (13)  i∈N fi j−  k∈N fjk=qj,

∀

j∈N

\

{

0

}

, (14) qjxi j≤ fi j≤

(

Q− qi

)

xi j,

∀

(

i,j

)

∈A, (15)  p∈P xp_{i j}=xi j,

∀

(

i,j

)

∈A, (16)  m∈M xpm_{i j} =xp_{i j},

∀

(

i,j

)

∈A,

∀

p∈P, (17)  p∈P f_{i j}p=fi j,

∀

(

i,j

)

∈A, (18) qjxi jp≤ f p i j≤

(

Q− qi

)

xi jp,

∀

(

i,j

)

∈A,

∀

p∈P, (19)

ϑ

pm 0j =Lx pm 0j,

∀

j∈N

\

{

0

}

,

∀

p∈P,

∀

m∈M, (20)

ϑ

pm i j +

ζ

pm i j =w pm i j ,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (21)

ζ

pm i j ≤ ui,

∀

i∈N,

∀

p∈P,

∀

m∈M, (22) bpm_{i j} xpm_{i j} ≤ wpm i j ≤ b p,m+1 i j x pm i j ,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (23)  i∈N  p∈P  m∈M



 r∈R

θ

pmr i j w pm i j +  r∈R

η

pmr i j x pm i j

+ i∈N  p∈P  m∈M w_{i j}pm+sj=  k∈N  p∈P  m∈M

ϑ

pm jk ,

∀

j∈N

\

{

0

}

, (24)  p∈P  m∈M



 r∈R

θ

pmr i0 w pm i0 +  r∈R

η

pmr i0 x pm i0

≤ U,

∀

i∈N

\

{

0

}

, (25)

xi j∈

{

0,1

}

,

∀

(

i,j

)

∈A, (26) x_{i j}p∈

{

0,1

}

,

∀

(

i,j

)

∈A,

∀

p∈P, (27) x_{i j}pm∈

{

0,1

}

,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (28) fi j≥ 0,

∀

(

i,j

)

∈A, (29) f_{i j}p≥ 0,

∀

(

i,j

)

∈A,

∀

p∈P, (30) w_{i j}pm≥ 0,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (31)

ϑ

pm i j ,

ζ

pm i j ≥ 0,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M. (32)

Inthe objectivefunction,expressions (7),(8),and(9)arethe engine,speed,andweightmodules ofthefuel cost,and expression(10)isthedistance-basedvehicledepreciationcost.cf_andcd_{arethecorrespondingcostcoefficients.}

Constraints(11)–(15)arestandardconstraintsofthetwo-indexsingle-commodityflowmodelofthecapacitatedvehicle routing problem (Gavish andGraves, 1981). Specifically, constraint (11) imposes that there are no more than K vehicles, constraints(12)and(13)arethevehicleflowconservationconstraints,constraints(14)arethecommodityflowconservation constraints,andconstraints(15)guaranteethecapacityofeachvehicleisnotexceeded.

Constraints (16)restrict thevehicletoselectexactlyonepath,iftravelingonarc(i,j)∈A,andconstraints(17)restrict the vehicle to depart at node i within exactly one time interval on the selected path _{Pi j}p. Constraints (18) and(19) en-surethatthe commodityiscarriedonexactlyone pathandvehiclecapacityisnotexceeded, ifthe vehicletravelson arc (i,j)∈A.

Constraints(20)specifythatthevehicledepartsfromthedepotnoearlierthanthebeginningofthetimehorizon. Con-straints(21) calculatetheactualdeparturetime atthedepot(i=0) oracustomer nodei ∈Nc_{ࢨ{0},iftimeflexibility (i.e.,}

flexibledeparturetime atthedepotand/orpost-servicewaitingtimeatcustomersites)isallowed.Constraints(22)ensure theflexibledeparturetimeatthedepot(i=0)andthepost-servicewaitingtimeatcustomersitei∈Nc_{ࢨ{0}donotexceed}

their correspondingupperboundsu_i.Constraints(23)restrict thedeparturetime atnode ion theselectedpath_{Pi j}p tobe within thetime interval asspecifiedin (17).If pathPp

i j isnot selected,w pm

i j =0,

∀

m∈M.Constraints (24)requirethat, if

the vehicle travelson path Pp

i j of arc(i, j), theearliest possibledeparture time atcustomer node j isthe summation of

theactualdeparturetime atcustomer nodei,thetraveltime onpath_{Pi j}p_,andtheservicetimeatcustomernodej,which ensures thetime continuity ofa vehicletravelinginthe network.Constraints (25)requirethat the delivery iscompleted withinthedeliverytimehorizon.

Last, constraints (26)–(32)are variable specification constraints. Note that constraints (15) become redundant due to constraints(18)and(19),butwekeepthemhereforclarity.

Whentime flexibilityisnot allowed,ui=0inconstraints(22),we have

ζ

i jpm=0andconstraints(21)reduce to

ϑ

pm i j =

wpm_{i j} .Wecanreplace

ϑ

_{i j}pmwithwpm_{i j} inconstraints(20)and(24),eliminatedecisionvariables

ϑ

_{i j}pm,andeliminateconstraints

(21),(22),and(32).

WhenthesetofcandidatepathsPijofeacharc(i,j)∈Acreducestoonepath,themodelTDVRP–PFreducestothe

time-dependent vehicle routing problem without path flexibility (TDVRP). We dropthe superscript pin the TDVRP–PFmodel. Thatis,x_{i j}p, f_{i j}p,x_{i j}pm,wpm_{i j} ,

ϑ

pm

i j ,and

ζ

pm

i j becomexij,fij,xmi j,wmi j,

ϑ

i jm,and

ζ

i jm,respectively.Consequently,we caneliminate

constraints (16), (18), (19),(27), and(30), and then drop the superscript p in the objective function and the remaining constraints.

5.2. StochasticTDVRP–PF

Todeal withstochastic traffic conditionsinthe roadnetwork,we formulatea two-stagestochasticmixedinteger pro-gram,combiningthevehicleroutingdecisionsbasedonstatisticalanalysisoftheexpectedtrafficpatternsandtheenroute pathselectiondecisionsbasedonreal-timetrafficconditions.

We model the stochastic traffic conditions as a finite set of scenarios



. Denote

ω

∈



as a random sce-nario and denote

ξ

(

ω

)

₌

(

θ

_{i j}pmr

(

ω

)

_,

η

_{i j}pmr

(

ω

)

_,

v

pr

i j

(

ω

)

,b pm

i j

(

ω

))

as the vector of the random variables representing the

traffic condition under scenario

ω

. In the two-stage stochastic program, the vehicle routing decisions x₌

(

x_{i j,}f_{i j}

)

are the first-stage decisions. Given the scenario

ω

∈



, the path selection and related traveling decisions y

(

ω

)

=

(

xp_{i j}

(

ω

)

,xpm_{i j}

(

ω

)

,f_{i j}p

(

ω

)

,w_{i j}pm

(

ω

)

,

ϑ

pm i j

(

ω

)

,

ζ

pm

i j

(

ω

))

are the second-stage decisions. The recoursedecisions thus involve the

Thetwo-stagestochasticmixedintegerprogram(TDVRP-PFS_{)isformulatedasfollows:} min E_ξ[h

(

x,

ξ

)

] (33) s.t.  j∈N x0j≤ K, (34)  i∈N xi j=1,

∀

j∈N

\

{

0

}

, (35)  j∈N xi j=1,

∀

i∈N

\

{

0

}

, (36)  i∈N fi j−  k∈N fjk=qj,

∀

j∈N

\

{

0

}

, (37) qjxi j≤ fi j≤

(

Q− qi

)

xi j,

∀

(

i,j

)

∈A, (38) xi j∈

{

0,1

}

,

∀

(

i,j

)

∈A, (39) fi j≥ 0,

∀

(

i,j

)

∈A, (40)

whereforeach

ω

_∈



, h

(

x,

ξ

(

ω

))

=min y(ω)c f



_ (i, j)∈A  p∈P  m∈M

ψ

e



 r∈R

θ

pmr i j

(

ω

)

w pm i j

(

ω

)

+  r∈R

η

pmr i j

(

ω

)

x pm i j

(

ω

)



(41) +  (i, j)∈A  p∈P  m∈M

ψ

s



 r∈R

(

v

pr i j

(

ω

))

3

θ

pmr i j

(

ω

)

w pm i j

(

ω

)

+  r∈R

(

v

pr i j

(

ω

))

3

η

pmr i j

(

ω

)

x pm i j

(

ω

)



(42) +  (i, j)∈A  p∈P

ψ

wd_{i j}p

(

μ

xp_{i j}

(

ω

)

+f_{i j}p

(

ω

))

(43) +cd  (i, j)∈A  p∈P d_{i j}px_{i j}p

(

ω

)

(44) s.t.  p∈P xp_{i j}

(

ω

)

=xi j,

∀

(

i,j

)

∈A, (45)  m∈M xpm_{i j}

(

ω

)

=xp_{i j}

(

ω

)

,

∀

(

i,j

)

∈A,

∀

p∈P, (46)  p∈P f_{i j}p

(

ω

)

= fi j,

∀

(

i,j

)

∈A, (47) qjxp_{i j}

(

ω

)

≤ f_{i j}p

(

ω

)

≤

(

Q− qi

)

xp_{i j}

(

ω

)

,

∀

(

i,j

)

∈A,

∀

p∈P, (48)

ϑ

pm 0j

(

ω

)

=Lx pm 0j

(

ω

)

,

∀

j∈N

\

{

0

}

,

∀

p∈P,

∀

m∈M, (49)

ϑ

pm i j

(

ω

)

+

ζ

pm i j

(

ω

)

=w pm i j

(

ω

)

,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (50)

ζ

pm i j

(

ω

)

≤ ui,

∀

i∈N,

∀

p∈P,

∀

m∈M, (51) bpm_{i j}

(

ω

)

x_{i j}pm

(

ω

)

≤ wpm i j

(

ω

)

≤ b p,m+1 i j

(

ω

)

x pm i j

(

ω

)

,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (52)  i∈N  p∈P  m∈M



 r∈R

θ

pmr i j

(

ω

)

w pm i j

(

ω

)

+  r∈R

η

pmr i j

(

ω

)

x pm i j

(

ω

)



+ i∈N  p∈P  m∈M w_{i j}pm

(

ω

)

+sj=  k∈N  p∈P  m∈M

ϑ

pm jk

(

ω

)

,

∀

j∈N

\

{

0

}

, (53)  p∈P  m∈M



 r∈R

θ

pmr i0

(

ω

)

w pm i0

(

ω

)

+  r∈R

η

pmr i0

(

ω

)

x pm i0

(

ω

)



≤ U,

∀

i∈ N

\

{

0

}

, (54) x_{i j}p

(

ω

)

∈

{

0,1

}

,

∀

(

i,j

)

∈A,

∀

p∈P, (55) x_{i j}pm

(

ω

)

∈

{

0,1

}

,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (56) f_{i j}p

(

ω

)

≥ 0,

∀

(

i,j

)

∈A,

∀

p∈P, (57) w_{i j}pm

(

ω

)

≥ 0,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M, (58)

ϑ

pm i j

(

ω

)

,

ζ

pm i j

(

ω

)

≥ 0,

∀

(

i,j

)

∈A,

∀

p∈P,

∀

m∈M. (59)

Theconstraintsaresimilartothedeterministicmodel(7)–(32),exceptthat thescenariosareintroducedtothe second-stage subproblems. The first-stage master problem(33)–(40) is a capacitated vehiclerouting problemwith theobjective functionconsideringtheexpectedcostofthesecond-stagesubproblems(41)–(59).Thestochasticprogram(TDVRP–PFS_)has

mixed-integervariablesinbothstages,andtherecoursematrixisrandom,whichmakesitverydifficulttosolve.

Whentimeflexibility(i.e.,flexibledeparturetimeatthedepotand/orpost-servicewaitingtimeatcustomersites)isnot allowed,thesimplificationoftheTDVRP–PFS _{issimilartothatinthedeterministicmodel.Whenthesetofcandidatepaths}

Pijofeacharc(i,j)∈Acreducestoonepath,themodelTDVRP–PFSreducestothestochastictime-dependentvehiclerouting

problemwithoutpathflexibility(TDVRPS_{),wherewedropthesuperscriptpinthesecond-stagesubproblems,similartothe}

conversionfromtheTDVRP–PFtotheTDVRP.However,thereisnoactualrecourseactioninthesecond-stagesubproblems, duetothelackofpathflexibility.Therefore,thesecond-stagesubproblemssolelyserveastheevaluationofthefirst-stage routingdecisionsunderstochastictrafficconditions.

InTable2,wegivetheproblemsizesofbothTDVRP–PFandTDVRP–PFS_{.Tobe specific,inthedeterministic TDVRP–PF,}

theintegervariablesconsistofxij,xi jp,andx pm

i j ,whilethecontinuousvariablesconsistoffij, fi jp,w pm i j ,

ϑ

pm i j ,and

ζ

pm i j .Inthe

TDVRP–PFS_{,theintegerandcontinuousvariablesare similar, exceptthatthescenarios areintroduced inthesecond-stage}

variables, i.e.,y

(

ω

)

=

(

x_{i j}p

(

ω

)

,x_{i j}pm

(

ω

)

,f_{i j}p

(

ω

)

,w_{i j}pm

(

ω

)

,

ϑ

pm i j

(

ω

)

,

ζ

pm

i j

(

ω

))

.Observethat the problemsizeof thesubproblem

ofonescenariointheTDVRP–PFS_{issimilartothatofthedeterministicTDVRP–PF,whichisalreadylarge.Therefore,evenif}

welimitthenumberofscenariosintheTDVRP–PFS_{,thesizeofthestochasticprogramisstillhuge.}

Table 2

Problem sizes of the TDVRP-PF and TDVRP-PF S .

# constraints # integer var. # continuous var. Deterministic O(| N|2 ·_{| P| ·| M|) O(| N|}2 ·_{| P| ·| M|) O(| N|}2 ·_{| P| ·| M|)} Stochastic Master problem O(| N|2 _{) O(| N|}2 _{) O(| N|}2 ₎

5km

Expressways Arterials Customers Depot

Fig. 7. Road network of the urban area in Beijing and candidate customer locations.

Table 3

Average speed of different road categories.

Expressways Arterials Residential roads 06 :0 0–07:0 0 48 km/h 32 km/h 24 km/h 07 :0 0–09:0 0 36 km/h 24 km/h 24 km/h 09 :0 0–13:0 0 41 km/h 28 km/h 24 km/h

6. Numericalexperiment

Inthissection,wefirstdescribetheexperimental settingsinSection6.1andthenprovidetheexperimental resultsand analysisinSection6.2.

6.1. Experimentaldesign

Clearly, instance design is important. Similar to Qian andEglese (2016) (for a London street network), we combined thetraffic network fromtheBeijingsituation andfolloweda similar approachasin Solomon(1987)forconstructing the customer graph. In what follows, we describe the construction of the geographic graph, based on the road network of Beijing,andthegenerationofatestbedofinstances.Further,wedescribetheexperimentalsettingsandthecomputational environment.

6.1.1. Networkconstruction

We construct geographical graph (road network) based on the urban area of Beijing (within the 5th ring road). We classify theroads into three categories(Wang etal., 2008): expressways, includingfour ring roads(named fromthe 2nd tothe 5thring road inthe orderof radial distancefrom thecity center) andurban freeways connectingthe ring roads,

arterials,includingarterialandsub-arterialroadsinthecity,andresidentialroads,referringtorelativelysmallerroads,other thanexpresswaysandarterials.

IntheurbanareaofBeijing,mostroadsarenorth-southoreast-west.Therefore,weonlydepictexpresswaysandarterials inthenetworkanduseManhattandistancestocalculatethetravelingdistancesbetweenlocationsalongresidential roads. ThenetworkofexpresswaysandarterialsinBeijingisshowninFig.7,consistingof673nodesand1870arcs.

Theaverage speedsof differentroadcategories are showninTable 3.Note that,due tothe spatialdifference of con-gestion,speeds of the same roadcategory vary in different areas ofthe city. The detailed speed setting is describedin

AppendixB.Notealsothatthetrafficconditionofresidentialroadsisnormallystableoverthedayandthespeedis there-foreassumed to be constant over time. Peoplenormally drive on theseresidential roads asthe first mileand last mile connectingtoexpresswaysorarterials.

Toconstructthecustomer graph,we first generatea candidatecustomer setof111 customer nodes(circles inFig. 7), whicharethelocationsofvariousretailstoressuchashypermarkets,fruitchainsstores,andhomeappliancestores.Inour