An Adaptive Genetic Algorithm for a New Variant of the Gas Cylinders Open Split Delivery and Pickup with Two-dimensional Loading Constraints

This paper studies a combination of two well-known problems in distribution logistics, which are the truck loading problem and the vehicle routing problem. In our context, a customer daily demand exceeds the truck capacity. As a result, the demand has to be split into several routes. In addition, it is required to assign customers to depots, which means that each customer is visited just once by any truck in the fleet. Moreover, we take into consideration a customer time windows. The studied problem can be defined as a Multi-depots open split delivery and pickup vehicle routing problem with two-dimensional loading constraints and time windows (2L-MD-OSPDTW). A mathematical formulation of the problem is proposed as a mixed-integer linear programming model. Then, a set of four class instances is used in a way that reflects the real-life case study. Furthermore, a genetic algorithm is proposed to solve a large scale dataset. Finally, preliminary results are reported and show that the MILP performs very well for small test instances while the genetic algorithm can be efficiently used to solve the problem for a widereaching test instances. Keywords—Vehicle routing problem; split delivery and pickup; multi-depot; two-dimensional loading; genetic algorithm


I. INTRODUCTION
LPG logistics transportation, in particular the distribution of the Liquefied Petroleum Gas (LPG) cylinders is considered among the basic building blocks in the Oil and Gas downstream supply chain and also known to be a very complex supply chain [11]. A typical LPG downstream supply chain consists of filling plants, distribution locations, fleet of trucks and customer's depots. The filling plant is an industrial unit composed of several processes in order to fill a broad range of gas cylinders. It should also be noted that the production and storage capacity is different from one filling plant to another. Therefore, each filling plant has an attached distribution location which is a cylinder's storage unit supplied from the production unit, then serves a set of customers by using a fleet of trucks. Moreover, due to the flammable nature of the gas, it is usually stored in liquid form under a specific pressure, in steel or composite plastic gas cylinders carefully checked and protected against deterioration. In general, the companies do not sell the gas cylinders but just the gas contents inside and the packaging remains their property. This type of packaging management system is called consignment.
Motivated by a case of a Moroccan petroleum company, our paper considers a real-life application in the LPG distribution industry. The objective is to deal with a problem which corresponds to current industry practice and also to give solution methods able to solve real large-scale instances in very fast computation times. The current distribution policy of our case study links two main entities: the gas filling plants and customers depots. Furthermore, each customer depot is served fully by only one truck. So, deliveries are very difficult to manage and generates a huge cost related to the production capacity and the use of trucks. Thus, the GPL Company faces a real problem which is serving customers in a way to optimize both loading and vehicle routing.
Through this document, we will propose a new distribution policy which will remove the dependency constraint between customer depots and filling plants. Those customers can be served by a fleet of heterogeneous trucks based in multiple filling plants. We note that all customers are visited just once by the truck in the same route. Every truck visits a number of customers in a certain order along its route. Thus, the new approach requires the assignment of customers to filling plants. In addition, a fleet of trucks is based at each filling plant, then, each truck originates from one filling plant serves the customers assigned to that plant. The objective of this policy is therefore to offer a global distribution model at minimal cost. This model concerns both the operation of loading cylinder racks into each truck and the construction of the trucks routes.
The remainder of this paper is organized as follows: Section 2 presents the loading and distribution problem considered in our case study. In Section 3, we explore a related literature review. In Section 4, we present a mathematical formulation of our problem based on a mixed integer programming model. Then, we present and discuss computational results. Section 5 presents a developed genetic algorithm (GA). Finally we conclude this paper by discussing results obtained by the GA and compare it with Cplex results, pointing out a few remarks of practical relevance and giving some perspectives of our further researches.

