Location-Routing for Distribution Centers

Read this article. One objective of this paper is to determine distribution center locations. Compare and contrast the two cases presented.

Solving approach

To validate mathematical model, two different numerical cases are solved by optimization software GAMS. Also, due to linear nature of the model, CPLEX solver was further used. The numerical cases were run under different scenarios using a personal laptop equipped with an Intel® Core™ i3 CPU at 2.4 GHz and 6 GB of RAM.

A location-routing problem combines location and routing problems which both of them are NP-hard class problems, so LRP is NP-hard as well. To solve a NP-hard problem, a proper solving algorithm and method is required because as the size of the problem increases, processing time increases exponentially, as is obviously seen in the second numerical case.

In this study, a genetic algorithm with a new chromosome structure is presented to solve problems with different sizes. The proposed algorithm is composed of the following components.


Chromosome presentation

Proposed chromosome consists of two separate parts. The first part refers to the multimodal routes from supplier to DCs (on multimodal network), and the second part refers to DCs location and tours routing from DCs to retailers (location-routing problem).

(a) In the first part, a matrix is developed to represent the multimodal route. The matrix dimensions are determined by the number of potential DC and intermediate nodes (between supplier and potential DCs). Number of potential DCs is taken as the number of matrix rows and each row represents a multimodal route for the corresponding potential DC. Regarding the number of columns, for each intermediate node two columns added to the matrix, i.e., each node is represented by two elements in each row (multimodal route). The first element has a binary (0 or 1) value showing if the intended node is used in the route or not. The second element distinguishes the transportation mode via which the product travels from the current node to the next one (one for road transportation, two for rail transportation, and three for sea transportation). Also, an extra first column is inserted to show the transportation mode via which products travel from the supplier to the first intermediate node along the route. Here the subject is clarified by a numerical example. Suppose a multimodal network with one supplier, six intermediate nodes, and four potential DCs. The following matrix expresses typical routes from the supplier to DCs.

\left[ {\begin{array}{*{20}c} 1 & 1 & 1 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 1 & 2 \\ 2 & 0 & 0 & 1 & 2 & 0 & 0 & 1 & 2 & 1 & 2 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 & 3 & 1 & 3 & 0 & 0 & 0 & 0 \\ 2 & 1 & 1 & 0 & 0 & 1 & 2 & 0 & 0 & 0 & 0 & 1 & 1 \\ \end{array} } \right].

The first row represents a route from the supplier to potential DC number one (DC 1). The first element along the row is 1, that is, the transportation starts from the supplier on road mode. The second element along the row is 1, which indicates the first intermediate node along the route. The third element shows that the transportation proceeds from the node 1 via road transportation. Being zero, the next pair of elements reveals that the intermediate node 2 is not used in the route. The sixth and seventh elements along the row show that the route enters the intermediate node 3 and leaves the rout via road transportation. And finally the last two elements determine that the intermediate node 6 presents along the route, and the route continues towards potential DC 1 by rail transportation with a mode change happening at node 6. Applying this process to the next DCs, row 2 shows that route to DC number two (DC 2) is on rail network and crosses the intermediate nodes 2, 4 and 5. Furthermore, row 3 presents the route to DC number three (DC 3). Nodes 3 and 4 uses in this route and transportation mode change at node 3 from road to seaway. Lastly, row 4 corresponds to the route to potential DC number four (DC 4); it indicates that intermediate nodes 1, 3 and 6 are used in the route and mode changing happened from railway to road on node 1 and 6 and from road to railway on node 3.

(b) The second part of the chromosome is a string of numbers representing established DCs and sequence of retailers along routing tours. The string length is equal to total number of retailers plus potential DCs minus one. To demonstrate the string, a sample string with nine retailers and four potential DCs is presented as follows:

1 \; 2 \; 4 \; 12 \; 11 \; 8 \; 9 \; 7 \; 6 \; 10 \; 3 \; 5.

The numbers 1–9 present 9 retailers, and the numbers 10, 11, and 12 refer to the DCs. Retailers are assigned to DCs via the following process. From the start of the string to one of the numbers 10, 11, or 12 (the one appeared first) will be assigned to DC 1. Obviously, if the string was started with one of the corresponding numbers to DCs, no retailer would be assigned to DC 1. In this case, retailers 1, 2, and 4 will be assigned to DC 1, with the same sequence. The retailers falling between the first and second corresponding numbers to DCs (10, 11 and 12) will be assigned to DC 2. In this example, as number 11 comes right after number 12, no retailer is assigned to DC 2, i.e., DC2 is not going to be established. Similarly, the retailers falling between the second and third corresponding numbers to DCs (10, 11 and 12) will be assigned to DC 3. In this case, retailers 8, 9, 7, and 6 are assigned to DC 3, with the same sequence of meeting the retailers. And lastly, retailers falling within the third corresponding number to DCs (10, 11 and 12) to the end of the string (retailers 3 and 5, respectively, in this example) are assigned to DC 4.

