Profitable Vehicle Routing Problem with Multiple Trips: Modeling and Variable Neighborhood Descent Algorithm

In this paper, we tackle a new variant of the Vehicle Routing Problem (VRP) which combines two known variants namely the Profitable VRP and the VRP with Mult iple Trips. The resulting problem may be called the Profitable Vehicle Routing Problem with Multiple Trips. The main purpose is to cover and solve a more complex realistic situation of the distribution transportation. The profitability concept arises when only a subset of customers can be served due to the lack of means or for insufficiency of the offer. In this case, each customer is associated to an economical profit which will be integrated to the objective function. The latter contains at hand the total collected profit minus the transportation costs. Each vehicle is allowed to perform several routes under a strict workday duration limit. This problem has a very practical interest especially for daily distribution schedules with limited vehicle fleets and short course transportation networks. We point out a new discursive approach for p rofits quantificat ion which is more significant than those existing in the literature. We propose four equivalent mathematical formulations for the problem which are tested and compared using CPLEX solver on small-size instances. Optimal solutions are identified. For large-size instance, two constructive heuristics are proposed and enhanced using Hill Climbing and Variable Neighborhood Descent algorithm based on a specific three-arrays-based coding structure. Finally, extensive computational experiments are performed including randomly generated instances and an extended and adapted benchmark from literature showing very satisfactory results.


Introduction
Nowadays, effective management of the supply chain is recognized as a determinant key o f co mpet it iveness and success for manufactu ring and d istribution organ izations. The problems of routing goods from depots to consumers or b et ween d ifferent log ist ical s it es are v ery imp o rtant. Drawing adequate plans may produce significant savings for many distribution systems. Therefore, the Vehicle Routing Pro b lem ( VRP ) is a well-kn o wn p ro b lem st ud ied in Operat ional Research. It deals with find ing the opt imal delivery routes configuration fro m one or several depots to a number of customers using a fleet of capacitated vehicles, while satisfying some constraints. The solution of a VRP is a set of min imu m cost routes wh ich fu lfill the custo mers' requirements. In the classical version of the problem, only one depot is considered with an unlimited fleet of vehicles where each vehicle performs exact ly one circu it. Several operational constraints can be considered in more practical applications of the VRP.
In this research, we focus on two variants of the VRP: the Profitable VRP and the VRP with Mult iple Trips, wh ich we combine to obtain a new variant not addressed in the literature to our knowledge and involving a new challenge in term of solving approaches. This problem which can be named the Profitable Veh icle Routing Problem with Multiple Trips (PVRPMT) calls for the determination of a set of routes for a given heterogeneous set of vehicles visiting a selective subset of customers such that: i) it is difficult, if not impossible, to visit all customers for a lack of logistical means or for insufficiency of the offer; ii) each visited customer p rovides a fixed profit that is recognized in the objective function to maximize. The latter is the difference between the total profit generated by the all v isits and the resulting total transportation cost; iii) each vehicle may perform several routes in the same planning period but does not exceed a strict durat ion limit (temporal capacity); iv) for each route, the total load does not exceed the vehicle capacity (physical capacity). The motivation to study this problem co me not only fro m its theoretical interest but also fro m its practical relevance. Fo r examp le, the size of the vehicle fleet indicates the capital invested in logistics resources. That the reduction of this size is naturally mo re desirable, the problem counts as a strategic choice in the problem reading to guide the optimization towards reducing both capital and operating (scheduling) costs. In this situation, there might be a trade-off between higher scheduling costs and lower vehicle ownership-based costs (and eventually associated driver salaries-based costs). Otherwise, this problem is very practical where daily drivers' schedules must be achieved with a limited vehicle fleet and short-distance distribution networks. Furthermore, when we are unable to respond to all costumer requests, the problem reading should be different. A subtle selection of costumers to be served is required. It should be based on a clear and measurable criterion. It is therefore important to define a significantly manner to classify customers according to their importance, behavior and costs implications, and to translate this into a global criterion that can be optimized as the total cost associated with the transportation program. This issue of profitability will be raised when the instances are to be conceived and tested, because most companies do not have a carefully designed evaluation of each costumer's profit. The latter should not only be related to the volu me o f its current request but must take into account other considerations, related to the future potentials of this costumer. And the value must also be reasonably projected on a coherent scale with the transportation costs that appear in the objective function in order to ensure significant aggregation-based assessments for the computed solutions.
The addressed problem concerns the PVRPMT as a new variant of pract ical VRP. Note that this kind of problems would be attractive in particu lar for p lanning mechanisms of ERP software wh ich must integrate a miscellaneous range of situations of management and control. For a mo re devoted emphasis ma inly to the operational context, our objective is mu ltip le with this problem co mb ining at hand: the maximization of the total collected profit, the min imization of the total routing costs and indirectly the maximization of the use of resources (vehicles and drivers). The issue of different mathematical models for this problem is to be addressed and used to seek solution optimality for s mall-size instances. Giving the co mplexity of the problem wh ich is NP-hard, two constructive polynomial heuristics with two different greedy-based constructions will be proposed as a first approach to solve the large-size instances. Their effectiveness is tested and then enhanced, in a second solving approach, with an iterative amelio rative procedure using successive local searches with mu ltiple and sequential neighborhood structures. For instances generation, we point out a new and quick discursive approach for profits quantification which is more significant than those existing in the literature.
The remainder o f this paper is organized as fo llo ws. Section 2 reviews the related research literature. In section 3, we formally present the optimizat ion problem PVRPMT with a detailed theoretical graph description and we introduce the corresponding proprieties and notations. Furthermore, we determine this problem co mp lexity. In Section 4, based on two strategies of sub-tours elimination constraints, we provide four mathemat ical models of the problem including MILP and 0-1 ILP. Section 5 describes two constructive heuristics and an enhancement procedure based on the Variable neighborhood Descent algorithm (VND). The co mputational experiment is presented in Section 6 to obtain insights about the different performances of the optimizing proposed algorithms. This section begins with a new precision about the evaluation of profitability indices to be incorporated into the nu merical instances' files and subsequently into the computation of the objective functions in adequacy with the measures adopted for the transportation costs. We draw conclusions and discuss future research directions in Section 7.

