A methodology for planning a new urban car sharing system with fully automated personal vehicles

Abstract

Purpose

The paper concerns a transport system for pedestrian areas, based on a fleet of fully-automated Personal Intelligent Accessible Vehicles. The following services are provided: instant access, open ended reservation and one way trips. All these features provide users with high flexibility, but create a problem of uneven distribution of vehicles among the stations. A fully vehicle based relocation strategy is proposed: when a relocation is required vehicles automatically move among stations. The paper focuses on a methodology that allows to plan the proposed transport system for wide pedestrian areas. The methodology aims to determine the fleet dimension and the relocation strategy parameters which minimize the system cost. The system cost takes into account the level of service and the efficiency. Relocation strategy parameters define when and among which stations relocations should be performed.

Methods

The problem faced is an optimisation problem where the search space is defined by all the possible fleet dimensions and relocation strategy parameters. As this cost function could be a multipeak function and since the search space is discrete and extremely large, a random search algorithm has been adopted. Because of the characteristics of the problem, a parallel optimization technique was required. Given a fleet dimension and relocation strategy parameters, a microsimulator models the activity of each user, as well as the activity of each vehicle over time with the aim of finding the level of service and the system efficiency.

Results, conclusions and application

The methodology has been applied to planning the proposed transport system for the centre of Barreiro, Portugal.

1 Introduction

The paper concerns a new car sharing system for pedestrian areas, like historical city centres. The transport system is based on a fleet of personal intelligent city accessible vehicles (PICAVs). The PICAV vehicle is a one person vehicle which is meant to ensure accessibility for everybody, and some of its features are specifically designed for people whose mobility is restricted for different reasons, particularly (but not only) elderly and disabled people. The PICAVs will travel in traditionally pedestrian areas and allow people to travel around parts of urban areas which are otherwise difficult or impossible for them to access.

