Computers & Operations Research 36 (2009) 919 – 935

www.elsevier.com/locate/cor

The petrol station replenishment problem with time windows

Fabien Cornilliera , Gilbert Laportea, b , Fayez F. Boctora , Jacques Renauda,∗
a Interuniversity Research Centre on Enterprise Networks, Logistics and Transportation (CIRRELT), Université Laval, Québec, Canada G1K 7P4
b Canada Research Chair in Distribution Management, HEC Montréal, 3000 chemin de la Côte-Sainte-Catherine, Montréal, Canada H3T 2A7

Available online 4 December 2007

Abstract
In the Petrol Station Replenishment Problem with Time Windows (PSRPTW) the aim is to optimize the delivery of several
petroleum products to a set of petrol stations using a limited heterogeneous fleet of tank-trucks. More specifically, one must
determine the quantity of each product to deliver, the assignment of products to truck compartments, delivery routes, and schedules.
The objective is to maximize the total profit equal to the sales revenue, minus the sum of routing costs and of regular and overtime
costs. This article first proposes a mathematical formulation of the PSRPTW. It then describes two heuristics based on arc preselection
and on route preselection. Extensive computational tests on randomly generated instances confirm the efficiency of the proposed
heuristics. Finally, a performance analysis on a real case shows a distance reduction of more than 20% over a solution obtained by
an experienced dispatcher.
䉷 2007 Elsevier Ltd. All rights reserved.

Keywords: Fleet management; Fuel delivery; Replenishment; Routing and scheduling

1. Introduction

This article proposes a mathematical model and two heuristics for the Petrol Station Replenishment Problem with
Time Windows (PSRPTW). This problem consists of optimizing the delivery of several petroleum products to a set of
petrol stations, which must be supplied once by a limited heterogeneous fleet of tank-trucks based at a terminal, within
given time windows. This problem is motivated by a real case faced by a Quebec based transportation company. In
North America most petrol companies subcontract the replenishment of their outlets to private transporters who are
paid on a delivered quantity basis. The objective of the PSRPTW is to maximize the total profit, equal to the revenue
which is a function of delivered quantities, minus the sum of routing costs and of regular and overtime costs. Decision
variables specify how much of each product to deliver to stations subject to their minimum requirements, how to assign
products to vehicle compartments, and how to design and schedule delivery routes. In this study, we formulate the
PSRPTW as a mixed integer linear program. We then propose two heuristics based on the same formulation. The first
heuristic consists of preselecting promising arcs and of solving the associated mathematical program to optimality. It
can solve small instances of up to about 15 stations. The second heuristic makes a preselection of promising routes
through a geographical decomposition method, and can be applied to larger instances.

∗ Corresponding author. Tel.: +1 418 656 7029; fax: +1 418 656 2624.
E-mail addresses: [email protected] (F. Cornillier), [email protected] (G. Laporte), [email protected] (F.F. Boctor),
[email protected] (J. Renaud).

0305-0548/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved.
doi:10.1016/j.cor.2007.11.007
920 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

A major difference between the PSRPTW and most vehicle routing problems with time windows is the loading
component: in the PSRPTW one must simultaneously design vehicle routes and assign petroleum products to truck
compartments for each trip. Related but different problems have been studied by a number of authors. Brown and
Graves [1] have considered the problem of direct deliveries (i.e. single-customer trips) and time windows. Different
algorithms have been proposed for other versions of the same problem without time windows. Brown et al, [2] and
Malépart et al. [3] have generalized this problem by allowing the delivery of more than one station in a same trip.
Heuristics for a multiperiod version of this problem have been developed by Taqa allah et al. [4]. In Cornillier et al.
[5], an exact algorithm is developed for a similar problem without time windows, where only one or two stations can
be delivered within a same trip, while Avella et al. [6] have proposed a heuristic and an exact algorithm based on a
route generation scheme and a branch-and-price algorithm to solve a similar problem. More recently, a heuristic for
the multiperiod problem without time windows was put forward by Cornillier et al. [7] for the case where the number
of stations on any given route is limited to two.
The remainder of this paper is organized as follows. The problem is defined and modelled in Section 2. The route
generation procedure is described in Section 3. In Section 4, we propose two heuristics, one based on arc preselection
and the other based on route preselection. Computational results are presented in Section 5 and conclusions follow in
Section 6.

2. Problem definition and formulation

The PSRPTW can be formally defined as follows. Let G = (V ∗ , A) = (V ∪ {0}, A) be a directed graph where
V = {1, . . . , n} is the set of stations, vertex 0 is the terminal, and A = {(i, j ) : i, j ∈ V ∗ , i  = j } is the arc set. Denote
by cij and tij the travel cost and the travel time associated with arc (i, j ), and by si the service time of station i. The
PSRPTW consists of maximizing a profit related function by designing delivery routes to replenish stations with a
limited heterogeneous fleet of K tank-trucks based at the terminal. Service at station i must start and end within a given
time window [ai , bi ] satisfying bi − ai si . A working day contains H regular working hours which can be extended
by using H  overtime hours. A regular wage rate applies to regular working time while a higher rate applies to overtime
hours. Only the hours effectively worked are paid, i.e. the hours from the beginning of the vehicle first trip to the return
of its last trip. The total variable cost is the sum of travel costs cij , regular and overtime wages. All trucks are assumed
to travel at the same speed. Moreover, each truck is subdivided into several compartments of known capacities and is
not equipped with a flow metre. Thus a given compartment can only be used to deliver one product to one station. Each
petrol station has a given number of capacitated underground tanks. The PSRPTW consists of determining:

• the quantity of each product to be delivered to each station, which should lie between a minimum and a maximum;
• the loading of these products into vehicle compartments;
• feasible delivery routes to these stations;
• the selected routes assignment to available trucks;
• the departure time of each truck trip,

in order to maximize a profit function.

2.1. Assumptions

The following assumptions are made:

• only one working day is considered;