Literature Review
The problem we focus on is a co mbination of the two known variants such that: the vehicle routing problem with mu ltip le trips or with mu ltip le uses of vehicles (noted by VRPMT or VRPM) and the profitable vehicle routing problem (PVRP). To the best of our knowledge, these problems are studied separately in literature, excepting some special recent works [2,3,5], but with time windows considerations. For the latters, the concept of Multiple Trips is obtained indirectly by imposing a duration limit on each tour and not on the workday duration. Thus, the problem may look like a Capacitated VRP, where the multip le uses of vehicle is a consequence since the tours duration limit is fixed sufficiently small. So, the tours could be seen as parallel tours conducted by several vehicles. However, the integration of time windows constraint into the same problem gives the importance to the tours' order and consequently gives a special (restricted) Multiple Trips consideration. In this section, we devote two sub-sections to review the literature of the VRPMT and the PVRP.

The Vehicle Routing Problem wi th Multi ple Tri ps
The VRPMT is an extension of the classical VRP in which the same vehicle may perform several routes in the same planning period (workday). In this case, the fleet includes a fixed number of available vehicles. However, the duration of a vehicle workday which is made o f a set of successive routes does not exceed a certain limit .
This problem has a very practical importance. For example, in the home delivery of perishable goods like food, routes are of short duration and must be combined to form a complete workday [1]. So me real-world cases dealing with this kind o f problem appear in applicat ive research publications. For instances, note the paper of Brandao et al. 1997 [6] wh ich studies the VRP in Burton's Biscuit Ltd; and the case study of the logistical act ivities of Santa Fe Indonesia (precisely in the office o f Jakarta, a co mpany specialized on relocation services to individuals as well as companies) wh ich is considered in [7]. Derigs et al. 2011 [8] study some real cases dealing with this problem fro m a consultancy company in the air cargo transportation field M odeling and Variable Neighborhood Descent Algorithm (air cargo road feeder services). Planned and ad-hoc tasks, different resources and environment constraints are taken into account, such as driving time with breaks, working time with min imu m requirements in relat ion to the health and the safety of persons performing mobile road transport activities. These problems deals with different applicat ions, such as: the transportation of live anima ls to a slaughterhouse, the transportation of the meals or the flight food fro m a central kitchen to the aircraft, the transportation of goods from suppliers to retailers using one or mo re cross-docks. More additional constraints, such as time windows, split delivery, heterogeneous fleet, are considered respectively in [9, 10, 11and 12]. Grünert et al. 1999[13] develop a decision support system (DSS) which is designed in order to assist planners of the German postal service at the Deutsche Post AG. The ma in goal is to reduce the costs of the letter-mail transportation. In [14], a vehicle routing problem including mu ltip le trips in a health care organization which operates about 240 buildings of medical offices in Southern California is treated. Finally, we mention Hernandez et al. 2009[15] wh ich present a VRPMT applied in a context of optimizing the agronomic performance and its impact on the environment.
In the academic literature, the papers treating the VRPMT consider in general addit ional constraints main ly the time windows. Azi et al. 2007 [1] study this problem with time windows consideration and a single vehicle's use. The extension for mu lti-vehicle version is tackled by A zi et al. [2,3,4] and Macedo et al. 2011 [5]. Other works include also both multip le trips and time windows such as [6, 7, 8, 11, 12, 14, 15, 16, 17, 18, 19 and 20]. The overtime constraint is incorporated to the VRPMT in [14, 16, 21, 22, 23 and 24]. Heterogeneous fleet is considered with this problem in [7, 10, 12, 16 and 27]. So me VRPMT are studied by including a products compatibility constraint which forbidden to gather in a same vehicle route two or mo re different categories of goods [7, 8, 16, 24, and 27]. The mu ltip le trips constraint is used in VRP with backhaul in [7 and 12], co mbined with allo wed split delivery case in [12], with location problem in [25 and 26], with a planning horizon of several days in [14 and 27], with mult iple depots in [28], with meal breaks during the driv ing time in [11]. The problem is studied in a dynamic environment in [4 and 8].
In term of solving methods, few papers develop exact approaches for the VRPM T, because of the problem's complexity, such as [1, 2, 5, 15 and 17] using MILPs, column generation with branch-and-price algorith m and network flow-based models. Note that solving the VRPMT may consist on solving a routing and packing problem. In this level, different heuristics and metaheuristics are developed. For the routing phase, we can d istinguish some constructive heuristics used to obtain an initial solution, such insertion heuristics used in [3, 4, 10, 14, 20, 22, 27 and 29], the Clark and Wight algorith m [24, 25 and 31], the mod ified Clark and Wight algorithm [25 and 31], the nearest neighborhood method [25 and 29], clustering algorithms [7, 19 and 24], the sweep-based algorithm [20, 21, 28 and 32], set covering-based approaches [7 and 32] and the petal method [32]. The init ial obtained solutions are enhanced by different methods: we find the tabu search in [6, 13, 14, 18, 19, 20, 21, 23, 25, 27, 28 and 29], the simulated annealing in [3, 25, 26 and 31], neighborhood search based on insertion moves in [10, 18, 22, 24, 25, 29 and 31], on swap moves in [22, 24, 25, 26, 29 and 31], on 2-opt in [18 and 25], on adaptive me mory in [21 and 30]. Other specific neighborhood search algorith ms are proposed such as large neighborhood search in [3 and 4], guided neighborhood search in [8]. In [8], we find also a decomposition approach. In [18], a filter and fan procedure is developed. For the packing phase, some we ll-known heuristics for bin packing problems are used, specially the best fit decreasing in [16, 22, 24, 25 and 30]. In [8], a g reedy packing procedure is used and a fuzzy theory-based method is proposed in [33].