The PICAV unit is an electrically powered with limited autonomy range (some 60 km), at low speed (max 25 km/h). The vehicle is about 1.1 m wide and 1.3 m long. Its total mass is about 400 Kg (batteries included) and its cinematic characteristics have been described in Masood et al. [18] and Cepolina [11]. The vehicle can be fully automated or partially automated. The automation is made possible by a number of laser scanners, ultrasonic sensors and also cameras as well as a highly sophisticated control system, which is being developed by INRIA (http://team.inria.fr/imara/publications-2/papers/). When the vehicle is user driven, the level of automation depends on the user’s abilities. The vehicle because of its small dimensions, tiny footprint and on-board intelligence perfectly fits in pedestrian environments [18], disturbing the other pedestrians as little as possible.

Stations with charging devices are spread in the pedestrian area and its border, possibly closed to interchange points. PICAV are available at stations and should be returned to stations.

The PICAV transport system is a second generation car sharing system, and it is meant to provide users a high level of flexibility, similar to private car. One way trips, open ended reservation and instant access are the main characteristics of the proposed transport system. These three features provide users a high level of flexibility. But on the other hand, the risk of unbalancement in the number of available vehicles at stations is high: due to uneven demand, some stations during the day may end up with an excess of vehicles whereas other stations may end up with none. Therefore, relocation is necessary. In literature, several relocation schemes are proposed [7]. In operator-based relocation strategies some operators manually relocate a vehicle or a platoon of vehicles from stations having too many vehicles to stations having too few. Several car sharing systems have adopted this relocation scheme, but sometimes have turned out into a failure due to the heavy staff and management cost. Some user based relocation strategies have been proposed, in Barth and Todd [2], in Cepolina et al. [9] and in Cepolina and Farina [8], which successfully reduce these costs. On the other hand, the level of automation of vehicles is increasing so rapidly that a fully vehicle based relocation strategy will soon be possible.

The proposed transport system implements a fully PICAV based relocation strategy: when relocations are required, a system supervisor has the duty to redirect PICAVs from a station to another. According with the Supervisor hints, PICAV vehicles automatically relocate themselves, thanks to their high level of automation.

This paper focuses on the problem of optimising the fleet dimension and the threshold values that determine relocation occurrences. The study area is described by a street network where stations are localized and by the PICAV transport demand. The optimal fleet dimension and threshold values are the ones that minimise a cost function. The cost function takes into account both the level of service provided to users and the efficiency from the management point of view. Since the cost function could be a multipeak function and since the search space is discrete and extremely large, Simulated Annealing (SA) has been chosen for solving the minimization problem. Given a PICAV fleet and a set of thresholds the related cost function value is assessed by micro-simulation. This is achieved by inputting the street network, the PICAV transport demand and the fleet dimension and the threshold values into the microscopic simulator, which then simulates the behaviour of the transport system during the reference time period and assesses the waiting times of users and the total amount of time spent by vehicles in relocation.

The paper is organised as follows. Section 2 describes the proposed transport system. Section 3 outlines the architecture of the transport system microscopic simulator. Section 4 describes the optimization problem and the random search algorithm. Section 5 presents the case of study of the old part of Barreiro, Portugal. Major conclusions and future work follow.

2 Main characteristics of the proposed transport system and the relocation strategy parameters

A relocation is required when a critical situation occurs. A critical situation occurs when the number of vehicles in a given station at a given time instant goes below the station’s low critical threshold. This situation is referred as ZV i.e. Zero Vehicle [7, 10, 15]. When this condition occurs, the station has a shortage of vehicles and to avoid a possible queue in the next future a request for a vehicle is generated. In Fig. 1 this situation occurs in station j, since the number of vehicles n _j results lower than the station j’s low critical threshold at the time instant the figure refers to. The low critical threshold could be assumed constant in time or a function of time.

Fig. 1

Conceptual layout of the fully PICAV based relocation strategy

When a ZV situation takes place, the vehicle request is addressed only to stations where the number of vehicles is above the low buffer threshold. According with Kek et al. [15], the low buffer threshold is the minimum number of vehicles that a station needs to have in order to be able to send vehicles. In Fig. 1 station k is able to send a vehicle to station j, station z is not in this situation at the time instant the figure refers to.

Among the stations to which the vehicle request could be addressed, the providing one is selected according to two criteria: the closest one (shortest time criterion) and the station having the highest number of vehicles (inventory balancing). The shortest time criterion relates mainly to service levels, while the inventory balancing mainly focuses on cost efficiency. Therefore, an appropriate choice of relocation technique should be made according to the current system situation: in periods of low usage, the most appropriate relocation technique is by inventory balancing whilst in periods of high usage, then the shortest time technique performs best [7, 10].

Another critical situation is the so-called FP i.e. Full Port [7, 10, 15]. When this situation occurs, a station has reached its capacity and therefore vehicles arriving must be redirected to other stations because there is not more free space for parking. In this case, vehicles are redirected and the supporting station is always the nearest one, in order to reduce the disutility created to the user. Experimental data show that a relocation strategy based on preventing ZV occurrences, guarantees that FP occurrences do not take place.

The relocation strategy parameters are the stations’ low critical thresholds and the stations’ low buffer thresholds.

3 The microscopic simulator

The proposed transport system has been modelled according to an object-oriented logic, because it is the most suitable to our targets. The developed simulation model indeed aims to be: simple, modular, easy to understand, easy to modify (in case we wanted to take into account further aspects in our modelization). The language chosen for writing the code is Python 2.5.

The simulator receives as input: the simulation time period, the road network, the transport demand, the PICAV fleet characteristics and the relocation strategy parameters. The simulator allows to track the second-by-second activity of each user, as well as the second-by-second activity of each vehicle. In Fig. 3 a snapshot of the simulation in a given time instant is provided: the green circle is proportional to the actual number of vehicles available at each station; the red circle to the length of the queue at each station and the white circle is proportional to the maximum waiting time at each station. In violet the actual positions of the PICAV vehicles are shown. The simulator graphical user interface has been provided by UCL [20]. The simulator gives as offline output the transport system performances, in terms of level of service (LOS) provided to users and in terms of efficiency from the management point of view. A general scheme of the micro simulator is reported in Figs. 2 and 3.

Fig. 2

General scheme of the micro simulator. The simulator is object-oriented and it involves three classes: user, PICAV and station. This is just a simplified scheme and methods are not reported entirely

Fig. 3

A snapshot of the simulation at a given time instant [20]

3.1 Input data

The micro simulator input data are: the simulation time period, the road network, the transport demand, the PICAV fleet characteristics, the relocation strategy parameters. All these inputs are deterministic. The only stochastic input is the transport demand, as it concerns the user arrival time instant.

3.1.1 The simulation time period

The simulation time period starts when the car sharing system opens to users and ends when the last user returns the PICAV. In the following we refer to a daily simulation time period. The simulation time period could be characterised by peak and off peak phases: for each phase an average pedestrian density k and a PICAV transport demand should be specified.

3.1.2 The road network

The road network includes stations, provided with charging equipments, and the roads in which PICAV vehicles are allowed to travel. The road network is defined by OpenStreetMap.

Stations have been represented through nodes. Between each pair of stations, we take into account only one path, which could be the shortest one or the more attractive one since characterised by a high concentration of shops, museums and other activities. The characteristics of the path we are interested in are: its overall length and the average upslope. These data are necessary for the discharging law of the PICAV vehicle’s battery. The overall length of a path is also required in order to determining the trip duration.

For each couple of stations, the path length and the average upslope are assessed through a routing algorithm written in JavaScript, which interacts with OpenStreetMap. These data are given as data input in the simulator in the form of two matrixes. The matrixes are squared and the number of rows (or columns) equals the number of stations in the network. In the first matrix, the cell ij represents the path length between the origin node i and the destination node j. In the second matrix, the cell ij represents the average upslope of the path between the origin node i and the destination node j.

The travel time between each couple of stations is calculated in the simulator, from the distance between each pair of stations and the PICAV speed. The vehicle speed is assessed from the pedestrian density according to the following model:

PICAV user driven:

$$v=-1.45k+1.58$$

(1)

PICAV automatically driven, relocation trip:

$$v=-1.45k+1.38$$

(2)

where k is the pedestrian density expressed in pedestrians per square meters; v is the PICAV vehicle speed expressed in m/s. The model for assessing the vehicle speed from the pedestrian density is described in detail in Cepolina et al. [6]. The model has been implemented in the micro simulator.

3.1.3 The transport demand

The transport demand refers to a given phase of the simulation time periods and it is given to the simulator in the form of OD matrix. Each row refers to a station of origin, and each column to a station of destination. Each cell gives the hourly number of trips from the station the row refers to, to the station the column refers to.

We consider two trip typologies. A trip by PICAV could be a direct trip on board a PICAV, or a sequence of shorter trips (multitask trip) where one accomplishes a number of short tasks that require short term parking along the street, before finally returning the unit. In both cases what is of interest for the proposed study is the overall duration of the trip. Given an origin, a destination, the path between them, and an average pedestrian density, the trip duration changes according to the trip typology.

Therefore each OD matrix refers to a given phase of the simulation time period and to a trip typology. In the simulation PICAV users are generated with the following characteristics: the origin of their trip by PICAV, the destination, the time at which they appear in the origin and the trip typology. These data are assessed according to the OD matrixes. The time at which a user appears in their origin is randomly generated: if X users have to be generated between 8 and 9 a.m. in a given origin, X casual numbers are extracted within the given time interval and these casual numbers are the exact arrival instants of the X users in the origin.

3.1.4 The PICAV fleet characteristics

The PICAV fleet characteristics are: the fleet dimension, the number of PICAV vehicles at each station at the beginning of the simulation time period, the battery capacity, the charging and discharging laws.

A lithium-ion battery has been selected by MAZEL [19]. The battery is composed of 15 blocks connected in parallel, each composed of 27 cells connected in serial, and provides 202 Ah and 48 V DC.

As it concerns the discharging law, the quantity of discharge is assessed from the overall length and average upslope of the paths. The upslope contributes heavily to the resistances to motion encountered by the PICAV vehicle. If the path is instead descending the recovery in battery charge is so slight that it is neglected.

The battery charging technique is the opportunity charging. The term opportunity charging refers to the charging of the batteries wherever and whenever power is available. The minimum charge level is the quantity of charge necessary to the vehicle to perform the longest trip or relocation journey. Every time a PICAV is returned in a station, a check on its charging level is performed. If the vehicle has a level of charge which is more than minimum charge level, it is available to users and to relocations, otherwise it starts the charging process. The minimum charge level has been determined through simulation, considering the most consuming path, both for users and relocation trips.

3.1.5 The relocation strategy parameters

The relocation strategy parameters are described by two vectors. Their dimension equals the number of stations in the area, the value of each vector component is the station’s low critical threshold for the first vector and the station’s low buffer threshold for the second vector.

A high value of low critical threshold gives rise to a high number of required relocations and to low waiting times, if the fleet is consistent and therefore there are vehicles available for relocation.

The low buffer threshold is greater than the low critical threshold. If the low buffer threshold is much greater than the low critical threshold, the number of satisfied requests for relocations is low because often no stations can provide the vehicles required: this results in an increase of the users waiting times.

If the low buffer threshold is slightly greater than the low critical threshold, the number of satisfied requests for relocations is high; on the other hand, it may occur that at a given time instant a station provides a vehicle and at a following time instant the same station is in shortage of vehicles. This results in an increase in the number of required relocations. As a result, it is necessary to optimize the low critical and low buffer thresholds values for each scenario under study.

3.2 Output data

3.2.1 Level of service (LOS)

LOS measurement are assessed based on the statistical distribution of users waiting times. Castangia and Guala [4] proposed a new LOS measurement scale (shown in table 1) using as reference the 50th, 90th and 95th percentiles of waiting time. The LOS measurement scale ranges from LOS from A (perfect service) to F (completely poor service). All the constraints on the three percentiles of users waiting times should be met to achieve a given LOS. LOS measurements could be assessed for each station or for the overall area, referring specifically to the waiting time of users arriving in each station or to the waiting time of all the population of PICAV users.

Table 1 The LOS assessed according to the percentiles of users waiting time expressed in seconds [4]

3.2.2 Efficiency

An explicit expression to assess the transport system efficiency does not exist. However, according to Barth and Todd [1] and Kek et al. [15], we assess the efficiency according to the following variables:

• fleet dimension;

• number of required relocation trips;

• percentage of vehicles available, with reference to the total fleet dimension, at each simulated time instant.

The values of the first two variables are assessed offline at the end of the simulation; the value of the last variable is assessed online, i.e. during the simulation run.

3.3 Stochastic effects

As the input data are stochastic regarding the users arrival times, the output data, in terms of users waiting times and relocation time, are stochastic as well. According with the criteria given in Law and Kelton [17], 30 runs of the microscopic simulator resulted sufficient to reduce these stochastic effects.

3.4 The validation of the microscopic simulator

Microscopic model calibration and validation require field data collection, as it concerns output data and their relative input data. In the proposed micro simulator, the critical input data is the OD matrix since a transport system like the one proposed does not exist yet and therefore we have no way to assess the demand for it. As it concerns the output data, the critical ones are related to the PICAV positions, the user waiting times and the queue lengths at the stations. Again the collection of these data in the field is not possible since the PICAV prototype is under construction and the transport system is not yet operative.

However, some sub-models implemented in the microscopic simulator have been calibrated and validated. These sub-models refer to: the battery discharging and recharging laws and the relationship between PICAV speed and pedestrian density. The sub-model related to the battery has been calibrated and validated by MAZEL who is in charge to provide the PICAV power system [19]. As it concerns the last sub-model, it has been calibrated and validated with the results of experiments aiming at understanding the relationship between the speed of an electric scooter and the density of pedestrians [6]. The electric scooter was a standard electric scooter designed for disabled people, and not a PICAV unit as this is yet to be constructed. However the authors believe that the speed density relationship related to the electric scooter can be transferred to the PICAV since the cinematic behaviours of the two are similar.

4 The optimization problem for the picav transport system

The problem of assessing the fleet dimension and the relocation strategy parameters is formulated as a constrained minimisation problem. The function to be minimised is a linear combination of the transport system cost C ^s _d and of the user’s costs C ^u _d : both costs refer to the simulation time period and their dimension unit is €/day.

In the formulated problem, the independent variable is a vector s whose dimension is equal to 3:

• the first component is equal to the fleet dimension n _V ;

• the second component is equal to the low critical threshold. Low critical threshold values are taken constant for all the simulation time period and also are assumed to be the same in all stations. The value reported in the second component of the vector s is the low critical threshold in one of the stations;

• the third component is equal to the low buffer threshold. Low buffer threshold values are taken constant for all the simulation time period and also are assumed to be the same in all stations. The value reported in the third component of the vector s is the low buffer threshold in one of the stations.

High critical thresholds have not been optimised since they correspond to the station capacities. The high buffer thresholds are not optimized, since a relocation strategy based on preventing ZVT occurrences, guarantees that FPT occurrences do not take place, as stated before.

The search space has therefore three dimensions, each one is determined by a component of the vector s. The search space is limited as follows:

• the low buffer threshold of each station must be greater or at least equal than the low critical threshold,

• the fleet dimension must be greater than the sum of the low buffer thresholds,

• the fleet dimension must allow a vehicle-to-trip ratio between 0.03 and 0.06 [1].

The constraints of the optimization problem refer to the 50th, 90th and 95th percentiles of users waiting time from which LOS depends.

4.1 The transport system cost

The daily cost of the system, C ^s _d , is equal to:

$$C_d^s=C_d^f+C_d^{\mathit{run}}+C_d^{\mathit{setup}}+C_d^{\mathit{man}}-C_d^{\mathit{tickets}}+C_d^r$$

(3)

For C ^f _d , i.e. the daily cost of the fleet, we assume that the purchase price for each vehicle is 9,000 €. This cost includes the vehicle with a lifetime of 8 years and two lithium-ion battery packs with a lifetime of 4 years each. The cost of adding vehicles to the system increases linearly with each vehicle. The linear increase in vehicle costs is a rather simplistic assumption and does not take into account any economy of scale, but it is a quite frequent assumption [1]. The daily cost of amortization of the fleet has been calculated with a discount rate r, according to the following formula:

$$C_d^f=n_v{c}_v\left[\frac{r{\left(1+r\right)}^{\mathit{lt}}}{{\left(1+r\right)}^{\mathit{lt}}-1}\right]\frac{1}{365}$$

(4)

Where:

• n _v is the number of PICAV vehicles within the fleet: it is the first component of vector s;

• c _v is each PICAV vehicle purchase cost (it includes two battery packs);

• r is the discount rate, equal to 8 %.

• lt (number of years) is the PICAV vehicle lifetime.

C ^run _d , i.e. the daily cost of running the system, includes the maintenance costs and the electricity cost to run the PICAV fleet. Regarding the maintenance costs, and excluding the cost of batteries that we already included in the vehicle purchase cost, electric cars incur in very low costs and we neglect them. As it concerns the electricity needed to charge all the PICAV vehicles during an operative day, it is proportional to the average daily kilometres travelled by the PICAV users and to the length of the relocation trips. Since we assume a constant PICAV transport demand, the average daily kilometres travelled by the PICAV users do not depend on the vector s. Conversely, the length of the relocation trips depends on s. The electricity cost related to relocation trips is taken into account in C ^r _d and it is proportional, through the vehicles velocity, to the total amount of time spent in relocation.

C ^setup _d , i.e. the daily cost of system setup, includes infrastructure cost which is related to the number of parking spaces in each station. This cost depends on the fleet dimension. But due to the small dimension of PICAV vehicles (1.1 m long and 0.9 m wide) [18] we consider it as a flat cost. C ^man _d , i.e. the daily cost of system management, is again a flat cost.

C ^tickets _d is the daily revenue derived from the tickets. The ticket price is the same for all the PICAV users. If we multiply the number of PICAV users by the ticket price, this provides approximate ticket revenues. Since the number of users in a simulation time period has been assumed constant, C ^tickets _d is again a flat cost.

C ^r _d , i.e. the daily cost of relocation, is calculated from the total amount of time that PICAV vehicles have spent relocating during the simulation time period:

$$C_d^r=c_r{\displaystyle \sum _{j=1}^{n_v}{t}_{\mathit{rj}}}$$

(5)

Where:

• c _r is the cost of each minute of relocation, assumed equal to 0.01 €/min; it takes into account also the energy consumed in each minute of relocation.

• t _rj is the time spent by j^th vehicle in relocation during the simulation time period and it depends on the vector s: t _rj (s).

The total relocation time ${\displaystyle \sum _{j=1}^{n_v}{t}_{\mathit{rj}}}$ takes into account the time spent by all PICAV vehicles in relocation and it is a function of s: ${\displaystyle \sum _{j=1}^{n_v}{t}_{\mathit{rj}}}\left(\mathbf{s}\right)$.

4.2 The users cost

The cost supported by users has been measured in terms of the total users waiting time and the cost the users have to pay for the service. The total users waiting time is the overall length of time, in minutes, users have to wait for a vehicle when it is not immediately available. C ^u _d has the following expression:

$$C_d^u=c_w{\displaystyle \sum _{i=1}^m{t}_{\mathit{wi}}}+C_d^{\mathit{tickets}}$$

(6)

Where:

• C ^tickets _d is the overall cost paid by the PICAV users and it is given multiplying the number of PICAV users by ticket price. Since the PICAV transport demand has been assumed constant, C ^tickets _d is again a flat cost. Moreover this term equals the revenue derived from the tickets that appears in the transport system cost, therefore both terms appear in the cost function with opposite signs, and therefore delete each other.

• c _w is the cost of a unit of waiting time, namely 0.10 €/min [12];

• m is the total number of users who have been in a queue during the reference time period;

• t _wi is the waiting time in minutes of the i^th user. It depends on s: t _wi (s);

The total users waiting time in the simulation time period ${\displaystyle \sum _{i=1}^m{t}_{\mathit{wi}}}$ is again a function of s: ${\displaystyle \sum _{i=1}^m{t}_{\mathit{wi}}}\left(\mathbf{s}\right)$.

4.3 The constrained optimization problem

We minimize a function f which is given by the sum of C ^s _d and C ^u _d without the flat cost terms since they do not affect the minimization problem. The f function has the following expression:

$$\begin{array}{cc}\hfill f\left(\mathbf{s}\right)=n_v{c}_v\left[\frac{r{\left(1+r\right)}^{\mathit{lt}}}{{\left(1+r\right)}^{\mathit{lt}}-1}\right]\frac{1}{365}+c_w{\displaystyle \sum _{i=1}^m{t}_{\mathit{wi}}\left(\mathbf{s}\right)}+c_r{\displaystyle \sum _{j=1}^n{t}_{\mathit{ri}}\left(\mathbf{s}\right)}\hfill & \hfill \forall \mathbf{s}\in N^3\hfill \end{array}$$

(7)

N³ is the search space and N is intended as the set of natural numbers.

The constrained minimization problem is therefore

$$\begin{array}{cc}\hfill \mathrm{Minimize}\hfill & \hfill f\left(\mathbf{s}\right)\hfill \end{array}$$

(8)

subject to:

$$g\left(\mathbf{s}\right)=\left\{\begin{array}{c}\hfill t_w^{50\%}-4<0\hfill \\ \hfill t_w^{90\%}-8<0\hfill \\ \hfill t_w^{95\%}-10<0\hfill \end{array}\right.$$

(9)

Where:

• t ^50 %_w , t ^90 %_w , t ^95 %_w are the 50th, 90th and 95th percentiles of users waiting time, expressed in minutes.

We transform the constrained minimization problem into a single unconstrained problem using penalty functions. The constraints are placed into a new objective function h(s) via a penalty parameter $\widehat{\mu }>0$ in a way which penalises any violation of the constraints:

$$\begin{array}{cc}\hfill h\left(\mathbf{s}\right)=f\left(\mathbf{s}\right)+\widehat{\mu }\cdot {\Sigma }_i{\left[\max \left\{0,g_i\left(\mathbf{s}\right)\right\}\right]}^2\hfill & \hfill \forall s\in N^3\hfill \end{array}$$

(10)

Where: g _i is the i^th constraint.

If g _i (s) ≤ 0, then the [max{0, g _i (s)}]² = 0, and no penalty is incurred. On the other hand, if g _i (s) > 0, then [max{0, g _i (s)}]² > 0 and the penalty term $\widehat{\mu }{\left[\max \left\{0,g_i\left(\mathbf{s}\right)\right\}\right]}^2$ is applied.

The weight of the penalty function $\widehat{\mu }$ is unknown and we determine it with an iterative process, according with Bazaraa et al. [3], starting from a small value and increasing it progressively. In fact, as the value of μ increases, the optimal solution of h(s) gets arbitrarily close to the optimal solution of the original constrained problem: minimize f(s), subject to: g(s) ≤ 0. The initial value of the penalty parameter is taken as μ₁ = 0.1, and is updated as follows: μ_k + 1 = βμ _k , where the scalar β is taken as 10.0 [3] and k and k +1 are two successive iterations.

For a given μ _k , we have a function h _k (s) to minimise. If the vector s_* that minimise h _k (s) is the same vector for which g(s_*) ≅ 0 we stop the iterative process and $\widehat{\mu }={\mu }_k$.

For a given value for the penalty parameter μ _k and a given point in the search space s _o , we calculate the value of the h _k (s _o ) function simulating the transport system with the microscopic simulator that has been described in the section 3.

Since there is no analytical expression for h _k (s), we cannot exclude the need to deal with a multi-peak function and the risk of reaching a local minimum, without being able to find the global minimum, is high [5]. To combat this issue and the fact that the search space is extremely large, Simulated Annealing (SA) has been chosen to solve the minimization problem. The procedure for solving the minimisation problem through the Simulated Annealing is described in the following section.

4.4 The Simulated Annealing

The Simulated Annealing (SA) scheme is a stochastic method currently very popular for difficult optimization problems. The term Simulated Annealing is motivated by an analogy to annealing in solids searching for minimal energy states. This procedure starts with the metal at a liquid state and at a very high temperature. In this state the atoms are quite free in their movements. The temperature of the metal is then slowly lowered. If the metal is cooled slowly enough, the atoms are able to reach the most stable orientation. This slow cooling process is known as annealing and so the method is known as Simulated Annealing.

The method is an iterative process that searches from a single point moving in its neighbourhood and allows sometimes to accept worse solutions. This is meant to avoid to get stuck in a local minimum in the optimization procedure. Worse solutions are accepted according to a probability, which depends on a parameter, i.e. the temperature, which decreases with the number of steps.

The algorithm evolves through an iterative cycle, in which the search space is explored. This search depends on a control parameter called temperature T which decreases as the number of the iteration of the cycle increases. In each iteration, a new point s _n is reached from s _o , according to the transition rule. At the new point, the value of the cost function h is checked.

Since the cost function does not have an explicit formula, at each step of the Simulated Annealing algorithm, the microscopic simulator is recalled to calculate the users waiting times and the relocation times from which the cost function value depends. The diagram in Fig. 4 represents the cost function value assessment in the case s _o has only 2 components.

Fig. 4

For each point of the search space, the cost function value is assessed by micro simulation

The updating happens according to:

• a) if h(s _n ) ≤ h(s _o ) → s _n substitutes s _o , i.e. s _o : = s _n

• b) if h(s _n ) > h(s _o ) → s _n will become the current solution s _o with a probability given by:

  $$p=exp\left(-\frac{h\left({\mathbf{s}}_{\mathbf{n}}\right)-h\left({\mathbf{s}}_{\mathbf{o}}\right)}{T}\right)$$

  (11)