Generation of initial population

Two different populations with equal numbers should be generated for the two different parts of the chromosome. Each member of the matrix part of the chromosome shows a possible route from supplier to potential DCs. To generate an initial population for this part, the algorithm is fed by all inbound links into intermediate nodes along with their types. The algorithm starts from potential DC 1 and selects a random inbound link from the set of corresponding inbound links. Then starting node of the selected link is determined before a random inbound link is assigned to the nodes. This process continues until the node 0 (supplier) is reached. The obtained route is then transformed into the corresponding chromosome representation as described above. Repeating the process for other potential DCs, matrix rows and columns are formed.

The initial population for the string part of chromosome is created by a routine called "randperm" in MATLAB. For each matrix part of the chromosome, a string part is generated and their costs are summed up into a single response.


Parents selection mechanism

Parent selection mechanism is an important part of the genetic algorithm. In the present research, roulette wheel selection mechanism was used.


Crossover operator

In crossover operation, for the matrix part, in selected parents rows are replaced with each other, as follows:



For the string part of the chromosome, a random number is selected between 2 and length of string minus 1. This number divides selected parents into two parts. Combination of these parts generates two offspring. First part of first parent along with second part of second parent make offspring one and second part of first parent along with first part of second parent make offspring two. Following with the operation, modifications are done and repetitive members replace missing ones. An example is presented below to clarify the operation.


Note that, for both crossover operations, the operator works only in 50% of the time. So, under this condition, offspring with changes in just one part of the chromosome are possible.


Mutation operator

In the matrix part of the chromosome, mutation operator performs as follows. Based on mutation rate, number of rows is selected randomly and for each selected row, a new route generates using of same process in initiating the population. In the next step, selected and primary routes (rows) distances are computed and compared together. If the new route has less distance, it replaces in the related row, otherwise primary row reserves.

For the string part of the chromosome, mutation operator is randomly selected from a set of five different mutation functions. In the first function, two different members of the string are randomly selected to be replaced with each other.


In second function, in addition to replacing the members, the order of members between them is reversed.


In third function, two members are random taken and the first member is moved after the second member.


In the fourth function, the assigned retailers for each tour are determined. Then, their order is reversed in the routes.


In the fifth function, after determining the retailer tours, two tours are selected randomly and a random member of each selected tour is exchanged.


Similar to crossover operator, mutation operators apply at a probability of 50%. As mentioned, this allows for the generation of offspring wherein only one part of the chromosome is changed.


Fitness function

The model's objective function is used as the fitness function. For each matrix part of the chromosome, a string part is generated; coupled together, both parts are used as a single chromosome in fitness function. For each chromosome, four cost functions and one penalty function are developed. The first cost function addresses the cost of routing tours from DCs to retailers. As explained before, open DCs and their allocated tours are determined based on the string part of the chromosome. Then, DC numbers are added to the beginning and end of each related tour before determining the distance between each pair of nodes along the tour. Lastly, overall tour cost can obtain by summing up the calculated distances and multiplying the result by the product's unit transportation cost.

In second cost function, established DCs are distinguished using string part of the chromosome and summation of their establishment costs is calculated.

To determine transportation cost on the multimodal network (from supplier to potential DCs), the following process is applied. As a first step, the matrix part of the chromosome is decomposed into its rows and in each row the intermediate nodes used along the route are identified. In the second step, transportation mode between nodes is distinguished. In third step, distances between nodes are determined and multiplied by the corresponding transportation cost; sum of the results gives total cost for each link. Finally, total demand for each DC is calculated and multiplied by link costs. Sum of these costs form the multimodal network transportation cost. Last cost function is the cost of providing mode changing facilities at multimodal terminals. To have the cost calculated, nodes which mode changes are happened should be identified. For this mean, inbound and outbound modes are determined for each node along multimodal links. If the modes are identical, then node may host no mode change; otherwise mode changing facilities need to be established on that location. Fixed cost of providing mode changing facilities is aggregated to form the fourth cost function.

For this fitness function, a penalty is considered for violating vehicle capacity along the routing tours. If sum of assigned retailers' demands to a potential DC violated corresponding vehicle capacity, the penalty is applied.


Forming the next population

To form the next generation, first population along with offsprings generated by crossover operator and mutants are sorted on the basis of their costs. New population with the same size as first population is selected from the first members of the set (fewer costs are selected).


Stop condition

The algorithm is set to stop once a predefined maximum number of iterations is achieved.


Algorithm parameter setting

Parameter setting is an important part of coding the algorithm where different parameter configurations contribute to the quality of answers. In this study, running the algorithm with different parameter configuration arrived to following optimum values of GA parameters (Table 1).

Table 1 Genetic algorithm parameters

Parameter name Value Parameter name Value
Population size 400 Roulette wheel selection pressure 20
Mutation rate 0.8 Mutation implementation rate 0.8
Crossover rate 0.8