Multiobjective car relocation problem in one-way carsharing system

In this paper, we present a multiobjective approach for solving the one-way car relocation problem. We fix three objectives that include the number of remaining rejected demands, the number of jockeys used for the relocation operations, and the total time used by these jockeys. For this sake, we propose to apply two algorithms namely NSGA-II and an adapted memetic algorithm (MA) that we call MARPOCS which stands for memetic algorithm for the one-way carsharing system. The NSGA-II is used as a reference to compare the performance of MARPOCS. The comparison of the approximation sets obtained by both algorithms shows that the hybrid algorithm outperforms the classical NSGA-II and so solutions generated by the MARPOCS are much better than the solutions generated by NSGA-II. This observation is proved by the comparison of different quality indicators’ values that are used to compare the performance of each algorithm. Results show that the MARPOCS is promising to generate very good solutions for the multiobjective car relocation problem in one-way carsharing system. It shows a good performance in exploring the search space and in finding solution with very good fitness values.


Introduction
Carsharing is a new concept for using a car by sharing it between multiple users. This kind of systems turns the mobility into a service that allows its users to enjoy the advantages of owning a private car without the need to buy it. Carsharing operators maintain a fleet of cars that are available for use to their members who usually pay yearly registration fees, a monthly membership and are billed on the usage according to time spent and distance traveled.
What makes the carsharing more convenient is that members do not have to care about insurance, maintenance and even fuel price, which are included in the usage fees and shared between all members. Usually carsharing companies provide a short-term car rental. Members can have access to vehicles on the go or according to a reservation in advance via telephone, smartphones or website; cars are available 24 h a day, 7/7. The fleet of cars is usually spread in a way that every potential member (user) can have a car in his neighborhood at a convenient walkable distance.
There are multiple models of carsharing systems. For instance, the free-floating model allows users to pick up and return their rented car to any legal parking space within a defined area. This kind of system is convenient to users, since they do not have to stick to defined stations to park or to pick up a car. Unlike free-floating carsharing systems, the round-trip model is a station-based carsharing system where the user has to return the car to the station of departure. Finally, one-way carsharing provides an interesting compromise between two previous approaches. It is another model of station-based systems where a user can drive a car from a station and ends his reservation at any other station in the system.
In our paper, we focus on one-way carsharing system. The one-way flexibility makes the system more attractive to its users. However, this has a major drawback. It often leads to an imbalance in the number of cars in stations. On one side, some stations will reach their maximum parking capacity, forcing users who want to return their car in these stations to look around for other stations that have available parking places. On the other side, other stations will not have any available car for a new arriving user who wants to take a car from these stations. When the imbalance problem occurs often, the carsharing system becomes less attractive since it is unreliable for users. Eventually, it may lead to service failure.
In this context, carsharing operators are working hard to alleviate this problem. For this sake, carsharing companies recruit employees to redistribute cars on the stations based on the users demand for cars. In this paper, we refer to these employees by jockeys and the redistribution operations by relocation operations. Car relocation operations can improve the system reliability but in the same time, they increase the cost of operating the system. Here emerges the need for optimizing these operations in order to reduce the cost of relocating the cars and to increase the system reliability. A relocation operation consists of moving a car from a station that has unused cars or from a station that needs free parking spaces for near future demands, to another station that has a shortage in the number of available cars for users who want to start their trip starting from this station.
It is good to know that bikesharing systems have the same mobility concept and eventually have the same relocation problem. However, the size of a bike is much smaller than the size of a car. Consequently, the bikesharing operator can relocate many bikes using a truck, which is not possible in the case of cars. In this paper, we consider that a jockey relocates one car at a time.
Many research papers have treated the car relocation problem. In [1], a study is presented for the relocation problem in one-way electric vehicle sharing system in Praxitele, France. The authors suggested a fleet of finite capacity to redistribute cars over the stations. The problem was modeled as pickup and delivery and formulated as a mixed-integer programming. After that, Barth and Todd [2] developed a queuing-based discrete event simulation model to analyze the performance of multiple station carsharing system with one-way trips. In 2006, Kek et al. [3] developed a time-stepped simulation model to evaluate the performance of one-way carsharing system. In order to define the best combination of system parameters, Kek et al. [4] proposed a new approach that is based on the simulation model developed in [3]. The set of recommended parameters suggested by the three-phase decision support system could lead to a 50% reduction in the staff cost, an improvement in the zero-vehicle time (empty stations) and full-port time (full stations) and decrease in the number of relocations between 37.1% and 41.1%. In another study, Cucu et al. [5] proposed another approach to help decision makers to maintain a balanced number of cars in each station and to evaluate the potential of new locations for carsharing stations. The authors developed a fuzzy logic algorithm that is based on the preferences of potential customers (departure time, weather conditions, day of the week, and traffic conditions) to forecast clients demands. The authors concluded that a good estimation of client demands is crucial for improving the level of service and relocation operations of a carsharing system. Nair and Miller-Hooks [6] developed a stochastic, mixed-integer linear programming model with a joint chance constraint in order to plan the vehicle relocation operations. This model aims to minimize the cost of these operations and to guarantee the satisfaction of a proportion of all near-term demand scenarios. The study showed that the proposed relocation operations that consider the demand uncertainty could highly increase the level of service of the carsharing system. Smith et al. [7] presented a study for the problem of relocating the employees in charge of the relocation operations in a one-way carsharing system. The study that is based on a fluid model focused on the optimization of the relocation operations in order to reduce the number of relocation vehicles and the employees who drive these vehicles. The simulations suggested that the number of relocation employees required for the relocation operations should be between 1/3 and 1/4 of the number of the vehicles in Euclidean network topologies. It is good to mention that many papers studied the case where the client himself is used to relocate vehicles by offering them incentives such as additional free time or bonus credits [8] or cheaper dynamic prices for their trips [9][10][11][12]. Although client participation in the relocation operations can contribute in reducing the number of required relocation operations, this is often not guaranteed [13,14]. In more recent papers, Boyaci et al. [15] tackled the problem of the optimization of vehicle and personnel relocation operations in one-way electric carsharing systems with reservations. They used a multiobjective mixed-integer linear programming optimization and a discrete event simulation framework. Their results showed that efficient algorithms play a major role in relocation operations. Also, Huang et al. [16] presented a mixed-integer nonlinear programming model (MINLP) to locate stations and determine their capacity for one-way carsharing systems taking into account the car relocation operations with flexible travel demand. A gradient algorithm is used to solve the MINLP. When applying the model on a real-world system in China, results showed that the algorithm can be used to solve real-world large-scale problem. Moreover, Xu et al. [17] developed a mixed-integer nonlinear nonconvex model for the electric vehicle fleet size for one-way carsharing systems taking into account the vehicle relocation and personnel assignment. The model was tested on a one-way carsharing system in Singapore to prove its efficiency. In this paper, we first present the transition from mono to multiobjective formulation of the studied problem. Then, we propose a multiobjective formulation for solving this problem. After that, we present the resolution methods that we use to solve the car relocation problem. Next, results are shown and finally a conclusion and perspective are provided. To our knowledge, papers that deal with the multiobjective aspect of the car relocation problem in one-way carsharing system are really scarce.