This is the core of Simulated Annealing and is known as the Metropolis algorithm. T is the value of the temperature for the current cycle [16]. Given that r ∈[0,1] is a pseudo random number, updating happens according to the following:

• if r ≤ p → the new solution s _n substitutes s _o ,

• if r > p → the new solution s _n is rejected and therefore s _o will not be updated.

Therefore the algorithm needs the definition of the cooling schedule, the local search and the starting and stopping conditions.

4.4.1 The cooling schedule

The cooling schedule is defined by the initial temperature, the law of its decrease and the final temperature. The starting temperature has been determined according to Laarhoven and Aarts [16].

An initial acceptance ratio p ₀ of the worse solution, e.g. 0.5, is fixed at the first step of the algorithm. From this point, the initial temperature T ₀ is determined from the acceptance ratio p ₀ in this way, according to Laarhoven and Aarts [16]:

$$0.5=p_0=exp\left(-\frac{h\left({\mathbf{s}}_{\mathbf{n}}\right)-h\left({\mathbf{s}}_{\mathbf{o}}\right)}{T_0}\right)$$

(12)

The choice of the initial acceptance ratio has the purpose of performing a quite good exploration of the search space without slowing down too much the algorithm.

As in Cepolina [5], the geometric temperature reduction function has been used: T_k + 1 = α ⋅ T _k where T _k and T_k+1 are the temperatures in two consecutive iterations of the algorithm. Typically, 0.7 ≤ α ≤ 0.95. In order to have a good exploration of the search space but not a too slow algorithm, α has been assumed equal to 0.9.