• the fleet is heterogeneous and limited;
• each station must be visited once and only once during the considered working day;
• since compartments are not equipped with flow metre, they must be entirely emptied once replenishment has
started;
• several trips can be assigned to the same truck;
• each station requires delivery within a given time window;
• waiting time between stations is allowed;
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 921

• regular and overtime working hours are limited;

• regular and overtime wages are known and constant;
• only effectively worked hours are paid;
• the transporter is paid a given amount for each litre delivered which varies as a function of station location;
• the travel time between any two vertices (terminal and stations), service times at stations and loading times at the
terminal are known;
• each station requires between a minimum and a maximum quantity of one or more products that can be computed
from initial inventories, expected consumptions, and the capacities of its underground tanks.

These assumptions correspond to the actual delivery practices in eastern Canada.

2.2. Mathematical formulation

Our mathematical model is based on the generation of all feasible routes a truck can follow. A route is feasible if it
satisfies all time windows and constraints on delivered amounts. Given a route r, one can compute its earliest and its
latest departure times, denoted by r and r , minimizing its total waiting time.
We first define the following parameters:

 regular wage per hour;

 overtime wage per hour;
r earliest departure time for route r;
r latest departure time for route r;
r minimum duration of route r (including waiting time if any);
asr a binary parameter equal to 1 if and only if station s is delivered within route r;
rk the profit of route r if performed by truck k. This parameter is equal to −∞ if truck k is unable to carry out
route r.
The decision variables are

xrkv a binary variable equal to 1 if and only if route r corresponds to trip v of truck k;
dkv the departure time of truck k for trip v;
hk the number of regular working hours of truck k;
hk the number of overtime hours of truck k.
The model is then
  
(PSRPTW) Maximize rk xrkv −  hk −   hk , (1)
(r,k,v) k k

subject to

asr xrkv = 1 ∀s, (2)
(r,k,v)

xrkv 1 ∀(k, v), (3)
r
 
xrk,v+1  xrkv ∀(k, v), (4)
r r
 
  
r xrkv dkv  r xrkv + M 1 − xrkv ∀(k, v), (5)
r r r

dk,v+1 dkv + r xrkv ∀(k, v), (6)
r
922 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

Fig. 1. Illustration of the model with two tank-trucks and three routes.


dkv − dk,1 + r xrkv hk + hk ∀(k, v), (7)
r

hk H ∀k, (8)
hk H  ∀k, (9)
+
dkv ∈ R ∀(k, v), (10)
hk ∈ R+ , hk ∈ R+ ∀k, (11)
xrkv ∈ {0, 1} ∀(r, k, v). (12)
In this formulation, the objective function (1) maximizes the total profit. Constraint (2) stipulates that each station
is visited once and only once. Constraint (3) ensure that at most one route is assigned to the vth trip of truck k. By
constraint (4), trip v + 1 of truck k exists only if trip v exists. In constraint (5), M is a large positive number; these
constraints require that routes departure times lie within the computed windows [r , r ]. Constraint (6) states that
trip departure occurs after the arrival time of the preceding trip. Since we do not know the number of trips of truck
k, constraint (7) ensures that the total working hours, decomposed into regular and overtime hours, is equal to the
duration, calculated as the difference between its latest trip return time and its first trip departure time. Constraint (8)
and (9) ensure that regular and overtime hours lie within the allowable limits.
Fig. 1 illustrates a case with three routes r1 , r2 , and r3 , where [1 , 1 ] = [1, 2], [2 , 2 ] = [1, 3], [3 , 3 ] = [0, 3],
and 1 = 1, 2 = 2, and 3 = 3. There are two identical trucks generating the same profit, and the maximal working
hours are H = 3 and H  = 2, for a total of five hour per truck. Fig. 1a depicts an optimal solution using one overtime
hour, while Figs. 1b and 1c correspond to suboptimal solutions using two overtime hours. In Fig. 1d, we show that the
solution of Fig. 1a cannot be improved by starting route r1 at time t = 0 because r1 would then need one additional
waiting hour.

3. Route generation

In the above formulation, the number of potential routes is generally huge. We first propose ways of reducing the
number of routes through feasibility and dominance criteria. A route can potentially visit as many stations as there are
compartments in the truck. However, a two-stop limit per route is common practice in North America. This is explained
by the fact that most trucks have from four to six compartments, while stations generally require two or three products,
one of which frequently requires two compartments. In this study, we consider the case where routes can visit up to four
stations. If G is a complete graph, the number of feasible and infeasible routes visiting at most m stations is equal to
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 923

m
i=1 n!/(n − i)! and can be rather large. Instead of making an explicit enumeration of all feasible and infeasible routes
to be checked, we use an adaptation of Johnson’s algorithm [8] to generate candidate routes from G, or a subgraph of
G, from which all infeasible arcs are removed. These are arcs that cannot be included in a solution without violating a
time or duration constraint. Given a directed graph, this algorithm consists of enumerating all or some of its elementary
circuits. In our adaptation, routes are iteratively built starting from the terminal, station by station, until no more stations
can be added without violating time or quantity constraints.

3.1. Infeasible arc deletion

In our problem, some infeasible arcs of G can be removed since some station pairs are incompatible in terms of
time windows or requested quantities. Because these stations cannot belong to the same route, the number of feasible
routes is reduced. We define the subgraph G = (V ∗ , A ) of G where each arc of A corresponds to a pair of compatible
stations with respect to their time windows. Also, for each truck k we define the subgraph Gk = (V ∗ , Ak ) of G where
each arc of Ak corresponds to a pair of compatible stations with respect to their time windows and demand feasibility
constraints. Demand feasibility of a route for a given truck is checked as shown in Section 3.2.

3.2. Demand feasibility check and quantity determination