From mono to multiobjective
A typical formulation of the car relocation problem in oneway carsharing system with its Integer linear programming (ILP) model and related constraints can be found in [18]. Below, presented the related objective function (1): Several decision variables are used in (1) : • move e i t j tþt ij : When an employee e is involved in a moving activity, the associated binary variable is assigned with the value 1, while it remains 0 otherwise. • rel e i t j tþt ij : When an employee e is involved in one of the relocation activities, the associated binary variable is assigned the value 1, while it remains 0 otherwise.
• u e : When an employee e is used at least once during the day, its associated binary variable takes the value 1, while it takes the value 0 otherwise. • outR i t : This variable can be assigned with integer values. It represents the number of rejected user demands to rent a car from station i at time t. • inR i t : This variable can be assigned with integer values. It represents the number of rejected user demands to give back the rented car to a station i at time t.
On other side, the input parameters used in the objective function (1) are listed below: • cMove : Denotes the estimated cost of a relocation or moving activity from a station to another one. • cRejectEmpty : Stands for the estimated cost of the rejection of a demand to rent a vehicle from a station. • cRejectFull : Stands for the estimated cost of the rejection of a demand to give back a vehicle to a station. • cStaff : Denotes the estimated cost of an employee during a day.
In such approach, the objective function (1) minimizes the weighted aggregation of the number of rejected demands to take or to return a car, the number of employees and the number of move and relocation operations needed to reduce the number of rejected demands. Related constraints are used to make sure that each employee is involved in only one activity at a time, while others are used to make sure that an employee will not start a new activity until he finishes the previous one and to guarantee a continuity for the employee's activities. In addition, some constraints are used to calculate the number of available vehicles at each station at each time step and to make sure that the number of available vehicles at a station cannot be greater than the capacity of the station. And finally, there are constraints used to ensure that the number of rejected demands at a station cannot exceed the demand itself. This problem is addressed in an exact way in [3] on small instances. For bigger instances, a greedy algorithm [19] has proved its performance regarding the execution time and the quality of the generated solutions in a comparison with an exact solver (CPLEX). As a result, these methods return a planning for the operations that a jockey should perform in order to reduce the number of rejected demands. A relocation planning consists of a schedule for the car relocation operations that should be done during the day. Two main drawbacks can be noticed regarding this formulation: • First, when we solve the car relocation problem using the two approaches described earlier, we consider that a jockey is working all the day, regardless of the number of operations he is doing. However, while analyzing the planning, we notice that the schedule for some jockeys is scattered all over the day. For example, even if a jockey has to do two operations only, one at the beginning of the day and another one at the end of the day, the jockey will be counted as any other jockey that has a schedule full of relocation operations during the whole day. For that reason, when we present the solutions provided by the two developed approaches for the decision makers, we are missing one important aspect of the problem, which is the total time spent by jockeys for the relocation operations plan. Adding a new constraint for the integer linear programming model makes the problem harder and getting results for real dataset using an exact solver in reasonable time is not feasible currently. • Second, the mono-objective approach considers in the same objective function, the cost of a jockey (cStaff ) and the cost of a rejected demand (cRejectFull and cRejectEmpty); this latter which can be seen as the carsharing price. Hence, it supposes to be able to fix these costs. Despite the importance of the cost parameters, in [18] those costs were left behind for simplicity reason. Setting different cost values will change the entire results of an ILP since we have one objective to optimize. But it should be noted that while cStaff is somehow easy to determine, the cost of a rejected demand is very difficult to evaluate. We can calculate the loss of not satisfying a demand; however, we also have to calculate the negative impact on the image of the service and eventually the risk of not using the service in the future. For these two reasons, using a multiobjective approach can be seen as a powerful alternative to the existing algorithm. In this approach, each objective is treated separately and independently from other objectives. Thus, it is not necessary to take costs into consideration when they are the same for each rejected demand, for each relocation operation and for each jockey. However, if costs were not the same (from one trip to another, from a relocation operation to another, from a jockey to another), it will be necessary to integrate them in the objective functions (for example, weighted objective functions, with non-constant costs, may lead to obtain solutions that prioritize solving rejected demands having high carsharing prices).