The final temperature scheme of the cooling schedule is replaced by a stopping condition. The algorithm is stopped when 100 iterations without accepting any more new solutions is reached, according to the stopping criteria given in Laarhoven and Aarts [16].

4.4.2 The transition rule

The transition rule regards the exploration of the search space: from a given vector s _o , a new vector s _n is selected in the neighbourhood of s _o .

The transition rule is probabilistic: it passes from s _o to s _n changing only one component of the vector s _o . The algorithm randomly determines the component of the vector to modify. In our case, each component has the same probability to be selected. The algorithm also determines whether to increase or decrease the chosen component: it is increased with a probability of 50 % and it is decreased with the same probability. More specifically, the first component of s _o , i.e. the fleet dimension, if selected, is increased or decreased by m, where m is the number of stations in the intervention area. The second and the third component of s _o , the low critical and low buffer thresholds, if selected, are increased or decreased by 1. Moreover, the algorithm avoids the situation where, in a given iteration, the vector component to change is the same as the one that has been changed in the previous iteration. In this way, it is guaranteed that the new vector s _n is taken in the neighbourhood of the previous vector s _o .

4.5 Parallel optimization

The optimization procedure previously described has shown a serious problem: the objective function is heavily dependent on the fleet dimension (first component of s) and slightly dependent on threshold values (second and third components of s). This causes a sensible slowdown of the optimization algorithm.