A feasible route should allow the delivery of all minimal quantities required by its stations, and should visit these
stations within the required time windows. In this section, we solve a tank-truck loading problem (TTLP) defined as
follows. Let P be the set of demands of all stations on the route, and let gp be the revenue associated with quantity
qp delivered to station p. This revenue is a function of the distance between the station at which the delivery takes
place and the terminal. The TTLP consists of determining the quantity qp to be delivered to each station p of the
route in order to maximize the sum of revenues, while respecting the minimal and maximal requirements up and vp ,
and without exceeding the capacity of any tank-truck compartment. Related problems using compartmented vehicles,
generally referred to as Loading Problems, have been addressed with different objectives (Christofides et al. [9], Yuceer
[10], Smith [11], Bukchin and Sarin [12]), and in different applications: bulk ship scheduling with flexible cargo holds
(Fagerholt and Christiansen [13,14]), livestock transportation (Oppen and LZkketangen [15]), grocery delivery (Eglese
et al. [16]), and oil delivery (Brown et al. [2], Van der Bruggen et al. [17], Bausch et al. [18]).
In the PSRPTW, the loading problem can be formulated as follows. Let ypc be a binary variable equal to 1 if demand
p is assigned to compartment c, and 0 otherwise. Then the problem is

(TTLP) Maximize gp q p , (13)
p∈P

subject to

up qp vp ∀p, (14)


qp  Qc ypc ∀p, (15)
c

ypc 1 ∀c, (16)
p∈P

YP C ∈ {0, 1} ∀(p, c). (17)

In this formulation, the objective function maximizes the total revenue. Constraint (14) ensure that the delivered
quantities lie between the requested minimum and maximum. Constraint (15) states that delivered quantity of demand
p cannot be larger than the sum of compartment capacities in which it is loaded. By constraint (16), two distinct
demands cannot be loaded in the same compartment.
This model is used to check the feasibility of each route with respect to a given truck and to obtain an optimal
load by maximizing the corresponding revenue. Once the delivered quantities are known, one can compute the service
time of each visited station and the profit rk for route r and truck k, which is equal to the difference between the
924 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

revenue generated by the delivered quantities and the route travel cost computed as the sum of the corresponding cij .
By convention, we set rk = −∞ if truck k is unable to deliver the requirements of route r.

3.3. Route duration and departure window

Since the aim is to select a subset of feasible routes and to determine their optimal truck assignments and schedules,
we must compute for each of these a time interval within which any departure time from the terminal minimizes the
total duration including service time and waiting time, if any. Given a route r delivering all stations of the subset
Vr ⊆ V , we index its stations according to the sequence in which they must be visited. Denote by Vr∗ = Vr ∪ {0} the
set of vertices including the terminal and all stations of route r. We check whether we can satisfy the time window
constraint of each station and if so, we determine the departure window [r , r ] for which the sum of waiting times
is minimal. Savelsbergh [19] has proposed an algorithm for the determination of r . It computes a forward time slack
which indicates by how much the departure time of a vertex can be delayed without making the route infeasible, and
iteratively computes the waiting times at each vertex. In our case, we need to determine the time interval [r , r ] within
which any departure time from the terminal minimizes the total duration. We then have to compute the sum wr of all
its necessary waiting times in order to determine its duration r . First, for each vertex i, we define a normalized time
window [ai , bi ] representing the time interval within which the truck should leave the terminal in order to satisfy the
time window constraint of station i if waiting times were not allowed:
 
i−1 
i−1 
i−1 
i
 
[ai , bi ] = ai − tu,u+1 − su , bi − tu,u+1 − su . (18)
u=0 u=0 u=0 u=0

If waiting times were not allowed, the route would be feasible if and only if the intersection of all normalized time
windows was not empty. However, waiting times may be needed in our problem and a new feasibility criterion is given
by Proposition 1. The proofs of all propositions are given in Appendix.

Proposition 1. If waiting times are allowed, a route is feasible if and only if

max {aj } bi , ∀i ∈ Vr . (19)

0  j <i

When a route r is feasible, we compute the sum wr of all its minimal waiting times by means of Proposition 2. This
waiting time is added to the sum of travel and service times in order to arrive at the route duration r .

Proposition 2. If a route r is feasible, the sum of its minimal waiting times wr is


wr = max 0, max∗ {ai } − min∗ {bi } . (20)
i∈Vr i∈Vr

Finally, we need to determine a departure window for the route r from the terminal, such that the time window
constraint of each station is satisfied.

Proposition 3. If a route r is feasible, its departure time d0 from the terminal has a time window [r , r ] which
minimizes the total waiting time, where

r = max∗ {ai } − wr (21)

i∈Vr

and

r = min∗ {bi }. (22)

i∈Vr

Note that starting route r at any time before r only increases the embedded waiting times and does not allow the
truck to return to the terminal earlier.
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 925

Fig. 2. Three stations route with infeasible time windows.

Fig. 3. Three stations feasible route without waiting time.

Fig. 4. Three stations feasible route with waiting time.

Proposition 1 is illustrated in Fig. 2 where b1 < max{a1 , a2 , a3 } = a2 ; in this case there is no feasible solution. In
Fig. 3, max{a1 , a2 , a3 } < min{b1 , b2 , b3 } and consequently, as implied by Propositions 2 and 3, wr = 0 and r r .
Fig. 4 shows the case where wr > 0 and r = r . In this case, the truck should leave the terminal exactly at r in order
to minimize its waiting time wr .
Knowing the total minimum waiting time of a route r, we are able to compute its total duration r , including travel
and service times:

n−1
r = wr + (tu,u+1 + su+1 ) + tn,0 . (23)
u=0

3.4. Route dominance

When several routes visit the same subset of stations but in a different order, we only retain Pareto optimal routes.
More precisely, given two routes r1 and r2 , both visiting the same set of stations with the same truck, route r1 can be
eliminated if 1 2 , 1 2 , 1 2 , and 1 2 .

4. Heuristics

As the number of stations grows, the problem becomes more difficult to solve, even if there are fewer candidate
routes, since the number of vehicles increases proportionally. In practice, it is often difficult to solve problems to
optimality with more than 15 stations. For larger problems, we propose two heuristic procedures in which we solve
the proposed mathematical model with only a preselected subset of all feasible routes instead of the whole set. The
aim of the first heuristic is to reduce the number of routes by preselecting a subset of all feasible arcs of Gk . It can
be applied to relatively small instances. In the second heuristic, we decompose the geographical space in order to
iteratively construct a candidate set of locally optimal routes which is then used to solve the global problem. This
heuristic is more appropriate for larger instances.
926 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

