Aurélien Carlier
email: aurelien.carlier@renault.com

Alix Munier Kordon

Witold Klaudel
email: witold.klaudel@renault.com

Optimization of a one-way carsharing system with relocation operations

Keywords: carsharing, operational research, carsharing demand generation, transport optimization, mathematical modelling

Carsharing systems have gathered increasing attention and are now recognised as innovative and ecological solutions to transportation issues. Especially, one-way carsharing systems offer a high level of service because they exempt users from the obligation to return the vehicle at the same station where it was borrowed. Unfortunately, this flexibility comes with design complexity. This paper addresses a part of those design issues and deals with fleet dimensioning including vehicle relocation operations.

We propose mathematical programming oriented approach and we introduce a simple linear model based on integer flow variables. This model is organized around three optimization criteria: maximizing satisfied carsharing demands while minimizing the fleet of vehicles and the relocation operations. From random generated benchmarks, we study the scalability of our method and we show that computation times remain acceptable for representative problem complexities.

I. INTRODUCTION

In the last decade, the importance of adapting the way we move in our societies has been highlighted by a lot of researches. The greater accessibility to private car has more particularly resulted in serious negative externalities, although it significantly improved people mobility in urban areas. Pollution, excessive energy consumption and wasted time evolve largely from the road network congestion. In many situations, current public transport does not offer a valuable alternative to the private car and a lot of efforts are made to find a different and acceptable solution [START_REF] Mitchell | Reinventing the automobile: Personal urban mobility for the 21st century[END_REF]. One of them is carsharing.

A carsharing system involves a small to medium fleet of vehicles, available at several stations, to be used by a relatively large group of members [START_REF] Shaheen | A short history of carsharing in the 90's[END_REF]. This service constitutes a real alternative to private car because it releases the user from constraints related to individual property like insurance, maintenance, fuel (or electricity), taxes, vehicle depreciation, etc. Such systems became popular in the last couple of decades and nowadays it is not unusual to find one in most of the world's major cities in industrialized countries. Their numerous benefits fostered several researches on environmental impacts, social characterization of the users, understanding of urban trip patterns, demand modelling, vehicle imbalance problem, operator relocation operations policy, optimal system dimensioning, etc. Jorge and Correia [START_REF] Jorge | Carsharing systems demand estimation and defined operations: a literature review[END_REF] give a review of these topics.

Initially, in the so called round-trip carsharing systems, the users have to pick up and return the vehicles at the same place. From an operational point of view, this limitation simplifies greatly the dimensioning and the management of the vehicles fleet. In this paper, we will focus on one-way carsharing systems which allows users to pick up a vehicle from a station and return it at a different one. The greater flexibility of these systems helps capture more requested trips (referred further as demands) [START_REF] Firnkorn | What will be the environmental effects of new free-floating car-sharing systems? the case of car2go in ulm[END_REF] but comes also with serious operational issues, especially those from unbalanced vehicle distribution. As a consequence, the dimensioning of a one-way carsharing system must include vehicle relocations to capture and meet real demand as much as possible. Relocation is often not considered in car rental services [START_REF] Fink | Modeling and solving the short-term car rental logistics problem[END_REF], [START_REF] Yang | Car rental logistics problem: A review of literature[END_REF] since the price of one-way trips is usually fully supported by the customers.

Some studies have been conducted to find optimal vehicle relocations when a system already exists (e.g. [START_REF] Wang | Dynamic relocating vehicle resources using a microscopic traffic simulation model for carsharing services[END_REF], [START_REF] Cucu | Management of a public transportation service: Car-sharing service[END_REF], [START_REF] Nair | Fleet management for vehicle sharing operations[END_REF], [START_REF] Smith | Rebalancing the rebalancers: Optimally routing vehicles and drivers in mobility-on-demand systems[END_REF]), while others searched for optimal dimensioning of the vehicle fleet taking into account those relocation operations. In most cases, mathematical programming is considered as a good and pertinent approach to describe and solve those problems. Fan et al. [START_REF] Fan | Carsharing: Dynamic decision-making problem for vehicle allocation[END_REF] proposed for example a multistage stochastic linear integer model accounting for demand variations and indicating operator profit improvement. More recently, Correia and Antunes [START_REF] Correia | Optimization approach to depot location and trip selection in one-way carsharing systems[END_REF] developed three mathematical programming models to balance vehicle stocks while dimensioning the whole system in terms of location, number and size of stations. The objective was again maximizing the profit for the operator considering all the revenues and costs involved.