Therefore, it has been decided to split the search space into two components and to work in parallel on two processors: on a processor the objective function is kept dependent only on the fleet dimension and all threshold values are fixed, while on the other processor the objective function depends only on the low critical and buffer thresholds while the fleet dimension is kept constant. This technique is in accordance to the search space decomposition methodology [13].

Search space decomposition refers to the case where the problem domain, or the associated search space, is decomposed and a particular solution methodology is used to address the problem on each of the resulting components of the search space. The chosen solution methodology is the Simulated Annealing for both components. The two SA processes are not fully independent and data exchange occurs at the end of each run of the parallel optimization algorithm. This is a simple parallel optimization technique. In fact parallel/distributed computing means that several processes work simultaneously on several processors solving a given problem instance [13]. The application of the proposed methodology allows the reduction of calculation time by 70 %. Indeed, without the parallelism, the parameter α for the decrease of the temperature must be taken equal to 0.98 instead of 0.90 if we want that successive runs of the algorithm converge to the same solution. Instead, with about 3–4 iterations of the parallel optimization, the algorithm converges.

Therefore the transition rule described in paragraph 4.4.2 changes on the 1st processor only the first component of s while on the 2nd processor only the last two components of s. According with Fig. 5, the parallel optimization steps are the following:

Fig. 5