The Vehicle Routing Problem wi th Profit
VRP with profits are a generalization of the vehicle routing problems. Given a fixed-size fleet of vehicles, it might not be possible to serve all customers. Thus, a known profit is associated with each demand node and the customers must be chosen based on their associated profit minus the travelling cost to reach them in the solution. The idea to associate a known profit at each customer is made by Dell'A mico et al. 1995 [35]. The constraint to visit all the customers is relaxed, but for each unvisited customer a given penalty has to be paid. The objective function is to find a balance between these penalties and the cost of the tour. Profits are exp licitly considered both in the Vehicle Routing Problem (PVRP) and the Traveling Sales man Problem with Profits (TSPP), as stated by Feillet et al. 2005[34] in an excellent co mprehensive survey. The TSPP can be formu lated as a d iscrete bi-criteria optimization problem where the two goals are referred: maximizing the profit and min imizing the traveling cost. It is also possible to use one of the goals as the objective function and the other as a constraint. These problems can be div ided into three categories according to Feillet et al. 2005[34]: a) the Profit Tour Prob lem (PTP), b) the Selective TSP (STSP), and c) the Prize-Collecting TSP (PCTSP).
For the first problem (PTP), the objective is to maximize the difference between the total collected profit and the traveling cost. Feillet et al. 2005[34] survey lists various modeling approaches to TSPP and exact methods as well as heuristic solution methods. Archetti et al. 2009[36] study the capacitated version of the PTP and propose exact and heuristic procedures for it. More recently, the authors [37] develop an exact approach based on a branch-and-price algorith m. A restricted master heuristic is applied at each node of the branch-and-bound tree in order to obtain primal bound values.
For the second problem wh ich is known through three equivalent names in the literature: the Orienteering Problem (OP), the Select ive TSP (STSP), or the Ma ximu m Collection Problem (M CP), the objective is to maximize the collected profit such that the total traveling cost or distance does not exceed a certain upper bound. The mu lti-vehicles version of the OP is very similar to the VRP, except that vehicles are assumed to be with unlimited capacity and there is a time constraint on each tour. This problem is called the Team Orienteering Problem (TOP).
Recently, Vansteenwegen et al. 2011 [43] provides a survey of the existing OP and TOP. In this paper, the applications and the solving approaches about the two problems is reviewed listing various modeling approaches and exact as well as heuristic solution methods. In addition, many relevant variants of these problems are formally presented.
For the third problem which is named the Prize-Collecting TSP or also known as the Quota TSP concerns the determination of a tour with the minimu m total traveling cost where the collected profit or prize is greater than a lower bound.
Note that the PVRP is also a kind of PTP when considering comparing to TSPP that vehicles have physical loading capacities. In [38], The PVRP is applied in the reverse logistics where a firm aims to collect cores fro m its dealers. The problem is an extension of the classical mu lti-depots vehicle routing problem (M DVRP) in which each visit to a dealer is associated with a gross profit and an acquisition price to be paid to take the cores back. First, two mixed-integer linear p rogramming (MILP) are presented. Then, a Tabu Search-based heuristic is proposed to solve med iu m and large-sized instances. In [39], the research deals with a VRP in which the total profit is to be maximized subject to market competit ion. The PTP in a dynamic environment is considered in [40]. In this problem, the rewards (profits) are unknown for the customers which are not yet served. Indeed, the rewards depend on competitors' prices and auctions. In another extension of the problem, the profit is associated to each edge non to vertices. The problem is then called the Profitable Arc Routing Problem (PARP). Archetti et al. 2010 [41] study the capacitated undirected version (UCA RPP) and develop a branch-andprice algorith m, several heuristics based on Variab le Neighborhood Search (VNS) and t wo Tabu Search heuristics. Zachariadis et al. 2011 [42] propose another local search approach for the UCA RPP. Two solution neighborhoods are considered and the overall search is coordinated by the use of the promises concept.
The problem studied in this paper is a new extension which co mbines the two previous variants.