Integer linear programming formulation
The relocation problem can be modeled as a two-dimensional time-space matrix of size N Â M, where N is the total number of stations S ¼ f1; 2; . . .; Ng and M is the number of time steps T ¼ f1; 2; . . .; Mg in the day. Each element of the matrix represents a station s i at time t. For each station s 2 S we generate M nodes to represent that station at each time t. Then we put all the S Â T nodes in one set V ¼ f1 1 ; . . .; 1 T ; . . .; N 1 ; . . .; N M g. During the day, we consider that an employee is always involved in one of these three types of activities: 1. Relocating: is the action taken by the jockey to move a car from a station i to another station j. 2. Moving: is the action taken by the jockey to move himself from his current station to another station in order to begin a relocation activity. 3. Waiting: when the jockey is not involved in relocating or moving activities we say that the jockey is waiting.
Therefore, to represent these activities we generate three sets of arcs in the time-space network. For each node i t 2 V, we construct an arc wa that represents a waiting activity between i t and i tþ1 ; we call this set WA ¼ f. . .; waði t ; i tþ1 Þ; . . .g. Then, for each node i t in V, we construct N À 1 arcs ma to represent moving activities between station i and j, 8 i, j 2 S, i 6 ¼ j, from time step t to time step t þ t ij where t ij is the number of time steps needed to go from station i to station j; we name this set MA ¼ f. . .; maði t ; j tþt ij Þ; . . .g. In the same way of moving activities, we build N À 1 arcs ra to represent relocation activities for each station, and we denote this set RA ¼ f. . .; raði t ; j tþt ij Þ; . . .g. We represent the available staff that will be involved in doing these activities by a set E ¼ f1; . . .; e; . . .; Wg where W is the maximum number of available employees. The ILP model described in [18] uses six decision variables, the variables related to the present paper here are: • u e : Binary variable, that takes the value 1 if the employee e is used during the day and 0 otherwise. • rel e i t j tþt ij : Binary variable associated with the set of relocation activities RA. It takes the value 1 if employee e has been relocating a car from station i to station j, from time step t to t þ t ij and 0 otherwise. • outR i t : Integer variable to represent the number of rejected demand to take a car out of a station i at time step t. • inR i t : Integer variable to represent the number of rejected demand to return a car into a station i at time step t. The mathematical model for the car relocation problem in one-way carsharing system was solved with a greedy approach in [19]. The proposed greedy search method proved its performance regarding the execution time and the quality of generated solutions in a comparison with an exact solver (CPLEX). However, a simple analysis of the results shows that the proposed method did not treat the multiobjective aspect of the car relocation problem. In the greedy search algorithm, the objective was to minimize the number of rejected demands. However, it is important to integrate other objectives. In fact, the proposed methods return a planning for the operations that a jockey should perform in order to reduce the number of rejected demands. A relocation planning consists of a schedule for the car relocation operations that should be done during the day. On the other hand, the results provided by the aggregated ILP model do not show the Pareto optimal solutions. To our knowledge, few studies deal with the multiobjective aspect of the car relocation problem in one-way carsharing system. Therefore, we decide to develop a multiobjective optimization algorithm that tackles this problem from this perspective. In the following, we present a multiobjective formulation for the relocation problem in one-way carsharing system.