Although mathematical programming is an effective and suitable approach to deal with the optimization of carsharing systems, the studies mentioned above often underline limitations due to computation time and solver capability. This frequently led to model simplifications which makes results unrealistic. Thus, our first aim is to develop a simple realistic and scalable flow based model taking into account relocation operations. The scalability is confirmed by the acceptable computation times for the large instances we analysed.

Three optimization criteria are considered: maximizing the total number of satisfied demands and minimizing the total number of vehicles and relocation operations. We experimentally proved that the two last criteria are in opposite and that there exists a tradeoff area between them.

To the best of our knowledge, there is no available benchmark for testing our model. Most authors use real data given by a carsharing operator or extrapolated from real life observations. Hence, we developed a random data generator taking into account the demand variation over the course of a day.

This paper is organised as follow. Section II formulates on the one hand the problem in mathematical terms, and proves on the other hand that the number of vehicles can be inferred from any solution. Section III is dedicated to numerical experiments based on randomly generated data. Results specifically focus on computation times, 3-dimensional Pareto frontiers and scalability. Finally, section IV concludes the paper with directions for further research.

II. PROBLEM MODELING

This section aims at modelling our optimization problem using an integer linear program. The inputs and the outputs of the problem are first described. Second subsection is dedicated to the building of an oriented valued graph, namely the Timed Extended Graph. This graph, previously introduced by Ahuja et al. [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF], allows to express all the constraints of the problem following the time and the space dimensions. The third subsection introduces the decision variables of our optimization problem. More precisely, it is shown that vehicles can be equivalently aggregated into flows to express all the constraints and the criteria of our optimization problem. Last subsection presents its formulation using integer linear programming.

A. Inputs and outputs of the problem

In order to describe the mathematical problem, let us define first the required inputs. As discussed later, those data can be extracted from various sources such as simulation tools or real operating carsharing system data for example.

First of all, let H = {1, • • • , T } be the set of time-steps considered in the study. For relevance issues, it has to cover a representative period of time, as an average weekday or an average week for instance.

The set of carsharing stations is defined as N = {1, • • • , N } and comes with Z(i) the maximum size (i.e. maximum capacity, in terms of number of vehicles) of station i ∈ N . Then, the demand D(i, j, t) contains the number of passengers wishing to borrow a car from station i ∈ N at time t ∈ H to station j ∈ N . Note here that there is no condition on the station themselves, thus allowing the travel to be a "round trip" or a "one-way" travel. In this case, i is simply equal to j.

Finally, δ(i, j, t) represents the travel time it takes for a car-user from station i ∈ N to station j ∈ N , when departure time from i is t ∈ H. We suppose that for any triple (i, j, t) ∈ N × N × H, δ(i, j, t) < T . This assumption comes from the fact that the highest distance from different stations is quite low (usually less than 150 or 200 kilometres) and that the time-steps are covering at least a day.

A feasible solution of our problem consists on a set of vehicle tours, each of them modelling the situation of a car at each time step. The three criteria (the demand, the number of vehicles and the number of relocation operations) can be polynomially computed from any feasible solution.

B. Time Extended Graph

To deal with discrete-time dynamic networks, Ahuja et al. [START_REF] Ahuja | Network flows: theory, algorithms, and applications[END_REF] suggested the use of time-space network, also known as Time Extended Graphs (TEG). It is a static network constructed by expanding the original network in the time dimension, considering a separate copy of every node i ∈ N at every discrete time-step t ∈ H. Thus, from the data listed above, let us consider G = (X , A, u) as a TEG defined as follows:

1) Set of nodes: Nodes are couples (i, t) with i ∈ N and t ∈ H associated with station i at time t. Formally, X = N × H.

Let η and θ the functions which reciprocally return the station and the time-step associated with every element of X . More formally,

η : X → N with x = (i, t) → η(x) = i θ : X → H with x = (i, t) → θ(x) = t
2) Set of arcs: Edges are partitioned into 3 sets A 1 , A 2 and A 3 defined as follows:

• A 1 is the set of arcs representing the possibility for the vehicles to stay at a station between two consecutive time-steps. Formally,

A 1 = {(x, y) ∈ X × X | η(x) = η(y) and θ(y) = θ(x) + 1 mod T }. • A 2 are
the arcs associated with a demand: each demand D(i, j, t) > 0 corresponds to an arc (x, y) with x = (i, t) and y = (j, t + δ(i, j, t)). Set A 2 is then formally defined as

A 2 = {(x, y) ∈ X × X | D(η(x), η(y), θ(x)) = 0 and θ(x) + δ(η(x), η(y), θ(x)) = θ(y) mod T }. • A 3 represents
all the possible relocation operations over time. Any triple (i, j, t) ∈ N × N × H with i = j is associated to an arc from x = (i, t) to y = (j, t + δ(i, j, t) mod T ). Thus,

A 3 = {(x, y) ∈ X × X | η(x) = η(y) and θ(x) + δ(η(x), η(y), θ(x)) = θ(y) mod T }.
The total number of arcs is given by:

|A| = |A 1 | + |A 2 | + |A 3 | = N × T + M + (N × T ) • (N -1)
where M is the number of requested demands. As M N 2 , |A| = Θ(N 2 • T ).

3) Arcs Capacity: Associated with each arc a = (x, y) ∈ A, is given a capacity function u : A → N corresponding to a maximum number of vehicles allowed on a. It is defined as follows for any arc a ∈ A:

u(a = (x, y)) =    Z(η(x)) if a ∈ A 1 D(η(x), η(y), θ(x)) if a ∈ A 2 +∞ if a ∈ A 3
For any arc a = (x, y) ∈ A 1 , the maximum number of cars is the capacity of the station η(x) = η(y). It corresponds to the demand for any arc a ∈ A 2 , and it is not bounded for relocation arcs a ∈ A 3 .

4) Additional notations:

We denote by Γ -(x) and Γ + (x) respectively the set of immediate predecessors and successors of a node x ∈ X , i.e.

Γ -(x) = {y ∈ X | (y, x) ∈ A} Γ + (x) = {y ∈ X | (x, y) ∈ A}
For any couple of time instants ∀(t, t ) ∈ H 2 , the number of time-steps between those two instants is defined through the following function

ϑ : H 2 → N (t, t ) → ϑ(t, t ) with ϑ(t, t ) = t -t if t ≤ t T + t -t otherwise.
For each arc a = (x, y) ∈ A, let us define the boolean value a as:

a = 0 if θ(x) ≤ θ(y) 1
otherwise The time required for a movement from x to y is then given by the function

: A → N a = (x, y) → (a) = θ(y) -θ(x) + a • T By extension, if µ = (a 1 , • • • , a p ) ∈ A p
is a path of the TEG from x to y, the value (µ) = p i=1 (a i ) is the total time required for a vehicle going from x to y following µ.

For any time value t ∈ H, let us define the set C t (µ) as the arcs a = (x, y) from µ starting at time t or earlier but ending after t. Formally,

C t (µ) = {a = (x, y) ∈ µ | ϑ(θ(x), t) < (a)}.

C. Decision variables

The aim of our study is to compute the planning of each vehicles during the period. At any time, each vehicle is either parked in a station or in transit between two stations. Its position over the period can be modelled as a vehicle tour i.e. a circuit c = (a 1 , • • • , a p ) in the TEG.

The size of any feasible solution may be highly reduced if we only consider the number of vehicles passing through each arc. For each arc a = (x, y) ∈ A, we call ϕ(a) the flow of vehicles transiting through the arc a. It can be interpreted as:

• the number of vehicle staying in station η(x) between two consecutive time-steps θ(x) and θ(y), if a ∈ A 1 ; • the number of vehicle picked by users from station η(x) at time θ(x) to station η(y), if a ∈ A 2 ; • the number of vehicle relocated between stations η(x) and η(y) at time θ(x), if a ∈ A 3 . The total number of vehicles transiting to any node x ∈ X is clearly constant, thus

y∈Γ -(x) ϕ((y, x)) = y∈Γ + (x) ϕ((x, y)).
It is immediate that a feasible flow may be obtained from any feasible set of vehicle tours. We prove in the following that vehicle tours may be easily computed from a feasible flow.