Proble m Descripti on and Notati on
We consider a comp lete undirected graph G = (V, E), where V= {0… n} is a set of vert ices and E is a set of edges.
Vertex 0 represents the depot and a fleet of vehicles k = {1… m} is based. Each vehicle has a limited capacity Q (or Q k if with heterogeneous fleet) and a ma ximu m nu mber of trips L. Note that, the parameter L is defined to facilitate the problem formu lat ion and can be used as a real constraint, i.e. as a limit to be respected. However, the best value of L can be experimentally identified. An edge (i,j)∈E represents the possibility to travel fro m customer i to customer j. A non-negative demand q i , profit p i , and time service S i , are associated with each customer i (with setting p 0 =q 0 =0). A travel time t ij and cost c ij are associated with each edge (i,j)∈E. Each vehicle starts and ends its tour at vertex 0, and can visit any subset of customers with a total demand that does not exceed the used vehicle capacity Q. In addition, there exists a time horizon denoted by the duration limit T max which establishes the duration of a workday. It is assumed that all parameters are nonnegative integers and the environment is determin istic.
This problem named the PVRPMT consists on determining a set of routes and to assign each route to one vehicle, such that the same vehicle can be used by several routes while respecting the time horizon capacity. The objective is to maximize the difference between the total collected profit and the cost of the total traveled distance. Note that the following properties: • The optimal solution may be co mposed by a subset of customers.
• Each route starts and ends at the depot, • The total customers' demand in the same route does not exceed the physical capacity of used vehicle, • The duration of routes assigned to the same vehicle does not exceed Tmax.
• The profit associated at each customer is fixed and can be collected by any vehicle.

Proble m Complexi ty
Olivera and Viera [21] proved that the VRPMT is NP-hard as well as the PVRP [34]. The studied problem represents the combination of these two NP-hard problems. This makes it also NP-hard. In addition, PVRPMT is a generalization of the classical VRP. There are instances of PVRPMT in which there exist enough vehicles in the fleet able to optimally visit during the workday all customers using one vehicle for each tour. These cases are naturally reduced to a classical VRP. As the VRP is an NP-hard [22], the PVRPMT is also NP-hard.