Objective functions definition
As described earlier, carsharing operators recruit employees to relocate cars between stations in order to reduce the number of rejected demands for taking a car or returning it to a station. Each car relocation operation op(i, j) starts when a jockey goes from his current location to station i where he will take a car; then, it ends at a destination station j where the jockey will deliver that car, with i 6 ¼ j. The duration of each relocation operation depends on the time needed to move between the stations. For this sake, we have a time matrix that contains the time needed to move between each pair of stations during the day based on the distance between these stations. For a carsharing system that has N stations we have N Â ðN À 1Þ relocation possibilities (since moving a car from a station i to the same station i does not make sense). Then, we put all the possible combinations of operations in one array that we name OP ¼ ð. . .; op k ði; jÞ; . . .Þ. Each element of OP has an index k that refers to one relocation operation, k 2 ½1; N Â ðN À 1Þ. Relocation operations can take place at different times of the day.
In this section, we present a multiobjective optimization algorithm where we consider three objectives to minimize: 2.3.1 Minimize the total number of remaining rejected demands: f 1 As described earlier, the relocation operations are necessary to reduce the number of rejected demands. The impact of each relocation operation a k at time t can be measured by the number of reduced or generated rejected demands D a k t , D a k t 2 ½À 2; þ 2. In the best case, D a k t can reduce two rejected demands (À 2 rejected demands). However, D a kt can generate two rejected demands (þ 2 rejected demands) when we make a bad relocation decision. A relocation operation consists of moving a car from a station i and relocates it to another station j. If for example the employee driver moved a car from a station i that had just one car, and few seconds later, a client arrives at this station to rent a car, he will generate a new rejected demand (þ 1). Within the same example, if the employee driver moved a car to station j that had n À 1 cars where n is the maximum number of cars a station can hold. And few seconds later, a client has arrived to this station to return his car, and he will not find any available parking space and thus another rejected demand will occur (þ 1). So that relocation operation will generate two rejected demands (þ 1 þ 1 ¼ þ2). On the other side, if the employee driver moved a car from a station i that had n cars, and few seconds later a client arrived to return his car, he will find an available parking space, thanks to the employed driver because he reduced one rejected demand (À 1). In this same relocation operation, the employed driver moved the car to station j that had zero cars, and few seconds later, a client arrived to rent a car, the client will find a car to drive, thanks also to the employed driver because he also reduced one more rejected demand (À 1). So this relocation operation reduced 2 rejected demands (À 1 À 1 ¼ À2). And sometimes the relocation operation will not have any direct effect on the rejected demands, so it will reduce 0 rejected demands.
To calculate the total number of remaining rejected demands f 1 , we do the summation of the impact of each relocation operation D a k t in each list of relocation operations L t , then we add the initial number of rejected demands R i Since the objective of the relocation operations is to reduce the number of rejected demands, we make sure to remove any operation that generates additional rejected demands. Thus, we only keep operations that reduce the value of f 1 .
To relate to the ILP model described in Sect. 2.2, f 1 can be expressed as follows: To calculate the number of jockeys needed to reduce the rejected demands, we create an array of integers U to count the number of jockeys involved in relocation operations during each time step. We refer to this array by U ¼ u 0 ;:::; u t ;:::; u T . As aforementioned, a jockey starts a relocation operation at the departure time t and finishes it at the arrival time to the destination station at time t þ t ij where t ij is the time needed to move from station i to station j. For each accomplished relocation operation, we increase the correspondent u t values in the range ½t; t þ t ij .
For example, if a jockey starts a relocation operation at time t ¼ 1 and he needs two time steps to arrive to the destination station (t ij ), then we should increment the correspondent two elements of the array U. However, for reasons of simplicity we also increment one more element to count the time needed to move from the destination station of the current operation to the next departure station of the next operation (we fixed this value to one). When we get the updated vector U of needed jockeys at each time step, we can get the value of the second objective function f 2 by calculating the maximum value in the array U: To relate to the ILP model described in Sect. 2.2, f 2 can be expressed as follows: 2.3.3 Minimize the total working time of jockeys: f 3 We propose the algorithm 1 to calculate the total operation time of jockeys using the array U described earlier. It starts by the initialization of the used variables. Then, in each iteration, the algorithm calculates the operation time for a jockey and adds it to the total time. We consider that the working time of a jockey j starts when he begins his first operation t u t min of the day and ends at the arrival time of his last performed operation t u tmax . In the proposed algorithm 1, we start by looking for the first index of the first and the last operation of the day in the U vector that contains the number of operations that should be done at each time step. The difference between these two values represents the working time for a jockey. After that, we decrement the elements in U vector by 1 if u t [ 0, starting from the index of the first operation until the index of the last operation.
To calculate the working time w j of a jockey j, we use the formula: We calculate the total working time of jockeys by doing the summation of all jockeys working time: To relate to the ILP model described in Sect. 2.2, f 3 can be expressed as follows: Algorithm 1 Calculation of the total working time of jockeys Require: U /*U is a vector that contains the number of jockeys needed for the relocation operations during each time step */ 1: T ime ← 0 /*Set total time to 0 */ 2: t ut min = 0 /*Set initial value for t ut min */ 3: while t ut min = −1 do 11: end for 12: end while 13:

return T ime
In the next two sections, we propose two metaheuristics to solve the multiobjective car relocation problem: the classical NSGA-II and a hybrid memetic algorithm as a hybrid approach. We start by NSGA-II, which will be our reference algorithm to compare the memetic algorithm that will be described after.

Resolution methods
In the following, we provide an overview of the methods and approaches used to solve and evaluate the multiobjective optimization problems. After that, we propose two algorithms for solving this problem: the well-known genetic approach NSGA-II and MARPOCS which is a hybrid approach combining genetic and local search approaches.

Overview of the resolution methods of multiobjective optimization problems
There are multiple approaches proposed in the literature for solving multiobjective problems (MOP). These methods can be classified into three main categories: • Approaches that are based on the transformation of a MOP into a mono-objective problem. Here we can find methods that are based on the aggregation of the objective functions f i into one function f. Despite the fact that these methods are easy to implement, they require a good understanding of the studied problem in order to choose the desired values for the coefficient or weight of each objective function. • Non-Pareto Approaches. In this category, all the objectives are treated separately. The main drawback here is that these approaches find only supported solutions that belong to the convex envelop. • Pareto approaches: this class of algorithms treats all the objective functions at the same time, taking advantage of the non-dominance concept.
Many studies in the literature prove that Pareto approaches are generally more efficient to solve multiobjective problems. As in mono-objective, hybrid approaches which are based on different simple heuristics, can be also applied in MOP. In many cases, hybrid approaches outperform the pure metaheuristics, as in [20][21][22]. The two proposed algorithms in the present paper for solving the car relocation problem belong to the Pareto approaches.
These approaches use the dominance relation in the algorithm process. Thus, solutions are not evaluated toward just one objective function; instead, they are evaluated with trade-off solutions that imply all the objective functions. Many Pareto algorithms have arisen after the introduction of the concept of the rank dominance of a solution in the first paper that discussed the Pareto [23]. Since, Pareto concept implies the use of populations, we find many genetic algorithms (GA) that belong to this family. However, Pareto approach also includes local search algorithms that are based on the concept of the Pareto, as well as hybrid algorithms. More recently, indicators-based methods arise in this category. Instead of evaluating the solutions of the population, we evaluate the population itself based on the chosen indicators. In the following, we present a brief description for the ranking operators.

How multiobjective solutions are compared?
Since a multiobjective solution cannot be evaluated against its objective functions as we do in mono-objective problems, it is important to define a metric capable of representing the quality of a solution. Thus, several methods have been proposed and are often associated with the commonly used GA. The main three methods that are listed in [24] are presented below: 1. Dominance rank For a given solution s, the dominance rank is defined by the number of solutions that dominate s, incremented by 1 [25]. Thereby, solutions of rank 1 are non-dominated solutions. As much as the number of solutions that dominate a solution s increases, its rank is worsened. The dominance rank was originally applied in the MOGA method (multiobjective genetic algorithm) [ , non-dominated solutions in the population are assigned the value 1, and then they are removed from the population. The new nondominated solutions are assigned the value 2, and then they are removed from the population as well. And so on, until that all the solutions of the population are assigned a value. In order to distinguish solutions that belong to the same front, we use complementary mechanisms, as we will see in the next sections. This definition of the dominance depth is used in NSGA algorithms (non-dominated sorting genetic algorithm) [26] and NSGA-II [27] that will be used in this study.

NSGA-II for the car relocation problem
Although there are many genetic algorithms that can be used in a multiobjective context, NSGA-II is considered a very popular non-dominated sorting genetic algorithm for solving MOP. It is an improved version of NSGA that integrates elitism and with no sharing parameter needed to be chosen a priori. It has the advantage of simplicity and good performance. This algorithm is often used as a reference to compare the performance of new approaches. Algorithm 2 recalls the pseudo code for NSGA-II in the general case.
Algorithm 2 Pseudocode for NSGA-II Require: Population size N , number of generations nbGen /*a time limit can replace nbGen */ 1: P ← init(N ) /*Initialize the population P with N random individuals */ 2: Q ← ∅ /*Initialize an empty population for children */ 3: eval(P ) /*eval is used for the evaluation of each individual */ 4: for i=1 to nbGen do 5: P ← P ∪ Q 6: assignRank(P ) /* Based on Pareto dominance */ 7: for each front f ∈ P do 8: setCrowdingDist(f ) 9: end for 10: sort(P ) /* By rank and in each rank by the crowding distance */ In the following, we will describe the genetic operators that we defined in the NSGA-II. Firstly, we present an individual coding, then we describe the selection, crossover and the mutation operators.

Individual coding
A good solution representation is a key factor in developing efficient evolutionary algorithms. The encoding of the solution should enable us to ease crossover and mutation operations, as well as perform easy fitness calculation. Here, a solution s for the car relocation problem consists of a list of relocation operations that are carried out by jockeys during all time steps of the day, s ¼ ðl 1 ; . . .; l T Þ. Each list of relocation operations l t contains the indices of the relocation operations that are applied at time t. The size of each list l t varies between 0 and R max where R max is the maximum number of rejected demands for the whole day before starting the relocation operations. A possible representation for car relocation plan may look like as shown in Table 1.
A chromosome is composed of T genes where each gene corresponds to a list of relocation operations l t . A relocation operation is represented by its index from the array OP which is the array of all possible combinations of relocation operations. For example, the line l 1 in Table 1 corresponds to relocation operations to do at time step 1, which are moving two cars from stations 1 to 3 and another car from stations 1 to 2.

Selection operators and crossover
In genetic algorithms, we use crossover operator in order to generate a new population by combining individuals from the current population between each other, in an analogy of the biological crossover. When we want to do the crossover, we have to select individuals from our population for the reproduction in the hope of obtaining new chromosome (offspring) that holds the best characteristics of both parents. There are many selection methods such as roulette wheel selection (RWS), elitism selection, rank selection, binary tournament, etc. In NSGA-II, a binary tournament is used: two individuals are chosen randomly from the population then are compared. The best dominant individual is selected to be the first parent. We do the same process for the second parent. Binary tournament is not a very selective method. In that way, NSGA-II avoids the premature convergence.
When we get the pair of parents, we apply the crossover operator. Unlike the selection operator, there is no particular crossover operator defined for NSGA-II. It depends on the studied problem. In our case, we propose a hybrid method for crossover, which is a combination between two points crossover and uniform crossover as shown in Fig. 1.
First, we start by randomly selecting two points for the crossover then we exchange the parts between these two points. After that, we apply the uniform crossover on the remaining parts, which consists of exchanging bits randomly between the two parents. In this way, the proposed crossover operator offers a balance between the uniform crossover that breaks the sequence of operations (high diversification) and the two points crossover that maintains a sequence bloc of operations (weak diversification).

Mutation operator
After the crossover, the mutation operator is used to maintain a genetic diversity in the population by changing one or more genes to new values that we may not find in the parents genes. Usually, mutations occur according to a low mutation probability called mutation rate. In general, we have different methods to apply a mutation depending on the encoding of the solution. For example, in binary encoding, we can apply a simple mutation by inverting the selected bits as shown in Fig. 2. ...

l TÀ1
(2, 6) l T In the proposed model, a special mutation is adapted to the problem representation. For this sake, we choose two random genes from the selected chromosome then, for each selected gene, we apply our dedicated mutation operator, which consists of adding an operation if the gene does not have any operation; otherwise, we randomly select an operation that belongs to this gene and we replace it by another operation selected randomly. When we add or replace a relocation operation, we check if the new operation causes new rejected demands. In this case, we remove that operation which deteriorates each fitness value. That way, the proposed mutation operation can add, replace or remove a relocation operation. We propose to choose two genes for the mutation operator in order to increase the diversification effect, since we find that mutating just one operation is not enough regarding the total number of operations.

Memetic algorithm
A memetic algorithm (MA) was introduced for the first time by Moscato et al. [28]. According to its definition, a MA consisted of modified version of a genetic algorithm where a local search is integrated. Hybrid approaches have proved their efficiency in solving single objective optimization problems. As a result, researchers were interested in testing the effect of metaheuristics hybridization on the solving MOP. In general, the mechanism consists of adding a local search in a genetic algorithm. The genetic algorithm plays the diversification role. It enables the algorithms to discover more zones in the research space. For every found solution, local search is used to improve its quality and therefore pushing it to the Pareto front.
In a multiobjective context, an algorithm has the role of finding the ''best'' possible set of non-dominated solutions. In general, Pareto set refers to the feasible non-dominated solutions that cover the entire search space. However, in real-life optimization problems, it is not often easy to obtain all the solutions that belong to the Pareto front. Nevertheless, in most cases we use Pareto front when we want to describe the best found non-dominated solutions.
The general scheme of a memetic algorithm is presented in ''Algorithm 3.'' Algorithm 3 General form of a memetic algorithm Require: N the size of the population, nbGen the number of generations /*A time limit or a any termination condition can be used to replace the nbGen */ 1: P 0 ← init(N ) /*Initialize the population P with N random individuals */ 2: P 0 ← ∅ /*An empty population for children */ 3: eval(P 0 ) /*Evaluate objectives for each individual */ 4: for i = 1 to nbGen do 5: Select P i from P i−1 based on fitness value 6: Apply genetic operator on P i → P i 7: Apply local search on P i → P i 8: eval(P i , P i ) 9: Select P i from (P i−1 ∪ P i ∪ P i ) 10: end for It is important to note that there are different methods to apply a local search in a GA. Therefore, local search can be used with a probabilistic factor just like a mutation, or when the algorithm reaches a given number of generations. Likewise, local search can be applied on every created individual. The choice of how to use the local search and how to maintain an equilibrium between the genetic operator and the local search is crucial to obtain good results.
When designing an algorithm for solving an optimization problem, it is important to study how it processes the search space based on the intensification and diversification aspects. In a genetic algorithm, the crossover and the selection operators are considered as intensification operators while the mutation is used for diversification. However, in a memetic algorithm, roles are exchanged between mutation and crossover: the mutation operator (LS) becomes an intensification operator while the crossover operator becomes a diversification operator. From this perspective, a memetic algorithm can be considered as an improved version of population-based local search where  Multiobjective car relocation problem in one-way carsharing system 305 the crossover operator is used for diversification and the local search is used to ameliorate the found solutions.

Genetic operator in memetic algorithm
There are many genetic approaches that are used in multiobjective context (multiobjective evolutionary algorithm, MOEA) that can be hybridized with local search approaches. SPEA and NSGA-II are frequently used for this type of hybridization [27,29], but all GA can be hybridized. A list of these approaches can be found in [30]. In general, MOEA constitutes a good base for applying a local search in population-based MOPs. In hybrid approaches, we still find the traditional genetic operators: selection and crossover. However, the mutation operator is often used to apply the local search.
In the present paper, the genetic algorithm that is used as a base for MARPOCS is NSGA-II [27], for its simplicity and its good performance. In the next section, we will describe the local search that will be used in MARPOCS.

Local search procedure to improve solutions
Dominance-based local search approaches [31] deliver very good results. These LS can be integrated in memetic algorithms. For example, M-PAES is a modified version of PAES, which is adapted to memetic algorithm [20]. However, adding a sophisticated local search metaheuristic to a genetic algorithm can make the memetic algorithm complicated in terms of implementation and execution. Hence, a simple local search algorithm combined with genetic operators can be more efficient in solving MOPs.
There are different approaches for applying a local search algorithm in memetic algorithm. Each objective can sequentially processed separately, or all the objectives can be treated together by random weighting to push the found solutions to Pareto front. When the algorithm processes solutions during local search procedure, it might encounter new non-dominated solutions. These new solutions can be simply copied to the archived set or they can be added to the current population.
To adapt the problem in the present paper to local search procedure, it is turned into a mono-objective problem where the fitness value is equal to the weighted sum of the different criteria with respect to the chosen search direction, which is randomly generated and defined by the weights vector x.
Each non-dominated solution found during the local search is added to the archive. In that way, the archive is built directly using the local search procedure. The last solution found in the local search procedure (local optimum) is added to the population used in genetic algorithm in order to be used to build the next generation.
The local search strategy used in this paper belongs to the first improvement hill climbing (FIHC) category. This type of local search successively accepts the first neighbor that has a better fitness. This way, a partial neighborhood exploration is done during each LS iteration. The algorithm stops when it cannot find any new improved solution, which is the local optimum. A local search is formally defined by the couple ðX; VÞ where X is the set of feasible solutions (decision space) and V is the neighborhood structure V : X ! 2 X that assigns for each solution s 2 X a set of neighbors V(s). Therefore, the algorithm stops when the algorithm finds the local optimum s Ã , that is reached if and only if 8 s 2 Vðs Ã Þ; f ðsÞ ! f ðs Ã Þ.
In MARPOCS, the neighborhood of a solution s is defined as the set of solutions s 0 that can be obtained by adding or removing an operation from the solution. In a more formal way, this relation can be expressed as: where dðs; s 0 Þ is the Hamming distance between s and s 0 . When all the operations in s and s 0 are the same then dðs; s 0 Þ ¼ 0. While dðs; s 0 Þ ¼ 1 when there is one operation in difference between s and s 0 . We also note that changing the time of an operation goes into two steps: removing the operation then adding it in a different time step. In MARPOCS, the local search procedure is performed through three steps as shown in ''Algorithm 4'': In the first step, the algorithm uses a list of the relocation operations that can decrease the number of rejected demands, obtained using the algorithm described in [19]. The algorithm calculates the weighted fitness of the solution after adding each operation from this list until finding an improvement. The algorithm will keep any operation that improves the fitness of the solution while it disregards the others.
Then in the second step, the local search procedure calculates the weighted fitness after removing each operation of the solution. If the fitness is not improved, the algorithm restores that operation.
In addition, in the third step the algorithm evaluates the fitness of the solution after changing the time of each operation (forward or backward through the time). If the fitness is not improved, the algorithm restores the old value of time.
Every time the fitness is improved, the solution is added to the archive if it is not dominated by any other solutions while all dominated solutions are removed. However, just one solution will be modified in the genetic population.

Algorithm 4 Local search procedure for the proposed memetic algorithm
Require: s the offspring that we want to improve /* s the child that comes after a crossover operation */ 1: /* Start Step 1 */ 2: O ← getGreedyOperations(s) /* Get the operations that reduce the number of rejected demands using the policy 3*/ 3: Ω ← generateRandomWeights() /* ω the vector of random weights of objective functions, such as n i=1 s.add(op) /*Integrate the operation in the solution */ 7:

Memetic algorithm process
In a memetic algorithm, there are three main steps: • Selection of parents p 1 and p 2 : • Crossover between the two parents that gives a new solution s. • Local search process based on the randomly generated vector that represents the search direction v to obtain the local optimum solution s Ã . This vector is defined by the weights vector x ¼ ðx 1 ; . . .; x k Þ. During the local search, the algorithm accepts all solutions s i obtained through the neighbor structure and that have better fitness based on the aggregation function.
The vector that defines the progress of the local search allows a good exploration of the search space. In fact, the vector v and the current solution s i , define a hyperplane orthogonal to v and passes by s i . During the local search, all the solutions s iþ1 that are above of this hyperplane improve the aggregated functions and are accepted. Figure 3 illustrates this LS for two objectives, in that case the hyperplane is a line. The algorithm below describes the steps used in MARPOCS:

Algorithm 5 Proposed memetic algorithm
Require: N the the size of the population, nbGen the number of generations /*A time limit or any termination condition can be used to replace the nbGen */ 1: P ← init(N ) /*Initialize the population P with N random individuals */ 2: Q ← ∅ /*An empty population for children */ 3: A ← ∅ /*Initialize an empty archive */ 4: eval(P ) /*Evaluate objectives for each individual */ 5: for i = 1 to nbGen do 6: P ← P ∪ Q 7: assignRank(P ) /*Based on Pareto dominance */ 8: for each front f ∈ P do 9: setCrowdingDist(f ) 10: end for 11: sort(P ) /*By rank and in each rank by the crowding distance */ 12: Q ← buildChildren(P ) 14: improveChildren(Q) /*Following the local search FIHC */ 15: end for The main difference between this algorithm and NSGA-II, is the integration of the new function improveChildren. This function applies the local search on all the individuals of the population Q. The selection and the crossover operators are applied in the function buildChildren. However, the mutation operator is omitted since it is replaced by the local search.

Mobility data
The mobility data used for this study comes from a real context that consists of survey data and socio-economical information collected by professional for regional planning purposes. It describes people mobility flows in a region of 20 km Â 10 km in Paris. This area is divided into a grid of equal cell sizes. Each cell has two characteristics: • Terrain type: to describe the dominant structure type of the area associated to the cell (roads, buildings, houses, business center, commercial center, etc.) • Attraction weight: based on the terrain type and survey data, this information attributes a dynamic attraction weight to each cell for each 15 min of the day People mobility between different cells is represented by a 3D matrix F = ðf i;j;t Þ, where f i;j;t represents the number of people who want to move from cell i to cell j at time t. We consider t to be a period of 15 min during the day, which makes 96 time periods. Then the flow mobility data is plotted on a map using GIS shapefiles. As a result, 400 cells have been detected as a potential origin or destination point knowing that some cells are eliminated because of their geographical nature, e.g., lakes, plains, etc. The final flow mobility data consists of 400 Â 400 Â 96 elements, which makes 15,360,000 records to represent how people move during the day.
In [32], we developed a platform for locating stations for a carsharing system. This platform takes advantage of the mobility data that we described earlier in Sect. 4.1. We used this platform to generate the data for this study. For each generated dataset, we use four parameters: 1. Total number of cars in the system. The platform is used to find the optimal locations of stations in a new carsharing system within an urban area. Given a combination of the described parameters, this platform locates the stations in a way that covers the maximum user demands while respecting the predefined objectives. The platform is used to constitute the data that will be used in the optimization process. It is important to mention that the traffic condition is omitted in this study, because the used platform does not provide precise information about it. It is really very important to include this variable when implementing the study, and it can be estimated statistically in case of real system data.

JMetal framework
JMetal is an open source object oriented framework based on Java for multiobjective optimization using metaheuristics [33]. This framework is provided with many state-ofthe-art algorithms, benchmark problems, and different quality indicators that are used to evaluate the performance of MOP. We use this framework to develop, experiment and to study our algorithms in order to solve the multiobjective optimization carsharing problem. The framework also facilitates the development of experimental studies with special tools that are dedicated for this purpose. It also offers the feature of statistical analysis of the experimentation results. In addition, this framework allows running experiments in parallel to exploit the multi-core processors in order to reduce the execution time of the used algorithms.
In our study, we use JMetal 4.5 on a PC with Intel Core i5-3550 CPU (3.30 GHz) and 16 GB of RAM.

Evaluating the quality of a set of solutions
In mono-objective optimization, the evaluation of algorithms is easy since solutions have only one fitness value. However, this is not the case in MOP where each solution has k values for k objectives, and each algorithm generates trade-off solutions that constitute an approximation of the Pareto Front. Thus to evaluate two algorithms, we have to compare two sets of solutions generated by each algorithm. When all the solutions generated by an algorithm B are dominated by a solution of algorithm A, in this special case we can say that A is better than B. However, it is not easy to compare two set of solutions where some solutions of A dominate some solution of B and vice versa. Therefore, we generally use performance indicators that can be used to evaluate an approximation of Pareto front. These performance indicators evaluate the solutions based on two criteria: convergence and diversity. There are many quality indicators that are proposed in the literature such as hypervolume, generational distance, inverted generational distance, epsilon, spread, etc. Each quality indicator is intended to measure one or both of the criteria at the same time as shown in Table 2. In another perspective, quality indicators can be classified into two main categories: unary indicators that assign a score for a set of solutions, and binary indicators that compare two sets of solutions and assign a score to a set relatively to the other one.

Epsilon
The epsilon indicator has been presented in [34]. This indicator has two versions: multiplicative and additive. In its additive version, the epsilon indicator is used to Fig. 3 Selection and crossover followed by a local search calculate the minimum translation factor e that should be added to every solution of an approximation A in order to dominate solutions of the Pareto front. This indicator can be defined as follows: where z 1 0 e z 2 if and only if 81 i n : The smaller the value of e is, the better the approximation is. It is good to note that a normalization of the objective functions is required in order to obtain efficient results when the scale is not the same for all the objective functions.

Generational distance (GD)
The generational distance indicator was proposed by [35]. This indicator measures the average distance (in objective space) between the approximation set and the Pareto front as shown in the Eq. 11. The smaller the value of GD is, the better the approximation is. When all the approximation solutions are found in the Pareto front the value of GD is 0.
with n being the number of solutions contained in the approximation set and d i the Euclidean distance in the objective space between each solution of the approximation set and its nearest solution in the Pareto front.

Inverted generated distance (IGD)
IGD is a modified version of GD where this indicator calculates the average distance between the Pareto front and the approximation set instead of the average distance between the approximation set and the Pareto front.
where n 0 is the number of solutions contained in the Pareto front and d i the Euclidean distance in the objective space between each solution of the Pareto front and its nearest solution in the approximation set.

Spread (D)
The spread indicator measures the spread level of the solutions that belong to an approximation set according to the Pareto front. It is defined as below: where d i is the distance between two adjacent solutions, d is the average value of these distances, and d f and d l are the Euclidean distances between the extreme solutions of the Pareto front and the approximation set. A spread value of zero refers to an ideal distribution (Fig. 4).

Hypervolume (HV)
The hypervolume indicator computes the volume in the objective space of the space covered by the set of nondominated solutions A of an approximation set. This indicator is very popular in the literature since it reflects the quality of an approximation based on two criteria: diversity and convergence. In its unary form, the hypervolume is calculated based on a reference point Z ref (''anti-optimal'' or nadir point) which refers to the worst possible point in the objective space as shown in Fig. 5 where we deal with a minimization problem. A hypercube v i is formed of each solution s i 2 A and the reference point Z ref , as diagonal corners. The hypervolume is the union of all hypercubes:

Reference front
We use quality indicators to compare the performance of algorithms in MOPs. In most cases, we compare the approximation set with the Pareto front. However, in realworld problems, usually the Pareto front is unknown. To cope with this problem, we build an approximation of the Pareto by combining all the results of all the runs tested for our problem. Thereby we can compare the relative performance of the algorithms. This is our approach for the car relocation problem.

Algorithms parameters
In the present paper, we compare two algorithms for solving the multiobjective carsharing problem: NSGA-II and memetic. We have to choose the most suitable values for the genetic parameters such as mutation probability, crossover probability and population size. A good choice of these parameters highly affect the performance of the algorithm. For this sake, we define different combinations of these three parameters then we use a class that is defined in the JMetal platform for this purpose. The concept is to run the algorithms many times with different parameters combinations. After that, we compare the results using the values of performance indicators that are defined in the platform and precisely the hypervolume (HV). The configuration that is associated with the best performance indicators values is considered to have the best parameters settings. Since the proposed algorithms are stochastic, we run our algorithms 5 times for each combination before comparing the results. After that, we compare the results of each algorithm for each configuration based on HV.
The values of the parameters that are tested together for NSGA-II and MA are shown in Table 3.
Therefore we test 10 levels of population, 5 levels of crossover and 6 levels of mutation probabilities. That is 300 combinations (10 Â 5 Â 6) tested for NSGA-II algorithm and 50 combinations (10 Â 5) for the memetic algorithm, since the local search replaces the mutation for all the individuals and is applied in a systematic way to offspring.
We start by running the NSGA-II algorithm over the 300 combinations with 5 runs for each. The average and the standard variation of the hypervolume are used to rank the performance of these different parameters combinations based on the hypervolume indicator. The best 5 combinations are shown in Table 4.
On the other side, we run the 50 experiments of the memetic algorithm with 5 runs for each combination parameters. Here also, we rank the experiments based of the HV indicator. The best 5 combinations are shown in Table 5.

Comparison of algorithm results
In order to evaluate the performance of the multiobjective approach for the car relocation problem, we use the best combination of parameters obtained for each algorithm in Sect. 4.5. The results concern a carsharing system that is composed of 18 stations with 10 parking spaces each, having a fleet of 88 cars with an average of 12 trips per car during the whole day. In the same way as for fixing the parameters, each algorithm is run 5 times using the best combination of parameters found in Tables 4 and 5.
We compare the approximation sets obtained using both algorithms based on the best combination chosen for each algorithm in Sect. 4.5. Since it is not easy to compare solutions of multiobjective problems in a 3-D space as shown in Fig. 6, we propose to plot each two objectives on a 2-D chart in order to make the comparison more readable.
In Fig. 7, we see a two objectives comparison between the solutions generated by NSGA-II and those generated by the MARPOCS. The number of rejected demands is plotted against the number of jockeys. As we can see, most of the solutions obtained with the MA dominate the solutions obtained using NSGA-II algorithm. The initial total number of rejected demand without any relocation operation is 312. As shown, the same number of jockeys may be associated with different numbers of rejected demands. Obviously, this can be explained by the fact that the third objective is not shown. The lowest values of remaining rejected demands for each number of jockeys are associated with the highest values of working time for this number of jockeys. The comparison also shows that MA succeeds to generate solutions that cover the best values for the rejected demands objective. Starting the thirteenth jockey, the algorithm starts to generate solutions with zero rejected demands, which is consistent with the value of the number of jockeys required to solve all the rejected demands in [19] using an exact or a greedy approach.
In Fig. 8, the working time is plotted against the number of jockeys. Here also, the number of jockeys can be associated with different values of working time depending on the value of the third objective. This proves the ability of the algorithm to generate different solutions for the same number of rejected demands depending on the dedicated working time. For example, in Fig. 7, the highest values of remaining rejected demands ( [ 130 rejected demands for number of jockeys [ 21) are associated with the lowest values of total working time shown in Fig. 8 (total working time \ 100 for the number of jockeys [ 21). This is the interest of a multiobjective approach: Giving the possibility to the operators to choose between several equivalent solutions. Sometimes, some solutions are not interesting. For example, solutions with high number of jockeys, high rejected demands and low total working time, may not be useful for the operator. But those solutions will be listed by a multiobjective algorithm since they are not dominated on at least one criteria.
It is important to notice that the memetic approach, as shown in Figs. 7 and 8, failed to obtain solutions with low remaining rejected demands when the number of jockeys exceeds 21. Firstly, this is not a problem, since the algorithm was able to find a set of solutions with 0 remaining rejected demand by using a lower number of jockeys. Secondly, this can be explained by the fact that the algorithm did not have the needed time to find those solutions.
In the last combination of objectives, we plot the working time against the number of remaining rejected demands (Fig. 9). As we can observe, the solutions quality of MA is much better than those of the NSGA-II algorithm regarding the total working time objective. The MA can generate solutions that solve all the rejected demands using a lower value of working time. It is also important to note that even when NSGA-II uses a greater working time, it cannot go under 31 rejected demands.
From another point of view, In Fig. 9, there is a negative relationship between the rejected demand and total working time. This can be explained by the fact that if we have rejected demands, we need to do relocation operations in order to avoid them. Each relocation operation needs time to be done. For example, we initially have 312 rejected demands in the system; the total working time needed to keep them as they are is zero, and as much as we do relocation operations to decrease the number of rejected demands, the total working time will increase also. In this figure, we see that MARPOCS performs better than NSGA-II. For keeping the same number of rejected demands, MARPOCS finds always a solution with less total working time than the solution found by NSGA-II.
From indicators point of view, Fig. 10 shows a comparison of the evolution of the hypervolume indicator and the number of solutions in the archive for NSGA-II and MA. The hypervolume indicator is obtained by comparing the approximation obtained for each algorithm with the reference set obtained earlier by combining the results of all the runs of our algorithms using different parameters settings.
As we can see, the hypervolume of solutions obtained by the MA is always higher than the hypervolume of solutions obtained by NSGA-II. We can also notice that after the first two seconds only, the hypervolume of MA solutions has the same level of the hypervolume of solutions obtained by the NSGA-II solutions after 30 s of algorithm run. After that, the hypervolume related to the MA starts to stagnate after 40 s when it reaches the 0.6 level. While the hypervolume of the NSGA-II algorithm    Multiobjective car relocation problem in one-way carsharing system 311 starts to stagnate after the 125 s of execution when it reaches the 0.52 level.
On the other side, we notice that the evolution of the number of solutions in both algorithms has the same trend over the time. It is also clear that sometimes the number of solutions in the archive of the NSGA-II algorithm is greater than the number of solutions in MA. However, these solutions have of lower quality as indicated by the hypervolume and as will be shown in the next paragraphs.
In Table 6, we observe a comparison between NSGA-II and MA based on different quality indicators. The comparison shows that the MA performs much better than the NSGA-II algorithm with regard to almost all the quality indicators. The number of solutions in the NSGA-II archive is close to that number in MA archive. However, the quality of solutions in MA is better than those of NSGA-II as shown in the quality indicators values. The spread value in NSGA-II is sometimes greater than its value in MA. However, all other quality indicators show that the quality of the MA solutions is better than those of NSGA-II solutions.
It is also worth noting that we run the algorithms through different carsharing configurations, and they proved that in average they needs between few seconds to one minute to converge in case of medium networks (30 stations). This duration should be acceptable to meet the

Conclusion
In this paper, a new approach is proposed for addressing the one-way car relocation problem based on a multiobjective problem formulation. The proposed model allows to find best trade-off solutions without making assumptions about the cost of a rejected demand, which is not easy to determine. Three objectives are fixed and modeled for the presented problem. To solve this problem, two algorithms are implemented. The first one is a classic NSGA-II that has been adapted to the problem and used as a reference. The second one is a new memetic approach.
The comparison of the approximation sets obtained by both algorithms shows how much the solutions generated by the MARPOCS are better than the solutions generated by NSGA-II. This observation is proved by the comparison of different quality indicators values that are used to compare the performance of each algorithm.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.