II. PROBLEM DESCRIPTION
The distribution problem studied in this paper includes two sub-problems: a two-dimensional loading vehicle problem and an open multi-depot split delivery and pickup routing problem with time windows (MD-OSPDTW).
As customer requests are in the form of a set of twodimensional and rectangular racks, the goal is to find an optimal combination to load it into a set of dissimilar rectangular trailers of height H and width W, without overlapping. The two-loading vehicle problem can be described as a loading of a limited fleet of heterogeneous trucks of capacity P k and a loading surface with length L k and height H k . During loading, racks have a prefixed orientation with respect to the trailers. In addition, all of them should occupy the entire surface width of the truck. We note that in our studied problem three types of trucks are considered.
We consider N racks of several types of gas cylinders to be shipped to r customers. Each rack c is characterized by a width w c <L k , a height h c <H k and a weight q c known. A problem instance includes a fleet of heterogeneous trucks, initially distributed among filling plants and may differ in term of their cylinders loading capacities. Furthermore, each truck can pull two types of trailers: flatbed trailer and lowboy trailer. Since, the order of each customer is a combination of different types of racks in which we put specific types of gas cylinder; we assume that the loading racks combinations of each truck are known. In addition, the loading surface in each truck can be represented on a rectangular coordinate system (x,y) with origin (0,0), in which the height represents the y-axis and the length represents the x-axis. To summarize, the truck loading is feasible if and only if all racks are placed completely inside the tray surface, two different types of racks cannot overlap and all racks are placed vertically parallel to the truck head.

B. The MD-OSPDTW Problem
The problem considered in our paper can be defined on a directed graph G = (V,A) in which V = V 1 ,. . . V n ,. . . V n+r represent the vertices and A = (V i ,V j ) ; i =j represent the arcs. The set of arcs comprises: V c = V 1 ,. . . V n as vertices corresponding to the filling plants and V e = V n+1 ,. . . V n+r as vertices representing customer's depots. We indicate by c ij the transportation cost between two vertices i and j. Moreover, customers have a daily requirement D i of cylinders, a service time S i and an associated time window. We have to take into account the trucks service time in filling plants s = (t + w) which represent both, the duration t of the trucks processing and the loading or unloading time w. Furthermore, the requirement of each customer may exceed the capacity of the truck, so it has to be split into several routes and each customer can be visited several times by different routes. In addition, each depot has a specific time window and should not be visited outside of this interval time. The objective is to determine for each truck the loading combinations and the set of routes to achieve in order to minimize the total transportation costs.
The problem in our paper is illustrated in Figure 1 and can be considered as an open vehicle routing problem (OVRP). It is a generalization of VRP in which a truck may not return to its departure depot after shipping the last customer on a route. Moreover, instead of the single-depot, deliveries are made from several multi-depots. The considered transportation problems include also more constraints: The time window restriction (OVRPTW), the multi-depots constraints, the fleet of trucks is assumed to be heterogeneous (HOVRP); the consideration of multi-types of products (MPOVRP). In this paper the transportation problem as a whole can be named as: open split delivery and pickup routing problem with multi-depot, multi-Product, time windows, using a heterogeneous fleet of trucks and obviously with consideration of Two-dimensional loading constraints. To the best of our knowledge, the OVRP as we present it has not been considered in the literature.