The parallel methodology implemented in the optimization procedure. In the figure, the fleet dimension is referred as Fl, and the low critical and low buffer threshold values are referred as Th

• Initialization:

  • The 1st processor receives in input the thresholds Th ₀ and it optimizes the fleet Fl ₁ . The optimal fleet, given the threshold values, is calculated through Simulated Annealing.

  • When the 1st processor ends its optimization algorithm, the 2nd processor receives in input the optimal fleet dimension Fl ₁ from the 1st processor, and optimizes the threshold values Th ₁ through Simulated Annealing. The point s^* ₁ is given by the optimal fleet dimension Fl ₁ and the optimal threshold values Th ₁

• k^th iteration

  • the two processors work together; they both receive in input the point s^* _k-1 . The 1st processor optimizes the fleet Fl _k through keeping the thresholds constant (Th _k-1 ) and the 2nd processor optimizes the thresholds Th _k keeping the fleet dimension constant (Fl _k-1 ).

  • When both the processors end their optimizations, the values of the objective function in the two points h(Fl _{k ,} Th _k-1 ) and h(Fl _{k-1 ,} Th _k ) are compared. The departure point for the following iteration of the parallel optimization is the point s^* _k which gave the minimum value of the objective function

• Stopping conditions

  • if the values of the objective function in s^* _k-1 and s^* _k differ less than 5 %,

  • or if the two processors do not find any better solution than the initial one