Next lemmas characterize the total time of any circuit c and the exact number of vehicles required for a unitary flow on c.

Lemma 1: The total time of any circuit c = (a T . Lastly, the minimum number of vehicles to insure ϕ c is constant over H and thus ∀t ∈ H, the lemma.

1 , • • • , a p ) is (c) = T × p i=1 ai . Proof. Let x i , i ∈ {1, • • • , p + 1} be the sequence of elements from X such that, x p+1 = x 1 and ∀i ∈ {1, • • • , p}, a i = (x i , x i+1 ). The total time of c is then (c) = p i=1 (a i ) = p i=1 (θ(x i+1 ) -θ(x i ) + ai • T ) = T × p i=1 ai . the result.
Theorem 1: Any feasible solution ϕ can be decomposed into a set of circuits S such that, for any arc a ∈ G, ϕ(a) = c∈S ϕ c (a). Proof. The proof is by recurrence on n(ϕ) = a∈A ϕ(a). The theorem is trivially true if n(ϕ) = 0.

Let suppose now that n(ϕ) > 0, thus there exists at least one arc a = (x, y) with ϕ(a) > 0. Set µ 0 = (x, y) and let consider the sequence of paths µ i built as follows:

1) Stop the sequence as soon as µ i contains a circuit c; 2) Otherwise, let ã = (x, ỹ) the last arc of µ i . Since ϕ(ã) > 0, following the conservation low of ϕ over nodes of G, there exists an arc a starting at ỹ with ϕ(a) > 0. We then set µ i+1 = µ i • a.

As G has a finite number of nodes, the algorithm stops and a non empty circuit is then returned. Let now define the flow φ as follows

φ(a) = ϕ(a) -1 if a ∈ c ϕ(a)
otherwise φ is feasible with n( φ) < n(ϕ), thus the theorem.

Note that the number of flow variables is a polynomial function on the size of the problem. It is not the case for the vehicle tours, which number can be of exponential size. The consequence is that the determination of a flow is in N P, which is not the case for the determination of vehicle tours.

D. Modelling of the optimization problem

Three objectives are to be considered for our optimization problem: the main objective is to maximize the number of demands. The two other ones are to minimize both the number of relocations and the total number of vehicles.

As we shall see later, the number of demands and relocations are easily linearly expressed using the flows. Next theorem shows that it is also the case for the total number of vehicles:

Theorem 2: The minimum total number of cars required for a feasible flow ϕ equals a∈A ϕ(a) • a . Proof. Let S be a set of circuits obtained from the decomposition of ϕ following Theorem 1 and let V be the minimum number of cars associated with ϕ. By Lemmas 1 and 2, the total number of car of any circuit c ∈ S is Lastly, the total number of vehicles required at time T to reach ϕ is exactly a∈A ϕ(a) • a , the theorem.

The modelling of our optimization problem follows. R and C are fixed bounds for respectively the total number of relocation operations and vehicles. Equation (1) is the maximization of the total demand. Equation (2) expresses the bound on the total number of relocation. Equation (3) expresses these on the total number of vehicles. Equations (4), ( 5) and (6) are lastly flow constraints.

max a∈A2 ϕ(a) (1)
s.t.                            a∈A3 ϕ(a) ≤ R (2) a∈A ϕ(a) • a ≤ C (3) ϕ(a) ≤ u(a) ∀a ∈ A (4) y∈Γ -(x) ϕ((y, x)) = y∈Γ + (x) ϕ((x, y)) ∀x ∈ X (5) ϕ(a) ∈ N ∀a ∈ A (6)
The total number of equations is around

2|A| + N × T = Θ(N 2 • T ).

III. EXPERIMENTATIONS / RESULTS

In order to test the linear mathematical model, we first need data. Unfortunately, we do not yet have enough relevant and actual data that could be used for study purposes. To cope with that, we propose to generate data using a random generator described in the first subsection. Then, the second subsection gives some discussion about computation times and solutions analysis of small instances, based on generated data. Threedimensional Pareto frontier is especially presented. Finally, third subsection is dedicated to scalability experimentations, particularly the solver's performance when the problem increases in size.

A. Random data generator