4.1. A heuristic based on arc preselection

The first heuristic preselects an arc subset Ak of Ak . In the first version of the heuristic, this subset includes the arcs
linking each vertex to its  nearest neighbours, where  is a parameter to be determined. Note that  is at most n − 1
because a station has n − 1 neighbours. In the second version, the arc subset includes all arcs of at most  successive
minimum spanning tree, where  is a parameter; this way of reducing the arc set is inspired by the work of Helsgaun
[20,21] for the Traveling Salesman Problem and of Toth and Vigo [22] for the Vehicle Routing Problem. The procedure
first generates a minimum spanning tree on the initial graph. It then removes the selected edges and repeats itself as
long as the graph is connected. The value of  can be at most
n/2 because the spanning trees sequential generation
procedure uses n−1 of the n(n−1)/2 potential edges for each tree. In addition to these selected arcs, all arcs linking the
terminal to customers in both directions are included. Nearest neighbours and minimum spanning trees are not based
on distances, but on travel times. Since minimum spanning trees are constructed on undirected graphs, we construct
them while setting the value of each edge equal to the minimum travel time for each of the two directions. All routes
are generated from each Ak as described in Section 3 and the optimal routes selection is made by solving the proposed
PSRPTW model.
If the parameters  or  are too small, the problem may be infeasible. On the other hand, for large instances and larger
values of the parameters, the number of generated routes can become prohibitive and make the formulation unsolvable.
We found that this heuristic becomes inefficient as soon as n reaches 20 when tested on randomly generated instances.
However, a real case experiment where 42 stations need to be replenished suggests that the arc preselection heuristic
can be efficiently used on larger real instances.

4.2. A decomposition heuristic based on route preselection

To solve larger instances, we propose a decomposition of the geographical space into sectors. Since any given decom-
position would be arbitrary, we consider successive random partitions. Each generated random sector s corresponds to
a different subset of stations V s ⊂ V . The decomposition is such that no sector appears in two different partitions. The-
oretically, any partition of V could be used, but since in practice distances are Euclidean, we have only used partitions
induced by non-overlapping sectors centred at the terminal. Once the partitions are generated, the problem associated
with each sector is solved exactly as a separate PSRPTW and, each time, the corresponding optimal routes are added
to the preselected route set. Thus, the idea of this decomposition scheme is to generate a set of locally optimal routes
which will be used to define the decision variables of the whole problem model.

4.2.1. Sector generation

Each sector includes a random number of stations between 5 and 10, so that the associated problem can easily and
quickly be solved to optimality while allowing the generation of good locally optimal routes. Each time a partition is
generated, we make sure that none of its sectors has previously been selected. We iteratively generate new partitions
until a given limit of preselected routes is reached, or until a given number of iterations have been executed. Note
that in a non-Euclidean space, we would have to partition the set of stations in a different way, based for example, on
a measure of geographical and time windows distances. The rest of the method would otherwise be identical.

4.2.2. Optimal routes for a given sector

For each sector, we solve the corresponding PSRPTW to optimality by means of a branch-and-bound algorithm
in order to generate a set of preselected routes. Since identical routes can be generated from different partitions, an
exponential penalty on the number
r of times a preselected route r has appeared in previous partitions is added to the
objective function in order to prevent cyclic generation of routes from one partition to another. This penalty is equal to
r (exp(
r /2) − 1), where r is proportional to the length of route r. The penalized objective function of the subproblem
is then
  
Maximize (rk − r (e
r /2 − 1))xrkv −  hk −   hk . (24)
(r,k,v) k k
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 927

Fig. 5. The route preselection heuristic using three successive partitions.

To solve the subproblem associated with a sector V s , we generate all its feasible routes using a restricted fleet in which
K|V s |/n trucks are randomly chosen from the whole fleet. Each time a subproblem is solved, we add all of its
optimal routes not already included to the set of preselected routes, up to the given limit of routes.

4.2.3. Recomposition procedure

After locally optimal routes have been extracted for all generated sectors, the recomposition procedure consists of
determining the best routes in order to obtain a global solution to the whole PSRPTW. The resulting route selection
problems (1)–(12) is solved using the entire fleet.
Fig. 5 illustrates the route preselection heuristic. Fig. 5a shows a first partition. The problems associated with all
sectors (1.1)–(1.4) are independently solved to optimality, and the resulting optimal routes are added to the preselected
routes set. Sectors (1.1) and (1.2) each give three locally optimal routes, and sectors (1.3) and (1.4) each yield only one
route. We generate a second partition (Fig. 5b) and a third one (Fig. 5c) which, respectively, give five and three new
locally optimal routes to add the preselected routes set (r9 .r14 for the second partition, and r15 .r17 for the third). After
three successive partitions, there are 16 routes in the preselected routes set. The problem is then to select the best routes
from the preselected route set in order to visit each station once and only once. In this example, the solution (Fig. 5d)
uses five routes from the first partition (r3 .r7 ), and three from the third one (r15 .r17 ).

5. Computational results

The two heuristics just described were coded in Objective-C and used the CPLEX 10.0 Callable Library. They were
run on dual AMD Opteron 250 processors 2.4 GHz computers with Linux operating system. We present two sets of
928 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

Table 1
Daily sales distribution

Category Daily sales (litres) Percentage (%)

1 0–1350 21.7
2 1350–2700 22.6
3 2700–5400 29.8
4 5400–8100 13.6
5 8100–10 800 6.2
6 10 800–16 200 6.1

Table 2
Underground tank typical configurations as a function of daily sales

Daily sales (l) Tank Tank size (l)

0–2700 1 25 000
2 15 000
3 15 000

2700–8100 1 35 000
2 22 700
3 25 000

8100–16 200 1 50 000

2 25 000
3 35 000

Table 3
Configurations of the tank-trucks

Type Total capacity (1000 l) Number of compartments Capacities (1000 l)

1 60 6 17, 6, 10, 10, 7, 10