4.6 Validation of the optimization procedure

The calibration and validation of the overall optimization methodology consists in the calibration and validation of the cost function parameters and of the optimization methods (SA and parallel optimization).

Regarding the cost function, the r and c _w values have been taken from the literature. We choose a discount rate of 8 % as it is an average rate of return for investments and therefore reasonably represents the opportunity cost of the purchase.

Given these values, the value of c _r has been assessed in such a way that the transport system cost and the users cost assume comparable values when calculated in the optimal s^* point. Under this hypothesis, the value of c _r has been assumed equal to 0.01 €/min.

Regarding the optimization algorithms, since the optimization methods converge to the same solution, with different starting points and in a reasonable amount of time, the values assumed for the parameters could be considered correct.

5 Barreiro case study

The proposed transport system has been planned and simulated for the old part of Barreiro, a suburb of Lisbon, Portugal. This village is on the edge of a peninsula, on the left side of the Tago River. The old part is about 1 km². The area is almost flat. An air view of the intervention area is given in Fig. 6.

Fig. 6

Aerial view of Barreiro [14]

In the intervention area, i.e. the old part of Barreiro, roads form a grid and vehicles are allowed to circulate in the area. The road pavement is quite irregular and therefore the speed is limited. Several parking spaces exist in the most peripheral parts of the intervention area and the parking here is for free. However, in the most inner part of the intervention area, due to narrow roads, the parking is restricted. Some images of the intervention area are provided in Figs. 7 and 8.

Fig. 7

A street on the border of the intervention area

Fig. 8

A street in the inner part of the intervention area

Several bus lines cross the old part of Barreiro, and several stops are in the area.

The railway line passes on the border of the old part of Barreiro and two are the railway stops of interest: Barreiro centre and the one at the fluvial terminal, as shown in Fig. 9.

Fig. 9

The PICAV intervention area

We identified 5 internal centroids and 3 external ones. Their positions are shown in Fig. 9. The old part of Barreiro is in part a residential area (zones 1, 3 and 4 in Fig. 9), where most of the population (about 18,000 residents) is aged and in part a commercial area (partially zone 2 and mainly zone 5 in Fig. 9). Zone 5 is characterised by a big new commercial centre (Forum Barreiro). The three external centroids correspond to the fluvial terminal (centroid 6), the Barreiro A railway station (centroid 7) and a big car park space (centroid 8).

According with the data provided by Transportes Collectivos do Barreiro (TCB) [14], the trips that interest the area are: trips that have the intervention area as origin or destination and trips internal to the area, that have origin and destination inside the area. There are not trips that cross the area and have origin and destination outside it.

Trips to/from Barreiro happen mainly:

• by boat, to/from the fluvial terminal. The trips by boat connect Barreiro with Lisbon.

• by train: to/from Barreiro centre station. Trips by train connect Barreiro to the villages in the mainland, on the other side of Lisbon.

• by car: trips to/from Barreiro centre take place along three roads.

Internal trips take place mainly by bus and by foot.

5.1 The simulation of Barreiro case study

The O/D matrixes for this case study have been assessed from data collected on the field by TCB [14]. These data regard the total transport demand in the study area. TCB has also made some hypotheses on the modal shift of a quota of the total transport demand to PICAV system. The hypotheses made are the following:

• 10 % of commuters arriving in Barreiro centre by car will be PICAV users;

• 20 % of commuters arriving in Barreiro centre by boat and by train will be PICAV users;

• 20 % of trips internal to the intervention area take place by PICAV. In this case users are mainly elderly and disabled people;

• 0 % of commuters arriving in Barreiro centre by bus will be PICAV users;

The overall number of daily trips by PICAV results 1296.

The simulation time period is a working day and two phases have been identified: a morning period from 7 a.m. to 1 p.m. and an afternoon period from 1 p.m. to 7 p.m. The demand has been assumed constant in each phase. The overall number of trips is the same in the two phases but the OD matrix in the second phase is the transpose matrix of the first phase. In the simulation, pedestrian density has been assumed the same in the two phases and equal to the greatest one registered in Barreiro footpaths (0.2 pedestrians per square meter).