III. LITERATURE REVIEW
The vehicle routing problem (VRP) is the most studied combinatorial optimization problem in transport and logistics. It was firstly presented in the literature by [1] as the truck dispatching problem. It can be defined as the determination of an optimal set of routes for a fleet of vehicles which needs to serve a set of customers. Few years later, a considerable number of variants have been considered by researchers: Capacitated vehicles, multi-depots, time-windows, pickup and delivery, heterogeneous fleet, etc. In view of the problem statement, the main purpose of this paper is to focus on the combination of loading constraints and vehicle routing problems. Both problems belong to the NP-hard type of optimization problems. Thus, the combination of these two NP-hard problems gives us a more complicated problem.
In this section, we review only the closely related works which deal with the multi-depots open split pickup and delivery routing problems with time windows and with twodimensional loading constraints.
The variant of OVRP was firstly introduced in 1981 by [3]. Since then, it has gained much attention from researchers mostly in the studies addressing companies that avoiding being charged for the return trip of the vehicles to their departure depots. So, the important feature of the OVRP problem is that the vehicles are not required to return to the starting depot, but if they do, it must be through visiting customers in reverse order [10]. Authors in [6] presented the first heuristic algorithm to solve the OVRP. They developed a method based on a minimum spanning tree with penalties procedure. Afterwards, a tabu search algorithm was presented in [8] in which authors explore the structure of the OVRP problem and compare its performance with results obtained by other heuristics methods. One year later, the same heuristic was proposed by [10] for finding the routes that minimize both total travelling cost and operating cost while satisfying vehicle capacity and route length constraints. Many researchers have developed various heuristic approaches to efficiently solve the OVRP, such as ant colony optimisation by [19], simulated annealing by the researchers in [9] and genetic and evolutionary computing in [20].
To address the real-world business issues, the OVRP can be combined with more complicated constraints, such as split delivery (SDVRP). The SDVRP was introduced in the first time in [5]. The Authors studied this problem to prove that it can generate savings by allowing divisible deliveries. In [21] a study of the maximum possible savings achieved by enabling divisible deliveries is presented. Authors analyse the SDVRP algorithmic complexity with small capacity vehicles. Another work in which two hybrid genetic algorithms are used to solve a split delivery problem with respect to the total travel www.ijacsa.thesai.org distance is presented by [22]. In regards to the SDOVRP, we find only one interesting paper [30]. Authors developed a mathematical model of the OVRPSD that minimize both total fixed cost of vehicles and delivery workers and total transportation cost. Other features are also considered such as delivery and pickup (DPRP). This problem has received continual attention in transportation logistics field. The early work on the VRPPD was in [2], it was conducted for diala-ride scenarios. Later, researchers in [7] provided a lower limit where a demand of each customer does not exceed the vehicle capacity. In the VRPSPD, the vehicles are not only required to deliver goods to customer, but also pickup goods at customer locations. A general assumption in the VRPSPD is that all delivered goods must be originated from the depot and all pickup goods must be transported back to the depot. Very recently an iterated local search metaheuristic and a branchand-price algorithm are proposed to solve the multi-vehicle one-to-one pickup and delivery problem with split loads [27]. Moreover, Authors in [31] studied a barge transportation of maritime containers between a dry port and a set of sea port terminals. They develop a hybrid local search metaheuristic algorithm, combined with a branch-and-cut method to solve the VRPSPD problem. The objective is to find the best allocation of containers to barges in order to guarantee on-time delivery and meet capacity restrictions. In a previous paper we have introduced the split delivery and pickup vehicle routing problem with two-dimensional loading constraints. We attempt to give solutions for small instances in order to validate the model [26].
Customer's time window is also an important characteristic of the studied OVRP. The modeling of such problems requires consideration of the arrival time of each customer, which requires the integration of time windows constraints of each customer. The first paper introduced the open vehicle routing problem (OVRPTW) was studied by [15] Authors studied a real-world cases such as the delivery of school meals, school bus routing, the plans of passing through tunnels of trains, etc. The problem variants are extended to include the case where more than one depot is considered. The multiple depots open routing problem (MD-OVRP) is a variant of the capacitated vehicle routing problem CVRP in which a fleet of vehicles is based at each depot and customers are assigned to depots.
Recently, [29] proposed a two-phase algorithm to solve a lowcarbon multi-depot open vehicle routing problem with time windows (MDOVRPTW) with minimum total costs. In the literature, several papers have studied OVRP and VRPSDP with times windows and multi-depots constraints [8] [28] [29].
Another complicated problem presented in this paper which is the two-dimensional loading constraint. In recent years, a problem of loading customer items into vehicles has become a fertile ground for researchers in the field of combinatorial optimization. Nevertheless, the combination of the loading constraints and vehicle routing problems has received a little attention in the literature. The first work related to transportation problem with loading constraints was introduced by [4]. They address the Pickup and Delivery Traveling Salesman with LIFO Loading where a single rear-loaded vehicle serves a set of customer orders with LIFO policy. Next, [12] addressed the Two dimensional Capacitated Vehicle Routing Problem with Loading Constraints (2L-CVRP) using an Ant Colony Optimization (ACO). In another work they present an exact approach based on branch-and-cut for solving the same problem [14]. Researchers in [17] resolved a two dimensional Capacitated Vehicle Routing Problem using a tabu search algorithm. The 2L-VRP is combined to various VRP variants, such as time windows [16] [23][ [24]], Pickup and delivery [13] [18] and heterogeneous fleet [25].