2 54 5 16, 6, 6, 10, 16
3 50 4 16, 8, 12, 14

tests. We have first solved a set of randomly generated 15 stations instances in order to evaluate the performance of
the arc preselection and route preselection heuristics. The route preselection heuristic was then assessed on randomly
generated instances with 50 stations. Finally, we solved a real case where 42 stations need to be replenished. All
instances are available on the website http://www.fsa.ulaval.ca/personnel/renaudj/Recherche/PSRPTW.

5.1. Test instances

We have generated test instances similar to real-life problems using a set of real data extracted from Malépart et
al. [23]. From these data, we have determined a discrete random distribution on six categories of stations in function
of their daily sales (Table 1). Station categories are randomly drawn from this distribution and daily sales are then
randomly determined within the lower and upper limits of the obtained categories.
The sales of regular, intermediate, and super petrol grades are 76%, 7%, and 17% of the total, respectively. Because
daily sales and underground tank capacities are generally correlated, we present in Table 2 the typical observed tank
capacities as a function of the total daily sales. For our test instances, the underground tanks configuration for each station
is randomly selected among these three typical configurations. However, the probability of choosing the configuration
corresponding to the station daily sales is 80% and the probability of choosing each of the other configurations is 10%.
We consider three tank-truck configurations among those commonly used in practice (Table 3). The terminal co-
ordinates are (50,50) for all instances, while stations coordinates are integer and randomly drawn from a uniform
distribution in the 100 km × 300 km Euclidean space. The fleet compositions as a function of the problem size are
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 929

Table 4
Fleet compositions

Number of stations Type I Type II Type III Fleet size

15 2 2 1 5
50 8 5 5 18

Table 5
Per litre revenue as a function of the distance from terminal

Distance Revenue per delivered litre

0–50 $0.004
50–100 $0.007
100–150 $0.010
150–200 $0.013
> 200 $0.016

given in Table 4. In all tests, the number of stations that can be visited in a same route is limited to four, except in tests
in which we evaluate the impact of this limit.
We have used the following data for all instances:
driver wage per regular working hour: $15.00;
overtime hourly cost: $30.00;
variable travel cost per kilometre: $1.70;
average travel speed (km/h): 60;
truck loading time (min): 15;
station delivery time (min): 30;
daily regular working hours: 9;
daily maximum overtime hours: 3.

The revenue per delivered litre is a function of distance from the terminal. Rates are given in Table 5.

5.2. Performance of the proposed heuristics

In this section, we study the performance of the proposed heuristics. We first analyse the results given by the arc
preselection heuristic used. We then evaluate the impact of limiting the number of delivered stations per route. Finally,
we analyse the performance of the route preselection heuristic. Average results over 20 instances of 15 or 50 stations
are given.

5.2.1. Performance of the arc preselection heuristic

In Table 6, we evaluate the performance of the arc preselection heuristic which uses  ∈ {3, . . . , 6} nearest neighbours
on instances of 15 stations. The last row corresponds to the case where all arcs are selected and the solution is therefore
optimal. We can see that with three nearest neighbours the profit of 504.24 is 98.88% of the optimum, and an optimal
solution is found 13 times out of 20. With  = 3, the arc preselection heuristic generates only 299 of all 3060 feasible
routes, i.e. it eliminates 90.2% of all feasible routes limited to four stations; it also reduces computation time by 88.6%,
from 350 to 40 seconds. A tangible improvement can be observed when four nearest neighbours are considered: we
then attain 99.2% of the optimal profit while eliminating 81.5% of all feasible routes. An optimal solution is found in
90% of the cases. Further marginal improvements are obtained by using a larger number of nearest neighbours.
Table 7 shows the average results of the arc preselection heuristic using  ∈ {3, . . . , 6} successive minimum spanning
trees. We observe that the arc preselection heuristic based on the computation of three successive minimum spanning
trees yields a profit equal to 99.3% of the optimum while eliminating 79.7% of all feasible routes. Further improvements
are obtained if five successive minimum spanning trees are generated, yielding a profit equal to 99.83% of the optimum.
930 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

Table 6
Average results as a function of the number of nearest neighbours

 #CR Dist Qty Rev RT OT #Rtes %Rtes(s) Profit #O CPU

1 2 3 4

3 299 2067 568 485 4743 41.36 3.46 10.75 70.7 20.5 7.4 1.4 504.24 13 40
4 565 2051 564 712 4711 41.18 3.35 10.65 69.5 21.6 7.5 1.4 505.69 17 61
5 888 2052 564 705 4713 41.21 3.34 10.65 69.5 21.6 7.5 1.4 506.18 18 92
6 1271 2052 564 705 4713 41.21 3.34 10.65 69.5 21.6 7.5 1.4 506.18 18 155

All 3060 2047 563 282 4709 40.93 3.51 10.60 68.9 22.2 7.5 1.4 509.94 20 350

, number of nearest neighbours; #CR, number of preselected routes; Dist, distance travelled; Qty, delivered quantity in litres; Rev, revenue; RT,
regular hours used; OT, overtime hours used; #Rtes, number of selected routes in the solution; %Rtes(s ), percentage of routes visiting s stations;
Profit, profit corresponding to the best solution; #O, number of times the optimal solution has been obtained; CPU, computing time in seconds; all,
all arcs are selected.

Table 7
Average results as a function of the number of minimum spanning trees

 #CR Dist Qty Rev RT OT #Rtes %Rtes(s) Profit #O CPU

1 2 3 4

3 621 2052 564 705 4713 41.21 3.34 10.65 69.5 21.6 7.5 1.4 506.18 18 54
4 1317 2052 564 705 4713 41.21 3.34 10.65 69.5 21.6 7.5 1.4 506.18 18 156
5 2080 2050 564 530 4714 41.04 3.47 10.65 69.5 21.6 7.5 1.4 509.06 19 243
6 2589 2050 564 530 4714 41.04 3.47 10.65 69.5 21.6 7.5 1.4 509.06 19 287

All 3060 2047 563 282 4709 40.93 3.51 10.60 68.9 22.2 7.5 1.4 509.94 20 350