A generator has been then implemented to emulate real demand data over time. We decided to focus our study on representing an average weekday demand. As a consequence, all the data depending on time are generated over a 24 hours period, segmented into T time-steps. The total number of time-steps is user-settable and can vary from 24 to 1440 (representing respectively a 60 and 1 min time-step).

Basically, the generator is based on two phases: station and demand generation. The first phase positions N carsharing stations within a given territory. Maximum size for each station is randomly generated using a discrete uniform distribution over an integer interval [Z min , Z max ] given by the user as a parameter. The station positioning is made over two distinct zones: a central area (in general representing the center of the city) included in a larger one (representing the suburbs area), both defined as a square. The generation algorithm takes two additional parameters: the percentage of total area the center must represent and the probability that a station is contained in the center. Once the geographic division is made, every station is then positioned randomly in the area where it belongs.

Then, the second phase generates randomly M demands over time between stations. First of all, the generator has to schedule and position randomly each request over time, which means defining a probability distribution. In order to do so, the generator allows to specify the distribution profile the demand will follow. In other words, it consists on defining, for every couple of hours (thus 12 values), the relative level of demand Dem(t) at this time t. Then, the probability distribution P is obtained normalizing all the values i.e. P(t) = Dem(t)/ t Dem(t). Typically, most profile distributions are very similar to the one presented in Figure 1. Now, the next step is to identify the two stations concerned by this demand at that time: where it starts and where it goes. Usually, in dense urban area, there are two rush hour slots, also called traffic peaks (generally the morning between 7 and 10 o'clock; the evening between 16 and 19 o'clock) for which the demand goes globally in the same direction: from the suburbs to the center during the morning and from the center to the suburbs the evening. The generator allows to define such rush time slots as well as the proportion of total demand during those rushes. Finally, it's also possible to specify an average car speed and a penalty coefficient during rush time for the calculation of travel-time between stations.

B. Experimentations

1) Experimental context: All experimental results addressed bellow were made using an Intel(R) Core(TM) i3-3227U CPU running at 1.9GHz. The linear programs were expressed using the AMPL format [START_REF]AMPL[END_REF] and numerical results were obtained using GLPK v4.52 [START_REF]GLPK -GNU linear programming kit[END_REF].

2) Experimentation of small instances: First of all, we start our experimentations considering N = 10 stations, T = 144 time-steps and M = 500 demands. The upper bound of relocation operations and the number of vehicles belong to the set {0, 10, • • • , 80}. For those fixed values, 30 instances were generated using the random generator previously described for studding the computation time of the optimal solutions and the impact of model's building time.

Table I presents the average computation time µ and the standard deviation σ for both integer linear program ILP and its relaxation to linear program LP. The minimal and maximal values are also specified. The first observation is that computation times remains quite low, in the order of a half-second for LPs and two seconds for ILPs. The major part is taken by the building of the mathematical program which takes 34s (mostly by the conservation law) regardless of which model is built. Secondly, when the solver runs the LP model, almost 8 problems (exactly 7, 66 on average) under 81 admit a non-integer value of the objective function, over the 30 generated instances. This result represents almost 10% of all instances. However, every time we get a non-integer optimal value, the integer one has always the same integer part. Hence it comes that we conjecture the existence of an integer optimal solution for LP and that the difference between the two criteria values comes from rounding errors.

Figure 2 shows a 3-dimensional Pareto frontier for a particular instance. When all the demands are satisfied, the minimal number of vehicles equals 50. In this case, at least 70 relocation operations are necessary. If the number of vehicle increases to 70, 50 relocations are needed. This confirms that the two criteria are in opposition and that there exists a tradeoff area between them.

3) Scalability experimentations: The aim of this second experimentation is to measure the overall computation times following the size of the problem. We observed previously that its most important part comes from the linear program generation.

Table II presents the generation times depending on the number of stations and time-steps. Each measure was obtained from one instance since the time needed for generating the model only depends on the size of the time extended graph. Note that when N = 50, T = 288 and R = C = 80, GLPK founds an optimal solution within 1 20 minute and an integer optimal solution within 1 30 minute. We observe that the generation time grows linearly following the number of arcs of the time extended graph, which is around 700 000 for this instance. Moreover, this is confirmed by the correlation coefficient between generation time and the number of arcs in the graph standing at 95, 9%. This result is coherent since the size of the linear model and the number of variables are both in Θ(|A|).