Mathematical Models for the PVRPMT
The design of the VRP solution stands against the presence of the sub-tours. For that, different sub-tours elimination constraints are proposed. The most classical constraint and the most used in the literature can be written in this way: M odeling and Variable Neighborhood Descent Algorithm This formu lation is exponential and influences strongly the resolution time.
In this section, we propose four mathemat ical formulat ions for the PVRPMT. The difference between them is based on the strategy used to eliminate the sub-tours and the decision variables defin itions. The different indices, parameters and decision variables are g iven in Table 1.

Modeling wi th 0-1 Integer Linear Programmi ng
We start with the idea to specify for each assigned customer his order in the trip. This idea has been applied in the scheduling problem. For each job, we determine the position of the job in the sequence. We use the visit order of customer i in the trip in order to eliminate the invalid tours. Accordingly, we eliminate the subset S and the associated constraints. We consider the following binary decision variables: A new decision variable is added and informs about the visit order of customer i in the trip l.
Based on the choice of the principal variable in the formulation, we can distinguish two different mathematical models.
Here, the second variable is used just to make the connection between the edge and the vertices. The resulting formulat ion is the following: In this formulat ion, the objective is to maximize the overall collected profit minus the transportation cost. The constraints (3)(4)(5) guarantee that each customer is visited at most once. In (6), if the route l exists, it should start and fin ish in the depot. (7) represents the capacity constraint. The limit duration on a workday is restricted by (8). (9) represents the flow conservation constraint. The constraint (10) establish the relation between the edg (i,j) and the relevant position of i and j. (11) represents the initialization of the counter and (12) stands against the addition of a customer if the route is closed. In (13), each constructed trip should start in the depot. (14) represent the integrity constraints.
• Second formulation : 0-1 ILP2 with as principal In this model, we conserve the constraints which establish the relation between and such in ILP1. When it is possible, we integrate the variable in the objective function and in the constraints. The objective is to test the influence of this transformation on the upper bound and to know which exp ression can conserve much mo re informat ion if the integrity constraints are relaxed. The resulting model is as follow:

Modeling wi th Mixed Integer Li near Programming
To overcome the limitation of the classical sub-tours elimination constraint, Miller et al. 1960 [44] p ropose a new ones which are corrected by Kara et al. 2004[45]. In this level, we propose an adaptation of these constraints to our problem. For that, we remove the decision variable and we define U i . as an associated variable to customer i used to reformu late the sub-tour elimination constraints; Firstly, we present the ordinary model. Then an extension of this model with some proposed cuts.
• Third formu lation: M ILP1 (without cuts) The formulat ion is the following: The constraint (29) guaranties that the customer i is visited at most once during the workday. (30) represents the capacity constraint. The workday duration limit is respected in (31). (32) represents the flow conservation constraint. The adaptation of Miller et al. sub-tours elimination constraints, as it is mod ified by Kara, for our p roblem is given by (33) and (34). (35) represents the integrity constraint.
• Forth formu lation : MILP2 ( with cuts) For the previous mathemat ical model (M ILP1 without cuts), we add an optional variab le δ kl wh ich informs about the used vehicle. So, the correspondent constraints which establish the relation between the two decision variables will be adjo ined.
The formulat ion becomes: To integrate the new decision variable, three additional constraints which establish the relation between the t wo decision variables are adjoined. (43) means that if the route l exists, it should start and fin ish in the depot. (44) guaranties that if the edge (i,j) is assigned to the trip l, this latter should be constructed. The opposite case is presented by (45). This constraint prohibits the construction of empty route, i.e. if the route l is constructed, at least one edge must be assigned at this trip.