, maximum number of generated minimum spanning trees.

Table 8
Average results as a function of the maximal number of stations per route

#S/R #CR Dist Qty Rev RT OT #Rtes %Rtes(s) Profit #O CPU

1 2 3 4

2 140 2223 585 279 5025 42.67 4.83 11.05 64.3 35.7 – – 460.15 5 182
3 807 2062 564 979 4732 41.44 3.27 10.65 68.5 22.1 9.4 – 507.01 17 113
4 3060 2047 563 282 4709 40.93 3.51 10.60 68.9 22.2 7.5 1.4 509.94 20 350

#S/R, maximal number of stations per route.

This procedure uses 30.6% less than the CPU time needed to obtain an optimal solution. From Tables 6 and 7, we can
see that there is no tangible performance difference between the two versions of the arc preselection heuristic for a
similar number of candidate routes.

5.2.2. Impact of limiting the number of delivered stations per route

It is possible to reduce the number of generated routes by reducing the maximal number of stations per route. In this
section, we evaluate the impact of this parameter. Average results are presented in Table 8. A significant improvement
over the common practice discussed in Section 3, which consists of limiting to two the number of stations per route,
can indeed be obtained by increasing this limit. A relatively large profit improvement of 10.2% is obtained when we
increase the limit to three stations per route. A marginal profit improvement of 0.58% is obtained by further raising
this limit to four stations. We note that the optimal profit with a limit of three stations per route is slightly better than
that obtained by the arc preselection heuristic with four successive minimum spanning trees or six nearest neighbours,
when limiting the number of stations per route to four.
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 931

Table 9
Average results as a function of the number of preselected routes for the 15 stations instances

Dist Qty Rev RT OT #Rtes %Rtes(s) Profit #O CPU

1 2 3 4

45 2077 564 455 4752 41.15 3.80 10.65 68.1 23.9 7.0 0.9 490.56 11 131
90 2082 567 680 4764 41.25 3.81 10.70 68.2 24.3 6.5 0.9 491.08 12 241
135 2066 564 765 4742 41.08 3.69 10.65 69.5 21.6 7.5 1.4 502.72 14 308
180 2070 568 161 4753 41.25 3.60 10.70 69.6 22.0 7.0 1.4 507.32 18 472
225 2070 568 161 4754 41.25 3.60 10.70 69.6 22.0 7.0 1.4 508.54 19 579
270 2070 568 161 4754 41.25 3.60 10.70 69.6 22.0 7.0 1.4 508.54 19 792
315 2070 568 161 4754 41.25 3.60 10.70 69.6 22.0 7.0 1.4 508.54 19 954

3060 2047 563 282 4709 40.93 3.51 10.60 68.9 22.2 7.5 1.4 509.94 20 350

, maximum number of preselected routes.

Table 10
Average results as a function of the number of generated routes for the 50 stations instances

Dist Qty Rev RT OT #Rtes %Rtes(s) Profit

1 2 3 4

150 7003 1 863 858 16 075 141.27 9.66 34.55 60.3 34.6 5.1 0.0 1 760.54
300 6862 1 849 960 15 843 138.91 9.57 34.20 60.1 33.9 5.7 0.3 1 807.58
450 6686 1 830 191 15 531 135.64 9.81 33.80 60.5 32.0 6.7 0.9 1 835.91
600 6632 1 818 840 15 438 134.62 9.86 33.55 61.1 30.1 7.5 1.3 1 848.42
750 6542 1 806 543 15 256 133.81 9.11 33.30 61.0 29.6 7.8 1.7 1 853.71
900 6571 1 808 900 15 307 134.81 8.60 33.35 61.2 29.4 7.8 1.6 1 856.91

5.2.3. Performance of the route preselection heuristic

To evaluate the performance of the route preselection heuristic, we set as a multiple of the instance size: ∈
{45, 90, . . . , 315} for the 15 stations instances, and ∈ {150, 300, . . . , 900} for the 50 stations instances. In the last
row, all feasible routes are generated, yielding the optimum. Table 9 shows the average results obtained with the route
preselection heuristic for the 15 station instances. For three preselected routes per station (45 routes, i.e. 1.47% of
all feasible routes), the profit is about 96.2% of the optimum. The largest improvement is obtained between 6 and 12
preselected routes per station (90 and 180 routes, i.e. 2.94% and 5.88%) with a profit of about 99.5% of the optimum.
Table 10 shows the average results of the route preselection heuristic for instances with 50 stations, which corresponds
to an average transporter working day. For these instances, optimal solutions are not available as we never succeed
to solve the model with more than 25 stations. For each instance, the allowed computation time was limited to two
hours of CPU time. We can see that a major profit improvement of 4.9% can be obtained by increasing the number of
preselected routes per station from three to nine (150 to 450 routes). A slight improvement can be obtained by further
increasing this number, but we observe that the profit and the MIP best bound level off (Fig. 6). When going from 15
to 18 preselected routes per station (750–900 routes), the profit improvement is only 0.17%. These results show that
this route preselection heuristic is capable of generating a set of good routes and can solve much larger instances than
any of the two versions of the arc preselection heuristic. In Tables 9 and 10, the number of generated routes #CR is not
displayed as it is always equal to .

5.3. Performance analysis on a real case

We have tested the arc preselection heuristic on a real-life instance arising in eastern Canada with a depot located in
Quebec City. The delivery data correspond to 42 stations on a single regular day. We used the drivers’ worksheets to
determine the routes and delivered quantities of each product for that day. We define the profit rk of route r performed
932 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

Best solution
MIP best bound

1880

1860
Average profit
1840

1820

1800

1780

1760
150 300 450 600 750 900
Number of generated routes

Fig. 6. Average profits as a function of the number of generated routes for the 50 stations instances.

Table 11
Time windows repartition

Number of stations ai bi

11 0h00 24h00
1 4h30 23h00
1 6h00 22h00
4 6h00 23h00
2 6h00 24h00
1 6h30 23h00
1 6h30 24h00
2 7h00 18h00
1 7h00 21h30
6 7h00 22h00
5 7h00 23h00
7 8h00 18h00