IV. MATHEMATICAL FORMULATIONS
Two virtual plants are considered to formulate our distribution problem: a departure plant and an arrival one. We suppose that the transportation cost is null between these plants and the other considered plants. We define our problem as an oriented graph G = (V,A) with V=V 0 ,V 1 ,..V n+1 ,..,V n+r ,V n+r+1 ,. . . ,V 2n+r ,V n+r+1 and A = (V i ,V j ) ; i =j. The arcs A comprises: a virtual starting plant V 0 , a filling departure plants Vc =V 1 ,..V n+1 , a set of customers to be visited V e =V n+1 ,..,V n+r , the arrival filling plants V c =V n+r+1 ,. . . ,V 2n+r , and the arrival virtual plant V 2n+r+1 ( Figure 2). A. Assumptions.
• A truck cannot leave or arrive at customer depots outside a specific time window; • Each filling plan has a determined service time corresponding to unload racks of the empty cylinders and to load racks of the filled cylinders; • Each customer depot has a determined service time corresponding to unload racks of the filled cylinders and to load racks of the empty cylinders; • The travel time and the distance between the gas filling plants and the customer depots are well known; • At the beginning, all trucks are empty and supposed parked at the plants where it belongs; • The starting location of each tour which is the filing plant can never be a pickup location; • The last pickup location of the tour is the last served customer; • Each customer may be served one or more times; • We suppose the use of heterogeneous and limited fleet of trucks; • The racks should occupy the entire surface width of the truck; • The capacity of each truck is well known.

B. Parameters.
In order to formulate the above problem and solve it, we define the complete nomenclature and then state the model. The nomenclature is divided into parameters, sets and variables: • o : Virtual starting plant; • n : Number of Plants; • r : Number of customers; • 2n+r+1 : Virtual arrival plant; • K : Number of trucks in the fleet; • C ijk : Truck k transport cost between vertex i and j; • P k : Truck capacity; • H k : Maximum height of truck k; • L k : Maximum length of truck k; • C : Number of racks types; • Max ck : Maximum number of racks with type c on truck k: • D ic : Customer requirement of cylinders of the type of rack c.
• D i : Customer requirement of cylinders: • n c : Max number of cylinders in rack c; • h c : Rack c heght; • l c : Rack c lenght; • q c : Rack c weight when it is full; • t ij : Travel time between vertex i and j; • d ij : Distance time between vertex i and j; • M : Large number; • (e i , l i ) : Time windows Five decision variables are used in our mathematical model as follow: The arrival time of truck k at customer i.
b cik : Number of racks c delivered to customer i by truck k.
z ck : Number of piles c occupied in the truck k.
D. Objective and Constraints.
Using the above notation, the problem can be presented as follows.

5)
n+r j=n+1 The objective function (1) minimizes the total transportation cost. Constraints (2) impose the continuity of a tour; x ijk = 1 indicates that the truck k visits j after i and y ik = 1 indicates that the truck k visits i. So if a truck k visits a vertex it should necessarily leave it except for the vertex 2n + r + 1. Constraints (3) assure that the virtual departure plant 'O' can never be a destination, and the departure filling plant can only be a destination for this virtual plant. Constraints (4) enforce the virtual arrival plant 2n+r+1 to never be a provenance and the arrival filling plant to be a provenance only for the virtual arrival plant 2n+r+1. Constraints (5) ensure the passage by a departure filling plant and an arrival filling plant. Constraints (6) require the truck to leave the virtual departure plant if it is used. Constraints (7)-(8) enforce the truck to finish its tour at the virtual arrival plant and pass necessarily by the virtual departure plant. Constraints (9) assure that the quantity b cik delivered to customer i cannot exceed the requirement of this customer. Constraints (10) limit the overall load to not exceed the truck capacity. Constraints (11) impose the satisfaction of customer's demands. Constraints (12) determine the arrival time and represent subtour elimination constraints. Constraints (13)- (14) enforce the respect of time windows. Constraints (15) ensure that the total surface occupied by racks on each truck is less than the total surface of the tray space. Constraints (16) ensure that the sum of the widths of racks loaded in the tray is less than the length of the truck. Constraint (17) determines the number of piles occupied by each type c of racks. Constraint (18) indicates that each truck can only be started at initial parked filling plant.