The Heterogeneous Fleet Case
In this subsection, we proposed a formulation extension to solve the special case of PVRPMT with heterogeneous fleet. So, each vehicle k has its own capacity . Then, some modifications on MILP2 are done. The constraint (4.

Solving Algorithms
The present section is devoted to the description of the ma in steps of the imp lemented algorith ms. We start by explaining the solution coding scheme. Then, constructive heuristics used to build init ial feasib le solutions to the problem are described. Finally, imp rovement procedures based on Hill Climb ing and Variable Neighborhood Descent (VND) are proposed.

Solution Coding Scheme
To further speed up the computation, we use tree-array data structure to represent a solution. The information concerning the customer is stored into the first array (V1[j]). So, its size is equal to the number of customers n (we eliminate the deposit fro m the solution representation). First, the visited customers are stored adjacently in the order of visit. Then, the unvisited customers are inserted in the end of the array. It is important to know the start and the end of each trip. This information is done by the second array (V2[j]). This latter is a b inary vector. Finally, the third array (V3[j]) informs about the index of the used vehicle to visit the correspondent customer. Using this structure, the lecture of the solution is clear and easy, and changes can be performed very quickly and in a constant time. The solution is coded as follows: Let the matrix V[i][j] with i ϵ {1, 2, 3}and j ϵ {1, 2,…n } be the following:

Constructi ve Heuristics for the PVRPMT
The goal of this subsection is to construct good feasible solutions for the problem. We propose two greedy constructive heuristics. These heuristics use some local optimalit ies in certains steps of the algorithm.

Heuristic H1
The insertion heuristic is used to build an init ial feasible solution. In every iteration, the procedure evaluates all feasible insertions of unvisited nodes and selects the node representing the best insertion. An insertion is evaluated with the following criterion. Let i be some node in a tour and let j be a node candidate for the insertion. Let Chg act and Tm act denote the actual charge and the traveled time for the current used vehicle respectively. The insertion of j is feasible if the capacity constraint and the constraint of time duration are verified (i.e. (51) and (52) are verified): (51) (52) Then, the best insertion is determined by the pair (i, j) for which p ij is maximu m: Where Pij the profit generates by the insertion of the customer j after the customer i, p j the profit associated with customer j and c ij the cost transportation between i and j. The route is considered as completed if the vehicle can not receive any mo re other costumers due to its physical capacity. The vehicle returns to the depot and a new trip starts with the same used vehicle in the last iteration. The vehicle is changed only if the daily time horizon T max will be violated by the addition of a new trip with the same vehicle. In this case, the next trip is done by the next vehicle (see Figure 2).
Remarq 1: in the case of a heterogeneous fleet, the vehicles are sorted in the decreasing order of their physical capacities. This order is used as a priority rule for the vehicle choices.

Heuristic H2
This heuristic H2 is almost identical with heuristic H1 with the following basic difference : to construct the trip, we use the same prev ious procedure, but the used vehicle is not necessary kept in prior for the next trip. At each iteration, we choose the vehicle which has the longest remaining time service by breaking ties with the largest capacity order (see Figure 3).  Remarq 2: in the case of heterogeneous fleet, we use the same priority rule for choices used in H1.

Improvement Sche me
For the PVRPMT, it is essential to have a neighborhood that changes the visit combinations for customers. Three kinds of moves including insertion move, swapping move and Cross exchange operator move are used to define a set of neighborhoods that allow the exp loration o f increasingly distant solutions from the incu mbent to overcome local optimality and strive for global optimality.
Our improvement phase consists on first developping three Hill Climb ing algorithm (HCi, HCs, HCc) using the three neighborhood structures, then we test the influence of combin ing the three last algorithms in a Variab le Neighborhood Descent procedure.

Hill Climbing Procedure
The local search algorith ms show its performance to solve various variants of routing problems. Here we use the Hill Climb ing heuristic wh ich belongs to the family of local search methods which often built on neighborhood moves that make small changes to the current solution.
The insert move consists of remov ing one customer fro m a current position j (origin position) and putting him into another new position k. The destination route (i.e. the route that contains the new position) can be an existing route or a new one. To take into account all the possible moves, a new decision variable α is defined.
All depends on the positions of j and k and the variable δ, we can distinguish different possible moves.
＊ the swap move consists on exchanging two customers. In the new solution just the position of the two selected customers will be exchanged ( The cross-exchange operator consists on interchanging non-consecutive customer segments between the same or two different routes with the restrict ion that the orientation of them be maintained. Remarq 3: Note that , to elliminate the neglected move, for the two first neighborhoods (insertion and swap), at least one fro m the two selected customers should represent a visited customer in the initail solution. For the cross-exchange operator, the two selected customers should be visited in the initail solution.

Variable Neighborhood Descent Algorith m
The second idea is to test the influence of combining the    = trip next the in inserted be will j suiv trip last the in inserted be will j prec : : It exists a vehicle with a " remaining time" sufficient to be used (vehicle availability) Step1: select the available vehicle (according to the priority list) and take the depot position.
Step2: select the next customer according to the local optimality rule The VND is the simp lest variant of the variable neighborhood search heuristic (VNS) which performs several descents with different neighborhoods until a local optimu m for all considered neighborhoods is reached. Let N 1 , N 2 , … , N k denote a set of K neighborhood structures (i.e., N i (S) contains the solution that can be obtained by performing a local change on s according to the i th type). The VND works as follo ws: Step 1. Choose an initial solution s in S.

Computational Results
In the following, obtained results concernening the performance of the different proposed mothods are reported. First, the test instances construction, the parameters choice and the profit quantification are presented. Then, the performance of the equivalent presented models are tested and compared. Finally, the last subsection is devoted to the heuristics procedures performance (contructive, Hill Climb ing, VND). The M ILPs and solving algorith ms was tested on a Intel Core 2 Duo CPU 2.20 GHZ and 4.00 Go RAM. The codes are written in C++, using CPLEX lib rairies for the first part. A ll the algorith ms were stopped before a computational time of one hour at atmost.

Test Instances and Parameters Choice
The tests will be applied first on our own benchmark devoted for small-size instances and then on the benchmarks of Taillard et al. 1996[23] taken fro m the VRP lib rary with adaptations and some extended data. In our benchmark, twenty small-size instances are generated randomly. A certain setting is used to obtain interesting instance values according to the practice and real-case situations. For each instance, we indicate the number of vertices which ranges fro m 6 to 20. The customers are randomly distributed in two-dimension area, and the depot is set at point (0, 0). the fleet of vehicles accounts for 2 or 3 vehicles. The vehicle physical capacity limit ranges from 1000 to 3000 kgs. The horizon t ime limit for all the s mall size instances is equal to 480 minutes (i.e. 8 hours of workday).
Note that in the litterature according to [36], the profit of customer depends on the three parameters namely: cons, h and the customer's demand , where h is a random rat io number uniformly generated in the interval[0,1] and cons is a constant factor that measures the profit according to a sale turnover. In [36], the authers consider p i = (cons + h)qi and suppose that cons=0.5 indifferently with the level of greatness of the other values used for the instances data. This implies a lack of guarantee of a necessary coherence between the different values while the choices of the units and the level of greatness of the four used data are unambiguous: costomers' demands q i , distances d ij , transportation times tij and transportation costs c ij .
In our opinion, it is important to obtain meaningful and significant proportions of the profits according to the transportation costs. In our case, they are assimilated to d ij . Concerning the profit calculation, we take as reference a general realistic model where the logistical cost represents between 5 and 10 per cent of the sale turnover, and the gross benefit may represent between 20 and 50 per cent of the sale turnover. So, it is possible to have a g lobal profit which ranges almost between 3 and 5 times the gobal transportation cost. The distribution of this global profit on customers may respect the demand quantity q i of each customer i. Let q max and q min be respectively the greatest and the smallest value of all demand quantities. q max and q min can be associated respectively to α max =5 and α min =3 factors. We use a pro jetion to calculate the factor α i of the customer i according to his demand q i .
In order to estimate the transportation cost, we calcu late firstly the average distance, identify and α i and p i such as: Note that t ij are also assimilated to d ij and the level of greatness of the vertices coordinates values is chosen adequately to obtain significant t ij values suitable with T max .
In addition, we point out that L, wh ich is used in MILP formulat ions, represents the total number of trips wh ich can be made by a vehicle during T max . Its value can be estimated as : For the original instances of Taillard et al. 1996, we add all previous adaptations and extensions, as well as the service time S i and the maximu m nu mber of t rip wh ich can be assigned at each vehicle L.

Mathematical Models Experimentati ons
In order to test the proposed models, we used the commercial solver CEPLEX 10.0. of ILOG ®.
First, the four mathematical models are teted on the small size instances and the results are reported in Table 3.The optimal solution obtained by CPLEX is indicated under column f . The "_" means that we cannot obtain the optimal solution before the time limit. The colu mn UB represents the solution of the linear relaxed problem (the upper bound). In addition, (%) represents the gap between the optimal ( ) We observe on Table 3 that the optimal value can be obtained within the time limit just for small size instances for the four mathematical models. Fo r the instances where the number of vertices is superior to 17, the optimal solution cannot be determined. We can show also that computationnel time grow strongly with the increase of instance size. M ILP 2 is ab le to solve more instances than the other proosed models and it clearly outperforms them in terms of co mputing time.
In Table 4, we try to co mpare the performance of the four mathe matica l mode ls by confronting the obtained upper bounds. We choose as reference the upper bound of the first mathematical model (UB1) and we calculate the deviation of the different upper bounds in co mparison with the first. (%*)w =(UBw-UB1)/UB1× 100 .
Different remarks can be taken into account. First by comparing the two strategies, the 0-1 ILP and the MILP, it is very clear that the upper bounds obtained by the mixed integer programming are better than those obtained by the 0-1 integer programming fo r all the tested instances. For the 0-1ILP, the choice of the principal variable has not great influences on the upper bounds quality (for the majority of the tested instances the upper bound obtained by the two models are equal) but strongly affects the nu mber of iterations. For the MILP, the upper bounds of the mathematical models with cuts are better than those without cuts. That shows the efficiency of these cuts. MILP2 is ab le to solve more instances than the other mathemat ical models and outperforms them in terms of computing time and the upper bound quality. So, we can judge that the additional constraints represent valid cuts for our model.

Heuristics Procedure Performance
To test and compare the performance of our heuristics, we compute the obtained gaps of the obtained solutions comparing to the linear upper bounds. These gaps are given: first for the two constructive heuristics bounds, second for the three Hill Climb ing based solutions, and finally for enhancement by using the Variable Neighborhood Descent algorith m. Gaps are also computed comparing to the optimal solution for the s mall-size instances. Then, to know the contribution generated by the VND algorith m, we calculate the gap between the initial solution given by the constructive heuristics and the lower bound. As it is previously ment ioned, MILP2 represents the best formu lation of our problem. So, the lower bound given by our heuristic will be co mpared with the upper bound obtained with this formulat ion.
The test includes our twenty instances generated randomly with adequate settings and forty adapted instances selected fro m the bench marks of Taillard et al. 1996[23]. The results are shown in table 5 and 6.
The column f and UB represent respectively the optimal solution and the upper obtained by the linear relaxation of MILP2. The init ial solution of our constructive algorithms (Constructive heuristic 1 and constructive heuristic 2) is given in column H1 and H2. The Hill Climb ing improvement tested for the three neighborhood structures by the insertion move, the swap move and the cross-exchange are showed in the column HCi, HCs and HCc respectively. The column VND gives the lower bound obtained with the Variable Neighborhoods Decent algorithm. The column gap and gap* represent respectively the gap between the lower bound (the VND solution) and the optimal solution and the upper upper respectively. The column gap** calculates the enhancement generated by the VND algorith m by mesuring the gap between the VND and the initial solution.
Fro m the first observation and by comparing the three Hill Climb ing shemes, we can see that these procedures give several acceptable results for the most tested instances. Howerever, some of them still require further enhancement. On the 12 instances for which the optimal solution is determined, the VN D algorith m is managed to find the optimal solution in 50% of cases. For the others, the gap is generally tiny, except for few exceptional instances, where it is quite signeficant reaching for the worst 33%. For the large size instances, this gap (gap*) becomes important and grows continuously because it is co mputed compared to the linear bound and not to optima. The VND algorithm enhances clearly the init ial solution for the total of the tested instances. But, for few instances (i.e. instances 7 and 8) the improvement is small and the VND converges quickly to local ma xima.
To crown all, we can conclude that the MILP2 (with cuts) represents the best formulation of our p roblem. The two constructive heuristics and the enhancement procedure based on the Hill Climb ing and the VND algorithm produce acceptable solutions close to the optimal ones for s mall size instances. For the big size instances, the obtained solutions are the best till now. In a future works, we think about enhancing more and more the upper bounds with some polyhedral techniques to have a clearer idea about the performance of our proposed methods for the large size instances. To escape from a current local optimu m, we should think to add, in a third solving phase, some perturbation moves to strengthen the search diversification in our algorith ms. Best performances could be deduced in future works after some enhancements.

Conclusions
In this paper, we describe a new variant of the vehicle routing problem namely the Profitable Vehicle Routing Problem with Multip le Trips. Two d ifferent strategies of sub-tours elimination constraints are used and for each strategy two different cases are defined. Thus, four mathematical models are obtained looking for optimal solutions. For large-size instances, two greedy constructive heuristics are proposed in a first solving step. Three Hill Climb ing algorithms based on three neighborhood structures are developed. With a VND procedure, these methods are managed to obtain imp roved solutions in a second solving step of the problem. The design of these methods is based on elements of reasoning to obtain intrinsically the best possible solutions fro m the first iterat ions. Enhanced diversification was covered by a broad research approach to effectively improve the results. Experimental study shows satisfactory results for small-size instances with MILPs using some cuts. Two strategies of sub-tours elimination constraints are used representing a good idea to formu late this kind of constraints. The empirical results show the performance of the proposed constructive heuristics which provide quick solutions very close to the optimu m, and also a satisfactory enhancement by using the improvement p rocedures. In future work, with some adjustments and the introduction of well-studied perturbation moves, the obtained results could be further refined to ensure a better optimality of solutions.