by truck k as the negative value of its distance in kilometres, and we set both  and  to zero which enable us to directly
apply models (1)–(12) and the heuristic developed in this study. The time windows [ai , bi ] for servicing stations have
been obtained by the company and are displayed in Table 11.
The manual solution obtained by the company dispatcher (which includes routing and loading of the products) has
a length of 7827 km and is composed of 16 routes visiting two stations, and 10 routes visiting a single station each. We
solved this problem with the heuristic based on arc preselection presented in Section 4.1. Table 12 details the results
obtained by the company and by the heuristic when the arc preselection is done both with the  nearest neighbours and
with the  successive minimum spanning trees. In comparison, solving this instance exactly with CPLEX requires the
generation of all the 17 720 feasible routes and more than 3.5 hours of computing time, including 3 hours for the route
generation alone.
The best results are obtained with the arc preselection procedure based on the minimum spanning tree, irrespective of
the  value. In this case, the heuristic produces within a minute an optimal solution of 6108 km, and an improvement of
22% (1719 km) over the solution obtained manually. Not only is this solution shorter, but more products are loaded into
the vehicle compartments. For the minimum spanning tree with  = 3, the vehicles delivered 1 184 000 litres, 16 500
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 933

Table 12
Detailed solution of the real case

#CR Dist(km) Qty(L) L/km #Rtes %Rtes(s) CPU

1 2 3 4

Manual 7827 1 167 500 149.2 26 10 16 0 0

=3 167 6777 1 185 000 174.9 25 12 10 2 1 12

=6 338 6673 1 186 000 177.7 23 8 11 4 0 29
=9 1200 6352 1 184 000 186.4 22 7 11 3 1 215
 = 12 1581 6352 1 184 000 186.4 22 7 11 3 1 293
 = 16 2386 6348 1 184 000 186.5 22 7 11 3 1 663

=1 124 6936 1 183 000 170.6 27 15 10 1 1 5

=2 330 6094 1 181 000 193.8 23 7 14 1 1 19
=3 610 6108 1 184 000 193.8 23 7 14 1 1 59
=4 944 6108 1 184 000 193.8 23 7 14 1 1 123
=5 1584 6108 1 184 000 193.8 23 7 14 1 1 250
=6 2336 6108 1 184 000 193.8 23 7 14 1 1 411

All 17 720 6108 1 184 000 193.8 23 7 14 1 1 12 885

all, all feasible routes.

litres more than the manual solution, which also confirms that combining a good loading with a good vehicle routing
is not easy for a dispatcher. The quantity delivered per kilometre is 193.8 litres with the heuristic, compared to 149.2
litres for the manual solution. It is worth observing that in the heuristic solution routes visiting three and four stations
are used, which enables a reduction in the total number of routes from 26 to 23, while still making larger deliveries.
We did not test the route preselection heuristic on the real case because we only had the distance matrix as input
data, as opposed to the station coordinates which are necessary to generate sectors.

6. Conclusions

We have proposed a mathematical formulation of the petrol stations replenishment problem with time windows.
Based on this formulation, an arc preselection heuristic was developed in order to reduce the number of candidate
routes. Computational results show that this heuristic reduces computation time while yielding near-optimal solutions.
For larger instances, a decomposition heuristic based on route preselection was proposed. On small instances, it was
compared to an exact algorithm. Computational results show that this decomposition heuristic succeeds in finding
near-optimal solutions while using a very small proportion of all feasible routes. Moreover, the effect of generating
more routes was analyzed on larger instances.

Acknowledgements

This work was partially supported by the Canadian Natural Sciences and Engineering Research Council (NSERC)
under Grants OGP0036509, OGP0039682, and OGP0172633. This support is gratefully acknowledged. Thanks are
due to José Eduardo Pécora Junior, to two referees, and the editor for their valuable comments.

Appendix
 i i
Proof of Proposition 1. Let Ti = i−1 u=0 tu,u+1 , Si = u=0 su , and Yi = u=0 yu , where yi 0 denotes a minimal
waiting time between stations i and i + 1. Then ai = ai − Ti − Si−1 and bi = bi − Ti − Si . The route is feasible if and
only if for each station i there exists a departure time, denoted by di ∈ R, such that

∀i ∈ Vr∗ , ai + si di bi , (25)

934 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935

which is equivalent to

∀i ∈ Vr∗ , ∃di ∈ R : ai + Ti + Si−1 + si di bi + Ti + Si (26)

⇔ ∀i ∈ Vr∗ , ∃di ∈ R : ai + Ti + Si di bi + Ti + Si (27)

⇔ ∀i ∈ Vr∗ , ∃di ∈ R : ai di − Ti − Si bi . (28)

But as di = d0 + Ti + Si + Yi , we have

(28) ⇔ ∀i ∈ Vr∗ , ∃(Yi 0, d0 ∈ R): ai d0 + Yi bi . (29)

Then, the route is feasible if and only if there exists d0 ∈ R and Yi 0 such that for all i ∈ Vr∗ :

ai − Yi d0 bi − Yi . (30)

For all (i, j ) ∈ (Vr∗ )2 such that j < i, we need to show that there exists a sum of minimal non-negative waiting times
between stations j and i, and a departure time d0 ∈ R from the terminal such that

aj − Yj d0 bj − Yj and ai − Yi d0 bi − Yi , (31)

whenever aj bi .

We get

∃(Yi − Yj 0, d0 ∈ R):

[aj − Yj d0 bj − Yj and ai − Yi d0 bi − Yi ] (32)

⇔ ∃(Yi − Yj 0): [aj − Yj bi − Yi and ai − Yi bj − Yj ] (33)

⇔ ∃(Yi − Yj 0): ai − bj Yi − Yj bi − aj (34)

⇔ bi − aj 0. (35)

Thus bi max0  j <i {aj }. 