V. RESULTS
In this section we start by describing the problem instances used within the computational experiments. The test instances considered in our paper are derived from a real-life dataset provided by a major petroleum company in Morocco. For the exact resolution, the dataset consists of 4 problem classes; each problem class is divided into 3 dispersed configurations of different sizes and the model was tested on all 12 instances. In our problem settings we have at most fifteen customer's depots C = C1...C15 supplied from four gas filling plants D = D1...D4 and using at most forty trucks of three different configurations.
Exact solutions are obtained by implementing the model using the commercial solver IBM ILOG CPLEX 12.5 and tested on a workstation with Intel(R) Core(TM) i7 CPU @ 2.60Ghz (8 CPUs) and 16 GB RAM. Our goal is to determine how large an instance could be solved to optimality in a reasonable amount of time.
Furthermore, the requirement for each customer for three different types of gas cylinders (B03, B06 and B12), the transport costs between each vertex and the dimensions of the storage racks are all known Furthermore, in the latest test class CL4 in which we have used the real volume of customers' demand, we use four filling plants and limiting the number of customers to be fifteen and the number of vehicles to be forty. Having discussed the savings generated by recovering the packaging, racks should return charged by empty cylinders. Thus, let us now consider the following decision rules: The delivered depot of order D i represent the pickup location of the same number and types of racks. In what follows, we detail our findings for each input setting. Table I reports also the result for each instance. It shows the instances characteristics, the obtained objective value, the best bound, the real gap percentage and the CPU time. We observe that for instances of small size (CL1-01, CL1-02 and CL1-03) the resolution does not take too much time to obtain an optimal solution at less than the 107th iteration. Furthermore, when we increase the number of customers to be served, CPLEX reaches good feasible solutions for medium size instances (CL2-01 to CL3-01). After ten tests, we tried to solve problem instances (CL4-02 and CL4-03); the execution took more than one hour without obtaining results. So, the solver can be stopped sooner without significantly finding a feasible solution. Moreover, we observe that both, the number of iterations and the time of resolution increase extensively with the number of customers affected by the volume of requirement.
Summing up, the MILP model solves the small-size problems (1 gas filling plant and 1 customer) optimally within seconds. When the number of customers increases to 7 (instance CL3-02), MILP could not prove the optimality of solutions within time limits. We note also that the customers' demands used in our first test instances are in the order of a hundred cylinders only, while the actual demand varies between 2000 and 11000. By testing the model of distribution with actual customers' demands (instance CL4-02), the solver ran for a long time without generating any results. Solving our problem requires examining a very large number of combinations. If there are few customers to visit, the number of combinations to explore is low and the problem is resolved quickly. However, adding only a few clients can dramatically increase the number of combinations to explore so that the resolution time becomes excessively long; this phenomenon is called the combinatorial explosion.
An example of an optimal solution to the routing, as well as to the loading problem with five customers and two gas filling stations is depicted in Figure 3. It represents the test instance CL2-3 in which six vehicles spread over the two filling stations. As we have considered in the problem formulation two virtual stations (0 and 0') with transport cost between these virtual stations and the other stations is equal to zero. In the MILP formulation we consider that all vehicles are initially located at depot O. In Table II we report the visited sequence of each tour and the related loading characteristics. As customers' order in the test instance CL2-03 are small sizes in which total demand for each type of cylinder do not exceed 1000, the distribution of orders is balanced between all vehicles. Thus, we get the optimality for this test instance.
The addressed model is considered to be a NP-hard combinatorial problem. Thereby, the use of exact optimization method to solve this problem may be difficult in an acceptable CPU times, mostly when the problem involves large data sets. Therefore, the exact method is used only for small size problems.

VI. GENETIC ALGORITHM FOR THE DISTRIBUTION PROBLEM
By drawing inspiration from the mechanisms of the natural evolution of living beings, many optimization problems are resolved successfully by using genetic algorithms. The GA is a metaheuristic, introduced for the first time by [32]. The main concern is to find approximate solutions for these problems in a reasonable computation time. To implement our GA, it is necessary to have a solution coding, parameters, a fitness function to evaluate the solutions and a mechanism to obtain the initial population. In addition, we need a method of selecting individuals to reproduce and genetic operators adapted to the problem.
In this section we detail our solution approach based on a genetic algorithm (GA) for the resolution of the gas cylinder distribution problem. We present the adopted coding, the initialization as well as the different genetic operators.

A. Solution Encoding
In our genetic algorithm we represent a coding based on a permutation representation for the adopted solution of the studied distribution problem. Each solution is represented in the form of a one-dimensional table in which a sequence represents the departure plant, the order of the visited customer's depots and the arrival plant. However, this representation poses a problem of delimitation between the filling plants and the customers. Therefore, we use a number n-2 of permutations in which n represents the number of customers in the tour, then we assume that the filling plants start from c = a +1 with n ≤ a (in our study we take a = 1000). An example of this encoding for the vehicle routing problem is given in [33]. Authors propose a permutation which contains both customers and route separators with a set [0,. . . , 9] represents the customers while [10,. . . ,12] are the route separators. Figure 4 represents an example of routes made by 3 trucks (K1, K2 and K3) in which a sequence of visited customers is assigned to them. In the first route made by truck K1, the departure plant represented by 1001, visited customers 4-1-6-2-3-5 (in that order), and the arrival plant as 1001. In the second route, truck K2 starts at same filling plant 1001 to serve customers 1-6-2-4 but this time it terminates at another filling center 1002. In the third route, truck 3 starts its route at plant 1002, visits customers 2-5 -3-1-6 and returns to the same departure plant. Furthermore, this representation makes it possible to represent trucks to which no route is assigned, when the number of trucks in the problem is fixed and greater than 1. This case is represented by the last empty route of the truck K4.

B. Initial Population Construction
To generate an initial population of our distribution problem we decided to make a mix between two techniques, by generating in a completely random way a half of the population and the rest with a priori good solutions. To build a solution using the random method, we assign randomly selected customers to each truck until all trucks reach their maximum loading capacity. It should be noted that this assignment is made while respecting the various constraints of the problem (time windows, two-dimensional loading constraints, satisfaction of the customer's requirement, etc.). In addition, while building the solution, if the generated route cannot be added to the under construction solution since one of the constraints of the problem will be violated, the route will be rejected. On the other side, if no route can be added to the solution under construction or there is a risk that the random method will turn around the unfeasible routes, we risk that the algorithm will run infinitely without finding a feasible solution. Therefore, we introduce a stop criterion defined by the elapsed CPU time, so the GA stops when this variable reaches the defined limit. To add a randomly generated route to the solution under construction, we have to know the filling plants, the truck that will make the trip and its maximum capacity and the customers who will be visited. The rest of the tour characteristics, in particular the start time, end time and number of racks assigned to each customer, will be calculated based on the solution method. In addition, the calculation of the start and end time of each route must take into account that the truck cannot move for another route until after the end time of its last one.
To generate an initial population of our distribution problem we decided to make a mix between two techniques, by generating in a completely random way a half of the population and the rest with a priori good solutions. To build a solution using the random method, we assign randomly selected customers to each truck until all trucks reach their maximum loading capacity. It should be noted that this assignment is made while respecting the various constraints of the problem (time windows, two-dimensional loading constraints, satisfaction of the customer's requirement, etc.). In addition, while building the solution, if the generated route cannot be added to the under construction solution since one of the constraints of the problem will be violated, the route will be rejected. On the other side, if no route can be added to the solution under construction or there is a risk that the random method will turn around the unfeasible routes, we risk that the algorithm will run infinitely without finding a feasible solution. Therefore, we introduce a stop criterion defined by the elapsed CPU time, so the GA stops when this variable reaches the defined limit. Add customer c to route 10 Add route of k to P i 11 p ← p + 1 To add a randomly generated route to the solution under construction, we have to know the filling plants, the truck that will make the trip and its maximum capacity and the customers who will be visited. The rest of the tour characteristics, in particular the start time, end time and number of racks assigned to each customer, will be calculated based on the solution method. In addition, the calculation of the start and end time of each route must take into account that the truck cannot move for another route until after the end time of its last one. The Random routes generator operates in a way that at each time it produces a random potential route that can be added to the solution under construction. On the other hand, if this route is not feasible the generator reproduces another one randomly which is different from the first until a feasible potential route is found. Otherwise, this mechanism cannot continue, so it has failed to find a feasible solution randomly. The greedy random method is similar to the random method except we do not generate the routes randomly. In fact, we associate a fitness with each route in order to be able to compare them.

C. Mutation and Crossover Operators
In the case of our adopted genetic algorithm we will try to reproduce the same mechanism of changing in the sequence as in biology. As a result, the proposed mutation operators should allow us to guarantee the homogeneity of the population and thus avoid too rapid convergence towards a local optimum. In addition, the developed mutation operators in this paper do not apply to all individuals in the population and in each generation. Indeed, we will define a mutation rate Pm which indicates what average proportion of the population must undergo a mutation.
Moreover, the use of completely free mutations allows a better exploration of the solution space. However, as in the case of population initialization, if we have a priori information about our problem, we can restrict the possibilities in order to achieve faster convergence. In this case, there is always the risk, if the information is not relevant, of guiding the GA towards a local optimum. Furthermore, the mutation operator behaves in a completely random manner in choosing the individuals subject to the mutation. As a result, we have the guarantee of the quality of the solutions obtained. In this paper, we propose a mutation operator, based on the local improvements illustrated in Figure 5.
The proposed crossover operator in this paper ensures the mixing and recombination of parental genes to form descendants with new potential. It corresponds to form two new chromosomes (children) based an exchange of genes between two reproducers (parents). However, we apply a single point crossover operator proposed by [34], which consists in randomly choosing a crossing point. The first child is obtained by taking the first sequence from the first parent, while the rest of the genes are taken from the second parent, so that the genes which have already been taken from the first parent cannot be considered again. It should be noted that the genes from the second parent remain in the same order ( Figure 6). Since the crossover operator is random, it occurs with a probability Pc, where we set the rate at 75%.

D. Fitness Function
Because of the variation of trucks capacity in the fleet the crossover and mutation operators can build solutions that do not respect the constraint of satisfying customer needs. On the other hand, our problem is an OVRP where the trucks are not attached to any filling plant at the end of their route. This possibility allows them to end their route at any other plant. Therefore, we encounter a problem of balancing the number of trucks between filling plants at the end of the day which generates a cost that we must consider in our fitness function. Thus, to obtain solutions which respect all the constraints resulting from our problem and after having applied the crossover and mutation operators we are facing a development of very complex mechanisms. Therefore, the best solution is to evaluate each individual by construct a penalization function.
In order to balance the number of trucks between filling centers at the end of the day, we make the least expensive routes, respecting the following criteria: A truck can move for the second time after all truck are used in the filling plant; a truck cannot move without a load.

Algorithm 3: Genetic Algorithm
Input: P i : Initial population, P c : Crossover rate, P m : Mutation rate Output: Best solution with minimum travelling cost 1 2 P i ← Generate 50% using the random algorithm 3 P i ← Generate 50% using the greedy random algorithm 4 5 Evaluate the population using fitness function 6 7 while s ≤ M axSteps do 8 9 Random select two parents from P i 10 Crossover operator(Parents, P c ) 11 Mutation operator(Parents, P m ) 12 13 Evaluate solution using the fitness function 14 15 s ← s + 1 We adopt the following indicator parameter: "The lower the fitness value, the better the individual". The fitness function is equal to the objective function plus the trucks balancing costs.
A selection operator used in our paper is based on Tournament method in which we select two individuals from the population in a randomly way, we compare their fitness function and we select the one how has the lower value.

E. Computational Experiments of the AG
In this section, we present the statistics and numerical results from the developed genetic algorithm. A heuristic was implemented using the Java language which is an object oriented programming language. In addition, the tests are performed on a machine with an Intel (R) Core (TM) i7 processor at 2.60 GHz (8 processors) and 16 GB of RAM.
To test the performance of our genetic algorithm as well as the power of the mutation and crossover operators used, we present an illustration of the statistics obtained during the tests carried out. For this, we will use a larger instance than the one presented in the exact resolution. The test instance considered in this part is derived from actual data provided by the same oil company. It is made up of 120 customers C = C1...C110 supplied by 16 filling centers D = D1...D16 and using at most 50 trucks of three different configurations.
We associate with each filling plant a number of heterogeneous trucks in terms of capacity. It should be noted that the tests are carried out for customer needs for a single day, but these needs influence either the number of trucks used  or the number of customers visited by the trucks or both. The solution obtained with the genetic algorithm depends on a number of parameters which are: population size, mutation rate and crossover rate. To study the influence of the population size we set the stop criterion at 100 iterations, crossover rate at Pc = 0.85 and mutation rate at Pm = 0.50. Figure 7 shows that the quality of the fitness function improves with the increase in the population size. In addition, we find a positive linear relationship between population size and the computation time. The graph illustrated in Figure 8 represents the distributions of the solution values over 37 iterations: the algorithm goes from a population of very dispersed solution at the beginning to a population more centered on the optimum found at the end, thus the improvement of the population is very fast at the start (global search) and becomes slower and slower as time passes (local search). Throughout the evolution of the genetic algorithm the crossover and mutation operators can build solutions that do not respect the problem constraints. However, in most cases these solutions are rejected by the algorithm mechanisms. With the aim of observing this behavior we define a reproduction error rate according to the number of iterations.

Error rate = N umber of invalidate individuals T otal number of individuals
The point cloud in Figure 9 shows an absence of connection between reproduction error and number of iterations. Moreover, the percentage of reproduction error varies between 2% and 8%, so the choice of the crossover and mutation operators is good depending on the probability of reproducing valid individuals. On the other hand, a good choice of the parameters of a method is very important, the adjustment of the parameters is crucial and conditions the quality of the results obtained. However, optimal tuning remains a difficult task. Fig. 9. Error Rate According to the Number of Iterations.
As described above, appropriate adjustment of parameters in genetic algorithms can make a significant difference in terms of performance. Some values can provide very high performance in specific instances while giving premature convergence in others, even over the same kind of problems. The parameters used in the comparison tests between the proposed GA and the exact method are listed as flow: • The population size equal to 1400 • The number of iteration equal to 100 • The crossover rate Pc = 0.85 • The mutation rate Pm = 0.50 Table III shows computational results of the developed genetic algorithm, which also makes comparisons with the exact method illustrated in table I. This table reports the configuration of each instance, the objective and CPU obtained from CPLEX running, solution and CPU obtained with GA and percentage Gaps. Based on results reported in table III we remark that both methods could find optimal solution for the four first test instances. The CPU time performs better in the case of GA for all test instances. In addition, the results obtained from the GA are not too far from those obtained by Cplex. This is justified by the choice of the crossover and mutation operators. Unlike the exact resolution of the studied open split delivery and pick up problem, the developed genetic algorithm gives ambitious results. We have performed tests for real instances of a large scale and in each time we obtain good solutions in a good CPU time. An example of a solution with five customers and two gas filling plants is depicted in Figure  10. It represents the test instance CL2-3 in which six trucks spread over the two filling plants. For example, we have the first truck which leaves from station 1, serves customers 3 and 4 and returns to the station 1 while the second truck which leaves from the same center, serves customer 5 and finishes its trip in the station 2. This is explained by the fact that the trip is open. On the other hand, the first truck which leaves from station 2, serves customers 1 and 3 and returns to the station 2 while the second truck which leaves from the same center, serves customer 1 then customer 2 and finishes its trip in the station 1.
Regarding the loading of each truck, it is clear that the combination is in the form of stacks of n racks of type L, M and N. For example, the truck which leaves from station 1 to serve customer 4 is loaded with five racks of type L, eight racks of type M and eight racks of type N.

VII. CONCLUSIONS
In this paper, we have addressed the open multi-depot split delivery and pickup routing problem with two-dimensional loading constraints in the case of a real-world problem of an LPG distribution company with its real-life assumptions and data. First, we have formulated our problem as a mixed integer linear programming (MILP), then we have tested the problem for 12 tests instances grouped in 4 classes in which we have proved optimality for 10 instances. Also, we have discussed the effect of changing the setting of the input data on the quality of the solutions obtained and on the running time. Due to the difficulty in solving the combination of two NPhard problems, such as vehicle routing and two-dimensional loading; efficiently, feasible solutions for large instances can be obtained only by developing a metaheuristic method like genetic algorithm. We developed a genetic algorithm and execute it for large data by increasing the size of instances which we are not capable of solving exactly within reasonable computational times. Finally, we gave some statistics and results for the developed GA and discuss obtained results.