We identified 8 PICAV stations in the area: each one refers to a centroid. The exact localization of stations takes into account the available space and stations are placed as close as possible to interchange points. Their localization is shown in Fig. 7. The station capacity has been assumed equal to 15 PICAV for each station. The high critical threshold is equal to 14 vehicles for each station and the high buffer threshold is equal to 12 vehicles for each station. These values guarantee to keep high the number of supporting stations, and on the other hand to avoid that a station in a given instant accepts a vehicle and in the following instant reaches the FPT situation.

The charging technique is the opportunity charging technique and the minimum level of battery charge is equal to 6 %.

At the beginning of the simulation period, the fleet has been assumed equally distributed among stations.

The low critical and low buffer thresholds, as well as the fleet dimension, are the result of the optimization procedure. The optimum fleet dimension resulted equal to 72 vehicles, i.e. 9 vehicles per station. And this is consistent to the results by Barth and Todd [1] regarding the Coachella Valley transport system, for which the vehicle-to-trip ratio is comprised between 0.03 and 0.06 vehicles per trip. Indeed for the fleet dimension of 72 vehicles, i.e. 9 vehicles per station, the vehicle-to-trip ratio is equal to 0.055. The optimized low critical threshold is equal to 2, the optimized low buffer threshold is equal to 4. The optimised fleet dimension and relocation strategy parameters provide an average waiting time equal to 0.20 min. The 95th percentile of user’s waiting times results equal to 1.84 min, the 90th percentile equal to 0.04 min and the 50th percentile equal to 0 min. This provides a level of service C. The given values of the percentiles of waiting times have been obtained by averaging the related values over the 30 runs of the micro simulation.

As it concerns the efficiency of the optimised transport system, in Fig. 10 the number of PICAVs is each state (occupied by users, available, in charge and relocating) is plotted against time: the time is expressed in minutes starting from 420, that refers to 7 a.m., therefore a time of 480 refers to 8a.m, a time of 960 refers to 4 p.m. and a time of 1,200 refers to 8 p.m. Figure 11 shows the distribution of users waiting times: the total number of users is 1,296 and 1,082 users have not been in queue. The optimised threshold values result efficient since: on one side, relocations are possible in fact waiting times are low (as shown in Fig. 11); on the other hand the number of relocations is kept low as it results in Fig. 10.

Fig. 10

Number of vehicles in each state during the simulation period

Fig. 11

Distribution of users waiting times

The relocation cost is equal to 368.37 €/day, the cost of users waiting times is equal to 83.86 €/day and the objective function results equal to 769.48 €/day. It results evident that the cost of users is much lower that the system cost, this is due to the fact that the cost of waiting times decreases strongly with the fleet dimension; conversely, the relocation cost increases slightly with the fleet dimension. However it should be noted that the function values do not have any meaning since the function does not take into account flat costs, like for instance the ticket prices.

5.2 Sensitivity analysis

With reference to the Barreiro scenario, a sensitivity analysis has been performed. The total transport demand has been increased by 10 % keeping the same optimized fleet dimension and threshold values. The cost function value increases from 769.48 €/day to 806.98 €/day, with a limited percentage increase (4,6 %). However, the users cost increases from 83.86 €/day to 399.78 €/day while the transport system cost decreases from 685.62 €/day to 407.2 €/day. Since the optimised transport system is not able to cope with the demand increase, it results necessary to optimise the transport system parameters for the new value of the transport demand.

6 Conclusions

A new shared vehicle system for pedestrian areas is proposed in this paper; its main characteristics are: one way trips, open ended reservation and instant access. These three features provide users a high level of flexibility. But on the other hand, the risk of unbalancement in the number of available vehicles at stations is high, therefore, relocation is necessary.

A new management strategy is proposed in the paper, based on vehicles that are fully automated and therefore move without an operator among stations in order to relocate themselves. The new management strategy is defined by relocation strategy parameters that define when and among which stations relocations should be performed.

In order to plan such a vehicle sharing system for a given pedestrian area, an optimization procedure is presented in the paper which allows to assess the relocation strategy parameters that minimize the system cost, both in terms of level of service provided to users (that depends on waiting times) and the efficiency from the management point of view (that depends on relocation time and fleet dimension).

Since there is not an explicit expression for the cost function, the distribution of users waiting times and the total amount of time spent by vehicles in relocation, from which system cost depends, are assessed by microscopic simulation. The microscopic simulator follows an object – oriented logic. The simulator follows each users and each vehicle within the simulation period, and gives the actual users waiting times and the relocation time.

As illustrative problem, the proposed transport system has been planned for the old town of Barreiro, a suburb of Lisbon, Portugal. The results of the simulation clearly show the effectiveness of the proposed car sharing system, because, with low staff costs, it allows users a high level of satisfaction. The model has been validated through a comparison of the simulation output data with those available in the body of knowledge.

The automatic relocation of PICAV vehicles still cannot be applied on the field because of legal problems in case of accident. To reduce the impact of automatically driven vehicles, also at legal level, it could be explored relocation by platooning. The operator drives the first vehicle of a platoon and the other vehicles follow the leader through automatic distancing. This relocation technique, however, increases the staff costs, as some operators to perform the relocation are needed. This cost could be reduced by optimization of relocation trips.