Proof of Proposition 2. We have shown in the proof of Proposition 1 that a route is feasible if and only if there
exists for all vertices i and j > i a non-negative sum of waiting times Yj − Yi between i and j within the interval
[aj − bi , bj − ai ] (Eq. (34)). Then there exists a non-negative sum of waiting times w = Yn which can be decomposed
as follows:
w = Yn (36)
= (Yj2 − Y0 ) + (Yj1 − Yj2 ) + (Yn − Yj1 ), (37)

where j1 = arg maxi∈Vr∗ {ai } and j2 = arg mini∈Vr∗ {bi }.

Since the route is feasible, for each i < j2 , there exists a non-negative value Yj2 − Yi such that aj 2 − bi Yj2 −
Yi bj 2 − ai , and we have bj 2 − ai 0. Since aj 2 − bj 2 0 is true by definition, Yj2 − Yi can always take a zero value
for each i < j2 and a fortiori for i = 0. On the other hand, for each i > j1 , there exists a non-negative value Yi − Yj1
such that ai − bj 1 Yi − Yj1 bi − aj 1 , and we have bi − aj 1 0. Since aj 1 − bj 1 0 is true by definition, Yi − Yj1
can always take a zero value for each i > j1 and a fortiori for i = n. Then, wr = Yj1 − Yj2 and, because the route is
feasible, there exists w 0 such that

aj 1 − bj 2 w (38)

⇔ max {ai } − min {bi } w (39)

0i n 0i n

⇔ max∗ {ai } − min∗ {bi } w (40)

i∈Vr i∈Vr
 F. Cornillier et al. / Computers & Operations Research 36 (2009) 919 – 935 935

and
wmin = max∗ {ai } − min∗ {bi } (41)
i∈Vr i∈Vr
= wr . (42)

When maxi∈Vr∗ {ai } − mini∈Vr∗ {bi } 0, there exists a feasible time departure from the terminal which is common to
all stations without waiting time. In this case, we have wr = 0. 

Proof of Proposition 3. Since waiting time always delays the return to the terminal even if departure occurs  time
units before r , it is always preferable to start from the terminal as late as possible, i.e. at d0 = r = mini∈Vr∗ {bi }.
Any other departure time r −  just makes the route duration  longer. Then, if waiting is needed, we have r =
maxi∈Vr∗ {ai } − wr = r . 

References

[1] Brown GG, Graves GW. Real-time dispatch of petroleum tank trucks. Management Science 1981;27:19–32.
[2] Brown GG, Ellis CJ, Graves GW, Ronen D. Real-time, wide area dispatch of mobil tank trucks. Interfaces 1987;17:107–20.
[3] Malépart V, Boctor FF, Renaud J, Labilois S. Nouvelles approaches pour l’approvisionnement des stations d’essence. Revue Française de
Gestion Industrielle 2003;22:15–31.
[4] Taqa allah D, Renaud J, Boctor FF. Le problème d’approvisionnement des stations d’essence. APII-JESA Journal Européen des Systèmes
Automatisés 2000;34:11–33.
[5] Cornillier F, Boctor FF, Laporte G, Renaud J. An exact algorithm for the petrol station replenishment problem. Journal of the Operational
Research Society, 2007, doi:10.1057/palgrave.jors.2602374.
[6] Avella P, Boccia M, Sforza A. Solving a fuel delivery problem by heuristic and exact approaches. European Journal of Operational Research
2004;152:170–9.
[7] Cornillier F, Boctor FF, Laporte G, Renaud J. A heuristic for the multiperiod petrol station replenishment problem. European Journal of
Operational Research, 2007, doi:10.1016/j.ejor.2007.08.016.
[8] Johnson DB. Finding all the elementary circuits of a directed graph. SIAM Journal on Computing 1975;4:77–84.
[9] Christofides N, Mingozzi A, Toth P. Loading problems. In: Toth P, Christofides N, Sandi C, editors. Combinatorial optimization. Chichester:
Wiley; 1979. p. 339–69.
[10] Yuceer U. A multi-product loading problem: a model and solution method. European Journal of Operational Research 1997;101:519–31.
[11] Smith JC. A genetic algorithm approach to solving a multiple inventory loading problem. International Journal of Industrial Engineering
2003;10:7–16.
[12] Bukchin Y, Sarin SC. Discrete and dynamic versus continuous and static loading policy for a multi-compartment vehicle. European Journal of
Operational Research 2006;174:1329–37.
[13] Fagerholt K, Christiansen M. A combined ship scheduling and allocation problem. Journal of the Operational Research Society 2000;51:
834–42.
[14] Fagerholt K, Christiansen M. A travelling salesman problem with allocation, time window and precedence constraints—an application to ship
scheduling. International Transactions in Operational Research 2000;7:231–44.
[15] Oppen J, LZkketangen A. A tabu search approach for the livestock collection problem. Computers & Operations Research, 2007,
doi:10.1016/j.cor.2007.02.021.
[16] Eglese RW, Mercer A, Sohrabi B. The grocery superstore vehicle scheduling problem. Journal of the Operational Research Society 2005;56:
902–11.
[17] Van der Bruggen L, Gruson R, Salomon M. Reconsidering the distribution structure of gasoline products for a large oil company. European
Journal of Operational Research 1995;81:460–73.
[18] Bausch DO, Brown GG, Ronen D. Scheduling short-term marine transport of bulk products. Maritime Policy and Management 1998;25:
335–48.
[19] Savelsbergh MWP. The vehicle routing problem with time windows: Minimizing route duration. ORSA Journal on Computing 1992;4:
146–54.
[20] Helsgaun K. An effective implementation of the Lin-Kernighan traveling salesman heuristic. Datalogiske Skrifter No. 81, Roskilde University,
1998.
[21] Helsgaun K. An effective implementation of the Lin-Kernighan traveling salesman heuristic. European Journal of Operational Research
2000;126:106–30.
[22] Toth P, Vigo D. The granular tabu search and its application to the vehicle routing problem. INFORMS Journal on Computing 2003;14:
333–46.
[23] Malépart V, Renaud J, Boctor FF, La distribution des produits pétroliers au Québec: État de la situation. Technical Report, Université du Québec;
1998.