Roland Jachimowski, Micha Kodawski
Warsaw University of Technology, Faculty of Transport
SIMULATED ANNEALING ALGORITHM FOR THE
MULTI-LEVEL VEHICLE ROUTING PROBLEM
The manuscript delivered: May 2013
Summary: The paper presents the possibility of using simulated annealing algorithm for solving the
transportation planning problem in the multi-depot, multi-level distribution systems. The transportation organization issue in the multilevel distribution system was characterized and an optimization task adequate to considered decision situation was formulated. Also the idea of the simulated annealing algorithm and the method for its usage do solve the problem considered in this article was introduced.
Keywords: simulated annealing, vehicle routing problem, multi-level distribution
1. INTRODUCTION
Effective goods transportation planning is one of the key ways to minimize the cost and also to increase the level of customer service. While the planning of direct transport with the full truckload from suppliers to customers does not seem to pose a problem of combinatorial nature, in the goods distribution, especially those being a general cargo from manufacturer to large number of customers, transportation planning becomes a very complicated and hard to solve issue. This is mainly because of the many aspects that must be taken into account at the distribution planning stage. These aspects can for example refer to: limited vehicle capacity and their availability, time windows for the pick-up and the delivery, asymmetric travel costs or multilevel structure of the distribution system where the intermediary performing functions connected with the goods (cargo) movement to final customer determine the so-called distribution level [6]. Taking the number of intermediaries acting in the distribution system into account we can distinguish the system of the direct distribution (0-level) and system of the multi-level distribution [9].
The organization of goods flow in the supply chain, especially in the multilevel distribution system is a problem, that is very difficult to solve. This is due to the lack of appropriate tools to support the transportation planning process especially in distribution systems where there is a multiple number of intermediaries (different types of warehouses or transshipment terminals). In most cases, only selected elements of the distribution
system are planned with the usage of computer tools. Unfortunately in most cases the transportation services for a specific region is planned on the basis of certain fixed patterns, or based on the experience of those who are responsible for this planning. This way of decision making about transportation service realization does not guarantee its optimality with the criteria of cost and time minimization.
Solving this kind of a problem needs using methods from the operational research area. In operational research the problem of goods distribution from the manufacturer to many customers is called the Vehicle Routing Problem (VRP). This problem is an extension of one of the oldest optimization problems on network called the Traveling Salesman Problem (TSP) [3]. Because the very large computational complexity of the vehicle routing problem, only heuristics algorithms are able to solve it for large, real life instances. One of these heuristics algorithms, which is known in literature to be successfully used for solving the vehicle routing problem is the so-called Simulated Annealing Algorithm.
Simulated annealing algorithm is based on the observation of the technological process of metal alloys cooling. Thus the particles of substance in the liquid state at high temperatures have very high energy, what cause them move quickly and chaotically.
With the gradual decrease of temperature the energy of particles decreases too. Particles begin to move more slowly and at the same time strive to obtain a position in which they reach a minimum of potential. If the liquid metal cooling process is slow, then the particles reach a position where they will have the minimum of potential and they will form a solid crystal structure. Otherwise, if the cooling is done quickly, the particles cannot make it to achieve this minimum of potential and in the solid state they will be set up chaotically. [10].
Particles striving to achieve the minimum of the potential has its similarity in searching for the minimal route in the Traveling Salesman Problem. This analogy can be seen if we refer to the actual location of the particle in the cooling metal to one of the possible feasible solutions of the traveling salesman problem. In the TSP case or its variation (VRP), the algorithm starts its work with drawing any feasible solution. The algorithm calculates the value of the criteria function for this solution. In TSP this any feasible solution can state any permutation of all points that have to be visited by the salesman. The criterion is the salesman total route length after visiting all points. In the next step the neighborhood of this solution is being searched. If the solution is a permutation of points, then its neighborhood can be expressed as a set of permutations, which are obtained by performing on the previous permutation a single operation such as: swap, insert or invert.
x Swap operation is made by exchanging locations of two selected elements.
x Insert operation is made by choosing any element of the permutation and putting it on another location.
x Invert operation is made by choosing any two elements and inverting the order of the elements located between those two chosen.
In the next step algorithm chooses one of the neighbors. For this neighbor location the value of criteria function is calculated and marked as F(S). If the value of the criteria function for this neighbor is greater than or equal to the initial permutation, then this value is set as a new solution. If this value is smaller than the previous permutation, (initial permutation in first iteration of the algorithm), then this value is set as a new solution with
the some probability. Thanks to this operation, the algorithm can temporarily accept the worse solution to have a possibility to escape from the possible local minimum of the criteria function. The probability of a worse solution accepting is calculated with the following equation:
/T
P
e
' (1)where:
{= F(S’) – F(S) (2) T – temperature regulates the process of worse solutions acceptance. This temperature is
decreasing successively with every algorithm’s iteration, and the higher temperature is the more likely is that the worse solution will be temporarily accepted by the algorithm.
The idea of the Simulated Annealing algorithm was presented by N. Metropolis in 1953 [8]. Successful results obtained for different kinds of routing problems with the usage of this algorithm were for the first time reported in [7][4].
2. PROBLEM FORMULATION
The main problem considered in this article is to determine optimal delivery plans for the multi-level distribution system with the usage of the simulated annealing algorithm. Therefore for the purpose of researches it was necessary to identify the elements of the multi-level distribution system, which are: suppliers, central transshipment terminal (as a first level of distribution), transshipment terminals of the second and the third level and final customers. A characteristic feature of the system is that customers can be serviced by all the transshipment terminals with no matter to distribution level, but terminals of the lower level can be supplied only by the terminals of the upper level. Transport to suppliers and customers is performed by the vehicles belonging to transshipment terminals. The suppliers cargo volume and customers demand are known. In addition every customer is characterized by the time interval in which he demands to be serviced. Criteria function minimizes the overall cost of distribution in whole system, and this cost is influenced by:
x vehicles route length, x drivers working time,
x penalty costs for premature or late customer service. Data necessary to solve the problem is presented as follows: D ={d: d=1, ..., d, ..., d’,…, D} – set of suppliers,
O ={o: o=1, ..., o, …, o’, ..., O} – set of customers, ST ={st: st=1, ..., st, ..., ST} – set of vehicles,
M ={m: m=1, ..., m, …, m’,..., M} – set of transshipment terminals, N = {n: n=1, ..., n, ..., N} – set of distribution levels,
SD = [sdd,d’] – matrix of distances between suppliers,
SMD = [smdm,d] – vector of distances between suppliers and central transshipment
terminal,
TD = [tdd,d’] – matrix of travel times between suppliers,
TMD = [tmdm,d] – vector of travel times between suppliers and central transshipment
terminal,
SO= [soo,o’] – matrix of distances between customers,
SMO= [smom,o] – vector of distances between customers and transshipment terminal,
TO= [too,o’] – matrix of travel times between customers,
TMO= [tmom,o] – vector of travel times between customers and transshipment terminal,
kst – unitary cost for vehicle usage by a unit of a distance, ktst – unitary cost of driver working time,
qdd – supplier’s cargo volume,
qoon – customer’s demand,
+ – loading/unloading time,
dm,m’ – distance between terminal of level n and terminal of level n+1,
tm,m’ – travel time between terminal of level n and terminal of level n+1.
Other data necessary for the task, occurring only in the constraints are: Qst – vehicle’s capacity,
Tst – allowed driver’s working time,
tkm – consolidation time at transshipment terminal,
(Am, B m) – transshipment terminal’s working time,
(ao, bo) – customer’s time window.
Matrixes of the decision variables X, XMD, Y, YMO and Z for such data formulation are defined as follows:
X=[xstd,d’],
where xstd,d’=1 if arc (d,d’) belong to the route of vehicle st,
XMD=[xmdstm,d],
where xmdstm,d =1 if arc (m,d) belong to the route of vehicle st,
Y=[yst,no,o’],
where yst,no,o’=1 if arc (o,o’) belong to the route of vehicle st at nth level of distribution,
YMO=[ymost,nm,o],
where ymost,nm,o =1 if arc (m,o) belong to the route of vehicle st at nth level of distribution,
Z=[zstm,m’],
zstm,m’=1 if vehicle st travels from terminal of nth level to terminal of level n+1.
More over there is a variable used only in constraints and defined as follows: Ton – moment of the arrival at the customer.
,' ' , ' ' ,' ' , ' , , , ' , , ( , , ) n n st st st st d d st ST d D d D m M st st st st d d d D d D m M st ST st d d d d D d D n d d m d m d d d m d m d k x sd xmd smd F X Y Z kt x td xmd tmd qd x W § § ·· ª º ¨ ¨ ¸¸ « © © ¹¹» « » « § § ··» ¨ ¸ « ¨ © ¹¸» « ¨ ¸» « ¨¨ ¨§ ·¸ ¸¸» « » ¬ © © ¹ ¹ ¼
¦ ¦
¦
¦
¦
¦
¦
¦
¦
¦
1 1 ,, , ' , ' , ' ', ' , , , ' ' , , ' ' , ' , , , ' , n n n n n st st st st m m m m m m m m N m M m M st ST m M m M st ST st st n st n st o o st ST o O o O m M n N st st n s o o st ST o O o o m o m o o o m o z d k z t kt k y so ymo smo kt y to ymo ª º « » ¬ ¼ § § ·· ¨ ¨ ¸¸ © © ¹¹
¦ ¦ ¦ ¦
¦ ¦ ¦
¦ ¦
¦
¦
¦
¦
¦
,
, , ' ' , min n t n st st n o o o o O m M o O o O m o tmo qo y W ª ª ºº « « »» « « »» o « « §¨ §¨ ·¸ §¨ ·¸·¸»» « « © ¹ »» ¬ ¬ ©
¦
©¦
¹¦
¦
¹¼¼ (3)Due to the limited paper length, instead of writing analytical constraints, authors presented only content version of these constraints:
x every supplier can be serviced only once, x every customer can be serviced only once,
x vehicle arriving at the service point has to also leave that point, x vehicle’s capacity cannot be exceeded,
x overall vehicle’s route duration cannot exceed drivers maximum working time,
x service of suppliers has to be completed early enough to allow consolidation and a service of customers assigned to the central terminal,
x customers’ service has to be completed before the end of working time of the terminal they are assigned to,
x vehicle can service the customer only during the customer’s time window,
x travel time between terminals of the two successive distribution levels cannot exceed drivers working time,
x cargo has to be delivered to the terminal of level n+1 during terminals working time and early enough to realize consolidation and service of that terminal’s customers.
3. PROBLEM SOLVING METHOD
Transportation plans in such distribution system are obtained by using a hybrid heuristic cluster first/route second [2]. In this point clustering aims to define subsets of suppliers and customers called clusters separated from the set of all service points. Clusters are built with the criteria of not exceeding vehicles capacity, drivers working time etc. Each cluster is operated by a single vehicle, and there is a single traveling salesman problem performed in every cluster.
To solve the considered problem the service points clustering procedure was developed. Suppliers are clustered in the same way as customer, except that every supplier’s time window is as long as terminals working time and doesn’t influence the suppliers service order. This procedure is a development of the algorithm described in [5], and is implemented in the following steps:
1. Define the list of customers assigned to the given transshipment terminal and sort it by increasing beginning of customers time windows.
2. Define the list of vehicles assigned to the given terminal.
3. Set the algorithm parameters referring to maximum route length between any two customers and the maximum waiting time for service in every cluster.
4. Choose the vehicle, which will be the cheapest in term of its capacity compared to its fuel usage.
5. Assign the first unassigned customer (from the sorted list of customers) to the first (cheapest) vehicle.
6. Assign the second (next) unassigned customer from the sorted list of customers. 7. Check the constraints imposed on assigning of the second (next) customer to the
cluster.
8. Assign next customers to the given vehicle and check constraints every time the new customer is assigned. Close the cluster.
9. Choose the next vehicle and come back to step Bd! Nie mo na odnale róda
odwoania..
10. Stop the algorithm when all the customers are assigned to clusters. The above algorithm has been presented in detail in [1].
After the suppliers and customers clustering, collecting routes from suppliers to central transshipment terminal are performed. In the central transshipment terminal goods are consolidated and prepared for their shipment to the central terminal’s customers and to the transshipment terminal of the lower level. Collecting routes from suppliers and delivery routes to customers are obtained with the usage of the simulated annealing algorithm.
Simulated annealing algorithm, due to its specificity requires the determination of the parameters relating to the initial temperature Tmax which is necessary to calculate the
probability of choosing worse solution and parameters relating to the final temperature Tmin
which is necessary to stop the algorithm. In addition to decrease the temperature in every next iteration of the algorithm it is necessary to define a decreasing function (which will reduce that temperature). In the algorithm prepared for the purpose of this article, this function g(t) in every next iteration i+1 has a value that is a product of the temperature from the previous iteration ti and the value of the cooling factor . Thus, the function g(t) is
defined as follows:
1
(
i)
ig t
O
t
(4) where:ti – temperature value in the i-th iteration,
Other algorithm parameters: Best – best found solution,
P(t) – probability function of choosing new solution (choosing worse solution),
( ') ( ) ( ) F S F S t P t e (5) where: S – previous solution, S’ – new solution.
The general simulated annealing program code is presented as follows:
t = T_MAX; A = initial_solution(); BEST = A; while(t > T_MIN) { i = 0; while(i < MAX_ITER) { B = new_solution(A);
if( (f(B) < f(A)) or (losuj(0,1) < P(t)) ) { A = B;
if( f(A) < f(BEST) ) BEST = A; } i++; } t = g(t); }
The computational complexity of the simulated annealing algorithm is dependent on the number of iterations of two loops:
x external, where the number of iterations depends on the initial temperature Tmax, final
temperature Tmin and a cooling factor 8,
x internal, depending on the given number of iteration marked as MAX_ITER. The number n of external loop iteration satisfy the following assumption:
Tmax · 8n | Tmin (6)
It gives a conclusion that this number will be proportional to ݈݃ ቀ்
்ೌೣቁ. Therefore, the total complexity of the algorithm will be Oቂܸ ȉ ݈݃ ቀ்்
ೌೣቁቃ.
The steps of the simulated annealing algorithm that was used to solve the problem considered in this article are organized as follows:
1. Set the initial temperature Tmax.
2. Set the final temperature Tmin.
3. Set the cooling factor 8.
4. Set the temperature in first iteration ti as Tmax and go to step 6.
5. Set the value of the temperature in the iteration i+1 with the usage of g(ti+1) function.
6. Generate a random permutation of all service points in a given cluster.
7. Calculate the value of F(S’) for the permutation generated in step 6 (permutation of service points – S’).
8. Generate the MAX_ITER number of neighborhood of the permutation S’ with the usage of swap operation.
9. Calculate the value F(S) for any single neighbor S. 10. Check the constraint: F(S) > F(S’)
If the constraint is satisfied, go to step 11.
If the constrain isn’t satisfied, accept S with probability expressed by a function P(t). After that if the S got accepted, go to step 11. In other case go to step 12.
11. Save F(S)as a new solution.
12. Check if the algorithm is going to reach Tmin in the next iteration (i+1).
13. If algorithm has reached Tmin, go to step Bd!Niemonaodnaleródaodwoania.. In
other case go back to step 6. 14. End.
Schematically the algorithm steps are presented in figure 1.
Fig. 1. A simplified block diagram of a simulated annealing algorithm
START Set initial temperature Tmax
Set cooling factor8
Set ti as Tmax in 1st iteration Set temperature in i+1 with
usage of g(ti+1) Generate random permutation
of service points S’ Calculate F(S’) Generate neighborhood of S’
with swap operation Calculate F(S) Check F(S) > F(S’) Constraint satisfied YES NO Accept S with probability P(t) Save F(S) as a new solution
Check if algorithm has reached
Tmin in i+1
YES END
Set final temperature Tmin
Check if S accepted Constraint satisfied YES Constraint satisfied NO NO
4. SUMMARY
The algorithm presented in this article is a result of work carried out on the performance of metaheuristic optimization algorithms used in the literature to solve issues relating to transport and distribution planning. This work was carried out under the grant 504M/1160/0412 “Testing algorithms performance for transportation planning in multilevel distribution system“. The quality and the efficiency of the algorithm, was tested on the example of 3-level distribution system located on the polish territory, consisting of 98 suppliers, 559 customers and 22 transshipment terminals. The algorithm was implemented in author’s computer application. In a result the tabular plan of transportation from suppliers to final customers with the usage of intermediate transshipment terminals was obtained. Moreover graphical version of the transportation plan was generated for every single vehicle servicing suppliers and customers. This make it possible to quickly verify the correctness of the obtained results. Thanks to the characteristics referring to calculation duration and characteristics referring to the value of criteria function in subsequent iterations of the algorithm, it was possible to test the quality and efficiency of the developed simulated annealing algorithm.
Due to the limited size of the article, the quite extensive calculation example (3 levels of distribution) which was used for algorithm testing was not presented. Nevertheless, the results of the work are available in the grant report which was completed in April 2013.
Acknowledge
This work is carried out under the grant 504M/1160/0412 „Testing algorithms performance for transportation planning in multilevel distribution system”.
Literature
1. Ambroziak T., Jachimowski R.: Algorytm klasteryzacji w zastosowaniu do problemu trasowania pojazdów. Logistyka 2: 337-342, 2012.
2. Basley J.: Route First-Cluster Second Methods for Vehicle Routing. Omega 11: 403-408, 1983.
3. Bektas T.: The multiple traveling salesman problem: an overview of formulations and solution procedures. Omega 34: 209 - 219, 2006.
4. Cerny V.: A thermodynamical approach to the travelling salesman problem: an efcient simulation algorithm. Journal of Optimization Theory and Applications 45: 41-55, 1985.
5. Dondo R., Cerda J.: A cluster based optimization approach for the multi-depot heterogeneous fleet vehicle routing problem with time windows. European Journal of Operational Research 176: 1478-1507, 2007.
6. Jacyna M.: Modelowanie i ocena systemów transportowych. Oficyna Wydawnicza PW, 2009.
7. Kirkpatrick S., Gellat C.D., Vecchi, M.P.: Optimization by simulated annealing. Science 220: 671-680, 1983.
8. Metropolis N., Rosenbluth A.W., Rosenbluth M.N., Teller A.H., Teller, E.: Equation of state calculation by fast computing machines. Journal of Chemical Physics 21: 1087-1091, 1953.
9. Tsiakis P., Shah N., Pantelides C.C.: Design of multi-echelon supply chain networks under demand uncertainty. Industrial and Engineering Chemistry Research 40: 3585–3604, 2001.
10. Woch M., ebkowski P.: Sequential Simulated Annealing for the Vehicle Routing Problem with Time Windows. Decision Making in Manufacturing and Services 3: 87–100, 2009.
ALGORYTM SYMULOWANEGO WYARZANIA DLA PROBLEMU UKADANIA TRAS POJAZDÓW W WIELOSZCZEBLOWYM SYSTEMIE DYSTRYBUCJI
Streszczenie: W artykule przedstawiono moliwoci zastosowania algorytmu symulowanego wyarzania do
rozwizywania problemu planowania transportu w wielobazowych, wieloszczeblowych systemach dystrybucji. W tym celu scharakteryzowano problematyk planowania przewozów w systemach dystrybucji o strukturze wieloszczeblowej oraz sformuowano zadanie optymalizacyjne adekwatne do rozwaanej sytuacji decyzyjnej. W artykule przedstawiono take ide algorytmu symulowanego wyarzania oraz metod jego wykorzystania do rozwizywania zdefiniowanego problemu badawczego.