For the biggest instance (N = 50, T = 288), the relocation arcs represent 98% of the total number of arcs. Therefore, we suggest to increase the scalability of the method by considering for example only short distance relocation arcs or defining them at fixed time-steps.

IV. CONCLUSION

In this paper, we investigated the optimization of a carsharing one-way system which can be integrated in a multimodal transportation system. Specifically, in dense urban context, the demand of those systems may originate from any other modes of transportation such as trains, subways, bikes, etc. The optimization focus on fleet dimensioning when station locations and accurate demand are given. We proposed an integer linear program based on flows and we suggested three optimization criteria: maximizing the overall demand the system can absorbed, and minimizing the number of vehicles and relocation operations needed. We also described how to extract the total number of vehicles from any optimal flow.

Two numerical experiments were made using an opensource solver (GLPK). The fist one investigated computation times over small instances generated by a random generator. It turns out that the latter is negligible compared to the model building time. Furthermore, we observed that for each noninteger optimal value we get solving the LP model, its integer part is always equal to the optimal value of the corresponding ILP model. This let us think that the LP solution must contain rounding errors and might be integer. We intend to clarify this point in further research. A graphical representation of a 3dimensional Pareto frontier served to underline and to confirm that the optimizations criteria are in opposition.

The second experiment aims at specifying the time needed for building the mathematical model when the problem grows. As it is shown in the paper, the number of constraints of the mathematical model as well as the number of decision variables are both in the range of the number of arcs present in the time extended graph. We suggested to reduce this number in order to improve model generation time selecting only those at fixed time-step representing short distance relocation operations.

The choice of GLPK instead of a commercial solver is mainly motivated by our industrial context. The good quality and the scalability of the solutions we obtained using this open-source solver allow us to transfer easily our methodology in an industrial environment.

Lemma 2 :

 2 Let c be a circuit and ϕ c a feasible flow such that ϕ c (a) = 1 if a belongs to c 0 otherwise The minimum number of vehicles to insure ϕ c is (c) T . Proof. The total number of vehicles needed at time t ∈ H for c is clearly |C t (c)|. We first show that |C T (c)| = a∈c a . For that purpose, setting B(c) = {a = (x, y) ∈ c | a = 1}, we prove that B(c) = C T (c). • B(c) ⊆ C T (c): If a = (x, y) ∈ B(c), then as θ(x) ≤ T , ϑ(θ(x), T ) = T -θ(x). Now, since a = 1, (a) = θ(y) -θ(x) + T ≥ ϑ(θ(x), T ) and a ∈ C T (c). • C T (c) ⊆ B(c): Let consider now an arc a = (x, y) ∈ C T (c). Since ϑ(θ(x), T ) = T -θ(x) < (a), we get θ(y) + a • T > T . As θ(y) ≤ T , we necessarily have a = 1 and thus a ∈ B(c). Now, by Lemma 1, |C T (c)| = a∈c a = (c)

  c (a). Now, from Theorem 1, ϕ(a) = c∈S ϕ c (a). Thus, c∈S a∈A a • ϕ c (a) = a∈A a • c∈S ϕ c (a) = a∈A ϕ(a) • a .

Fig. 1 .

 1 Fig. 1. Classical average demand profile in dense areas over the course of a weekday

Fig. 2 .

 2 Fig. 2. 3-dimensional Pareto frontier

TABLE I COMPUTATION

 I

				TIMES IN SECONDS
		µ	σ	min	max
	LP	0, 56	0, 25	< 0, 01	1, 17
	ILP	1, 93	1, 86	< 0, 01	16, 77

TABLE II GENERATION

 II TIME DEPENDING ON THE SIZE OF THE GRAPH (IN SECONDS)

	H T H H H N	10	20	30	40	50
	72	14	59	138	147	310
				≈ 2 30	≈ 3	≈ 5
	144	53	234	343	632	1437
			≈ 4	≈ 5 30	≈ 10 30	≈ 24
	288	211	997	1417	2481	4735
		≈ 2 30	≈ 16 30	≈ 23 30	≈ 41 30	≈ 1h20

ACKNOWLEDGMENT

This research work has been carried out in the framework of the Technological Research Institute SystemX, and therefore granted with public funds within the scope of the French Program "Investissements d'Avenir".