Yousra Chabchoub
email: yousra.chabchoub@isep.fr

Rayane El Sibai
email: rayane.el-sibai@isep.fr

Christine Fricker
email: christine.fricker@inria.fr

Bike sharing systems: a new incentive rebalancing method based on spatial outliers detection

Keywords: Outliers detection, spatial data mining, Moran scatterplot, Bike Sharing Systems Data streams processing, Data summarization, Data quality, and Cloud Computing

Since its launch, Velib' (the Bike Sharing System -BSS-in Paris) has emerged in the Parisian landscape and has been a model for similar systems in many cities. A major problem with BSS is the stations' heterogeneity caused by the attractivity of some stations located in particular areas. In this paper, we focus on spatial outliers defined as stations having a behavior significantly different from their neighboring stations. First, we propose an improved version of Moran scatterplot to exploit the similarity between neighbors, and we test it on a real dataset issued from Velib' system to identify outliers. Then, we design a new method that globally improves the resources' availability in bike stations by adapting the users' trips to the resources' availability. Results show that with a partial collaboration of the users or a limitation to the rush hours, the proposed method enhances significantly the resources' availability in Velib' system.

Introduction

In order to reduce congestion and pollution in cities, Bike Sharing Systems (BSS) have been widely adopted around the world. The use in these systems is quite simple: a user takes a bike from a given station with a finite number of docks (called capacity of the station) and returns it at the arrival station, near his destination. BSS provides a 24/7 service, complementary to the other public means of transport, with high flexibility and a relatively low cost for the city. But the system is not so easy to manage.

Our work is motivated by the problem of lack of resources, both bikes, and available slots, in bike sharing systems. Indeed, some stations are often almost empty, without enough available bikes or almost full causing users dissatisfaction. In this context, we define outliers as stations with one of the resources lacking while the other stations in the neighborhood are globally balanced. First, we use an adapted version of Moran scatterplot to explore and characterize the neighborhood of such stations. The results show a local heterogeneity in the system: in a small area, bikes availability is often very variable, depending on the station. This local heterogeneity motivates us to propose a new incentive method which encourages users to improve the placement of bikes among the stations. This mechanism is based on a local small change in users trips. In this ecological regulation, users are redirected to another station in the neighborhood of their source or destination to locally reduce the heterogeneity of stations. Even if this problem is addressed for Velib, the BSS of Paris because we have access to data, we are convinced that our study could apply to other BSSs or electric car sharing systems [START_REF] Ait-Ouahmed | Relocation optimization of electric cars in one-way car-sharing systems: modeling, exact solving and heuristics algorithms[END_REF] and also [START_REF] Illgen | Literature review of the vehicle relocation problem in one-way car sharing networks[END_REF] for a survey on car sharing systems.

The main tool in this paper is the detection of anomalies in a spatial context, where an object is considered as an outlier if the values of its non-spatial attributes are significantly different from those of other objects in its surrounding. Spatial outliers detection is useful in many applications, such as the detection of abnormal highway traffic patterns [START_REF] Shekhar | Detecting graph-based spatial outliers: algorithms and applications (a summary of results)[END_REF], the identification of disease outbreaks [START_REF] Wong | Rule-based anomaly pattern detection for detecting disease outbreaks[END_REF], the detection of tornadoes and hurricanes [START_REF] Lu | Wavelet fuzzy classification for detecting and tracking region outliers in meteorological data[END_REF] and the identification of urban soils pollution [START_REF] Zhang | Use of local moran's i and gis to identify pollution hotspots of pb in urban soils of galway, ireland[END_REF].

Outliers are defined as a set of observations that are inconsistent with the remainder observations. Outliers identification has practical applications in many areas, such as intrusion detection, fraud detection, fault detection and medical informatics [START_REF] Chandola | Anomaly detection: A survey[END_REF]. Outliers detection is also an important task in the data analysis process. It aims to detect abnormal patterns and leads to the identification of unusual phenomena, and new knowledge about the monitored environment. To isolate outliers, it is necessary to first characterize the normal observations, which can be provided by the past values of the same object or by the current values issued from other objects in the neighborhood. In this latter case, the outlier is said spatial. In a spatial context, each data is defined with two categories of attributes: spatial attributes and nonspatial attributes. Spatial attributes include the shape, position, and other topological characteristics of the sensor, and they are used to define the neighborhood of the spatial object. Non-spatial attributes include the ID, manufacturer, age, and sensor measure (called behavioral attribute). A spatial outlier represents a local instability and is only compared to the surrounding dataset [START_REF] Shekhar | Identifying patterns in spatial information: A survey of methods[END_REF]. This is based on the rule: "Everything is related to everything else, but nearby things are more related than distant things" [START_REF] Tobler | A computer movie simulating urban growth in the detroit region[END_REF]. This paper is organized as follows. We present in Section 2 the existing related works. We describe in Section 3 the dataset for Velib' used in this work and we discuss the problem of allocation of resources in bike sharing systems. We describe in Section 4 the so-called Moran scatterplot technique and propose an original adaptation to Velib' context. We also detail in this section the experiments carried out on Velib' to illustrate the heterogeneity of the system. We present in Section 5 our new solutions to balance BSS and to improve the resource balancing between the stations and validate them on Velib'. Section 6 concludes the paper.

Related works

Several algorithms have been developed to detect the outliers in a spatial context. These algorithms can be classified into two categories: graphicalbased algorithms, and quantitative-based algorithms. Graphical-based algorithms use visualization. They present, for each spatial point, the distribution of its neighbors and identify outliers as points in specific regions. This category includes Variogram Cloud, Pocket Plot, Scatterplot, and Moran scatterplot methods. Quantitative-based algorithms perform statistical tests to distinguish the outliers from the rest of the data. These methods include Z-Algorithm, Iterative r Algorithm, Iterative z Algorithm, and Median Algorithm.

Scatterplot represents the data in a two-dimensional space where the X-axis represents the values of the nonspatial attribute (the observable) of each object and the Y-axis represents the mean value of the observable of the neighbors of this object. A regression line is used to identify outliers points [START_REF] Haining | Spatial data analysis in the social and environmental sciences[END_REF]. Variogram Cloud [START_REF] Haslett | Dynamic graphics for exploring spatial data with application to locating global and local anomalies[END_REF] is a scatterplot between the spatial distance (X-axis) and the difference of the observable values (Y-axis) for each pair of points in the data. Outliers points are identified as pairs of points having a small spatial distance and a big difference for the observables measurements.

The Z statistic approach [START_REF] Shekhar | A unified approach to detecting spatial outliers[END_REF] is one of the most known quantitative-based algorithms for spatial outliers detection. For each spatial object x, S x denotes the difference between the attribute value of x and the average attribute value of its spatial neighbors. Spatial outliers are simply identified using a threshold based on µ s and σ s which respectively represent the mean and the standard deviation of the attribute value of S over all the spatial objects.

In [START_REF] Lu | Algorithms for spatial outlier detection[END_REF], authors propose two iterative algorithms (iterative r and iterative z) for the detection of spatial outliers. These algorithms detect the outliers on several iterations. Each iteration detects a single outlier and modifies its value in order to reduce its negative impact on its neighbors in the next iteration.

We apply in this paper the spatial outlier detection approach to a particular case study: the performance evaluation of a balancing mechanism in Bike Sharing Systems (BSS). Nowadays, public authorities are more and more encouraging this ecological mean of transportation by expanding the BSS to the suburbs and building new bike paths. Since its launch in 2007, Velib' (the BSS in Paris) has emerged in the Parisian landscape and has been a model for similar systems in many international cities. Velib' provides a significant proportion of people travels as it daily ensures about 110, 000 trips. It involves about 1800 stations with an average distance of 216 meters.

Literature for BSS is now huge. BSS is widely seen as a healthy sustainable urban mobility mode. Indeed many papers present their positive impact on energy consumption and emissions, safety and economy and traffic conditions. See Zhang and Mi [START_REF] Zhang | Environmental benefits of bike sharing: A big data-based analysis[END_REF] and reference therein. A major problem in BSS is the problem of empty stations and full stations caused by the asymmetric attendance to the stations. According to the annual satisfaction survey of Velib', only 50% of users are satisfied with the availability of bikes and docks in the stations [START_REF]Tout sur Vélib[END_REF], [START_REF] Sibai | Using spatial outliers detection to assess balancing mechanisms in bike sharing systems[END_REF]. Despite the performed redistribution (moving bikes using trucks), users often find themselves in front of stations that are totally full or empty. The main approach used to study redistribution is optimization for the static one-vehicle problem in Chemla et al. [START_REF] Chemla | Bike sharing systems: Solving the static rebalancing problem[END_REF] (see also Cruz et al. [18] and references therein) or multi-vehicle problem in Forma et al. [START_REF] Forma | A 3-step math heuristic for the static repositioning problem in bike-sharing systems[END_REF] and more recently the dynamical problem in Dell'Amico et al. [START_REF] Dellamico | The bike sharing rebalancing problem with stochastic demands[END_REF]. And very few have been done for designing user incentives for rebalancing strategies. Fricker and Gast [START_REF] Fricker | Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity[END_REF] proposed explicit 2-choice algorithms, well-known for balancing queueing networks and Haider et al. [START_REF] Haider | Inventory rebalancing through pricing in public bike sharing systems[END_REF] propose user incentives to decrease redistribution and solve an optimization problem.

In most cities, operators provide open access to realtime status reports on their bike stations. Several studies show the interest of using these data (Froelich et al. [START_REF] Froehlich | Sensing and predicting the pulse of the city through shared bicycling[END_REF] and Borgnat et al. [START_REF] Borgnat | Shared bicycles in a city: A signal processing and data analysis perspective[END_REF], Vogel and Mattfeld [START_REF] Vogel | Understanding bike-sharing systems using data mining: Exploring activity patterns[END_REF]). The main objective in the first papers is to understand and characterize the behavior of the users in order to help in designing and planning policy in urban transportation. See a comparative study of the trip data in [START_REF] Zaltz | The structure of spatial networks and communities in bicycle sharing systems[END_REF] and docking data [START_REF] Obrien | Mining bicycle sharing data for generating insights into sustainable transport systems[END_REF] in different cities and [START_REF] Borgnat | Shared bicycles in a city: A signal processing and data analysis perspective[END_REF] for classification of the flows of trips in the Lyon bike sharing system, using spatiotemporal characteristics to perform clustering. See also Côme et al. [START_REF] Etienne | Model-based count series clustering for bike sharing system usage mining: a case study with the vélib'system of paris[END_REF] and Bouveyron et al. [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF]. But these data are now used for rebalancing purposes. The first paper as far as we know is Gast et al. [START_REF] Gast | Probabilistic forecasts of bike-sharing systems for journey planning[END_REF].

The same trends are present for free-floating BSS (FFBSS), recently deployed in many cities. They have to face a more serious imbalance due to less restriction on parking locations. Moreover, rebalancing is more expensive as bikes may be recovered independently. Furthermore, fleet dimensioning and location of the bikes appear as main issues as they can be fixed by optimization in Zhang et al. [START_REF] Zhang | Electric fence planning for dockless bike-sharing services[END_REF] or deep-learning using data in Pan et al. [START_REF] Pan | A deep reinforcement learning framework for rebalancing dockless bike sharing systems[END_REF]. And this issue is also studied for docked BSS. See Park and Sohn [START_REF] Park | An optimization approach for the placement of bicycle-sharing stations to reduce short car trips: An application to the city of seoul[END_REF].

Dataset description and problem definition

In order to promote innovation and collaboration with scientists, different kinds of data relative to the Velib' system are "Open Data" available for the research community. We performed all the experiments presented in this paper from these datasets.

First, we have static data describing the Velib' stations. They consist of spatial attributes: the geocoordinates of the station (latitude and longitude), and non-spatial attributes: Id of the station and its capacity (total number of docks). Then, the dynamic data are of two kinds: First, the number of bikes present in each station for each timestamp t are provided in real time. This parameter is varying during the day and is closely dependent on the activity of users. Second, Velib' users' data trips are also available (one file for the trips during a month). A trip is characterized by a departure and arrival timestamp, and a departure and arrival station. The analysis of several months of trips showed a very strong periodicity: trips can be divided into two main categories: the working days and the weekends. Two days of the same category are very similar. We focus in this paper on the working days and we choose to analyze trips during one day: trips that took place on Thursday, October, the 31 th , 2013. This duration includes 121, 709 trips, involving 1226 Velib' stations. 1.03% of the trips are related to maintenance (bikes taken for repair) and 1.48% are trips of regulation (bikes moved by trucks).

According to many research studies ( [START_REF] Chabchoub | Classification of the vélib stations using kmeans, dynamic time wraping and dba averaging method[END_REF], [START_REF] Fricker | Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity[END_REF] and [START_REF] Etienne | Model-based count series clustering for bike sharing system usage mining: a case study with the vélib'system of paris[END_REF]), the Velib' system has some weaknesses caused by the strong attractivity of some stations that can be explained by their location near a railway station or a monument or a business area. Such stations are very often almost empty (no available bike) or completely full (no available dock to put a bike). Despite the performed regulation (bikes moved by trucks), the system is still unbalanced. This unbalanced distribution of bikes among the different stations causes users dissatisfaction. The unbalanced stations are referred to as problematic stations. More precisely, we introduce the following definition: a station is said problematic at timestamp t if its occupancy ratio is under 10% or more than 90%.

The occupancy ratio of the station, at a timestamp t, is defined as follows:

occupancy rate t =
N umber of bikes present at t Capacity of the station

Spatial outliers detection with an improved Moran scatterplot

The objective of this section is to estimate the number of isolated problematic stations at a given timestamp t, which motivates the incentive method detailed in the next section. A good understanding of the current use of the Velib' system and the real needs of the users is mandatory to improve the performance of this system and to plan its future expansion and evolution. An isolated problematic station satisfies both following conditions: First, it is almost empty or almost full at timestamp t. Second, its occupancy ratio is significantly different from the average occupancy of the neighboring stations at the same timestamp t. Thus the isolated problematic stations are among the spatial outliers.

In this section, we consider the system at a fixed timestamp t. In order to detect spatial outliers, we opted to use Moran scatterplot [START_REF] Anselin | The Moran scatterplot as an ESDA tool to assess local instability in spatial association[END_REF] that we adapted to the specificities of our context.

Moran scatterplot

Moran scatterplot [START_REF] Anselin | The Moran scatterplot as an ESDA tool to assess local instability in spatial association[END_REF] illustrates the similarity between an observed value and its neighboring observations. It measures the global spatial autocorrelation over a geographical area, the well-known Moran's I. Let us denote by Z = {z i : 1 ≤ i ≤ n} the set of the different values of the considered observable at a fixed given time t, in n different locations. For each location, the neighborhood is defined based on the geographical distance. Moran scatterplot visualizes the relationship between the values z i and their neighborhood average

W i • Z,
Z i × j w ij Z j < 0
W is the contiguity matrix of weights. It indicates the spatial relationship between every couple of objects. W is also called the row-normalized neighborhood matrix. It is based on a threshold d of the geographical distance: i and j are considered as neighbors if and only if 0 ≤ d ij ≤ d, where d ij is the distance between i and j. Moreover, all the neighbors of i are equivalent and have the same impact on the calculation of the neighborhood average

W i • Z.
Thus, the contiguity matrix W is given by:

w ij = 1 N umber of neighbors of i , if 0 ≤ d ij ≤ d 0, otherwise. (1)
To apply Moran scatterplot to the context of Velib', one has to estimate the crucial parameter d, which represents the highest distance between two neighboring Velib' stations. The choice of d has to achieve the following trade-off: On the one hand, this distance has to be small enough to let the users slightly change their trips at a local scale, and on the other hand, it has to be high to make sure that most stations have a reasonable number of neighboring stations. Velib' stations are generally close to each other and concentrated in the center of Paris and near attractive locations whereas they are distant in the suburbs.

To get in an idea about the geographical proximity between the different stations we plotted, in Figure 1, the distribution of the geographical distance to the nearest station. In Figure 1, the mean distance to the nearest station is of 216 meters. It leads us to propose that the value of parameter d is larger than 300 meters.

To better estimate the choice of d, we plotted in Figure 2 the distribution of the number of neighbors for all the Velib' stations. We tested different values for the threshold distance d (300, 400 and 500 meters). We conclude that a distance of 400 meters is reasonable as, in this case, a given station has on average about 5 neighboring stations. Moreover, with d = 400, only 4.4% of the stations do not have any neighboring station.

However, when detecting spatial outliers, the assumption that all the neighbors have the same impact on the neighborhood average may lead to missing some true spatial outliers. In the dataset described in Section 3, there are 1226 stations. As plotted in Figure 3, the capacity of the stations is highly variable between 8 and 114 bikes, with an average of about 31 bikes. As Velib' stations have different capacities, we defined the occupancy ratio in order to compare normalized bikes availability in these stations. The key idea is that two neighboring stations should have almost the same occupancy ratio if they have similar capacities. That is why the capacity of the station has to be taken into account when calculating the neighborhood occupancy average.

Improvement of Moran scatterplot using Gower's coefficient

We will replace W with a new weight matrix W also based on the degree of similarity between the station i and the corresponding neighboring stations. This new matrix will take into account the distance and also the difference of capacities between a station and its neighbors. The set N i of neighbors of station i is defined as previously by the stations with a maximal distance d from station i.

In order to measure the similarity degree between two spatial objects, the Euclidean distance is most often used. However, in our case, the use of this distance is inappropriate since the location and capacity attributes are measured on different scales. Hence, we propose to use the Gower's coefficient [START_REF] John | A general coefficient of similarity and some of its properties[END_REF] to calculate the similarity between two stations. Gower's coefficient is a similarity measure which computes the distance between two instances on each attribute k, and then aggregates all of them to finally calculate the similarity degree. Gower's similarity degree GOW ER ij between two stations i and j is defined by:

GOW ER ij = n k=1 W ijk × S ijk n k=1 W ijk (2)
where • W ijk is the weight associated to the attribute k,

• S ijk is the similarity between two stations i and j for the k th attribute, given by

S ijk = 1 - |x ik -x jk | r k
where x ik is the observable attribute k in station i and r k is a standardization for the attribute k since each attribute is of different unit.

In the context of Velib' stations, we calculate the similarity S ij of the location SD ij and capacity SC ij between two neighboring stations i and j by:

SD ij = 1 - d ij d SC ij = 1 - |Capacity i -Capacity j | Capacity max -Capacity min
where • d ij is the distance between the two stations and d is the maximal distance.

• Capacity max and Capacity min are respectively the maximal and minimal stations capacities in the neighborhood of station i.

In this definition, W ijk = W ij , where W ij is previously defined by equation ( 1).

We propose in the following to modify the construction of the contiguity matrix of weights by incorporating the spatial and non-spatial attributes and in a weighted manner in the calculation of the weights associated with neighbors. For each neighboring station j, its new weight GOW ER ij regarding the station i is given by equation ( 2).

The normalization of the contiguity matrix of weights is done per line, so for each station i, the weight of each neighboring station j is divided by the sum of the weights of all the neighboring stations of i.

Thus, the new contiguity matrix W is given by:

wij = GOW ER ij , if 0 ≤ d ij ≤ d 0, otherwise.
We applied the improved version of Moran scatterplot to detect the isolated problematic stations. Recall that these stations are defined as spatial outliers with a critical occupancy ratio. We used the same data set described in Section 3.

Moran scatterplot representation for the occupancy data of the stations at a fixed timestamp: 10 : 00 am is given in Figure 4. At this time of day, we can expect that the system is highly unbalanced, as in general in a working day a lot of trips take place in the morning around 8 : 00 am. The spatial outliers stations (almost 300 stations) are located in the upper-left and lower-right quadrants. One can notice that there are fewer points in these quadrants compared to the locations with positive correlation.

The number of detected isolated problematic stations at 10 : 00 am, depending on the allowed distance, is given in Table 1. Recall that isolated problematic stations are defined as spatial outliers with critical occupancy ratio. According to this table, there are about 50 isolated problematic stations at 10 : 00 am. The allowed distance does not have a considerable impact on the number of outliers and the isolated problematic stations. Moreover, with a local change of their trips, Velib' users can enhance the occupancy ratio of about 300 stations (spatial outliers), which represents 24.48% of Velib' stations.

Rebalancing the bikes in Velib' system

Our objective is to improve resource availability in the Velib' system by reducing the number of problematic stations. For this purpose, we propose and test in this section a new incentive method, based on a natural and ecological regulation performed by Velib' users themselves. The main idea behind the proposed method is to balance the global system by performing small changes in the trips in small local areas. The preliminary study provided in the previous section proved the presence of isolated problematic stations. Then the aim of this part is to show that around a given problematic station (in a distance smaller than 500 meters), there are many balanced stations (with an occupancy ratio around 50%), which make it possible for Velib' users to balance this problematic station by slightly changing their trips (in practice with an award, extra-time for example).

Using the dataset described in Section 3, we plot in Figure 5 the evolution of the number of current trips during the day (on Thursday 10/31/2013), in order to understand the usage of the Velib' system. One can easily identify two peaks at about 8 : 00 am and 6 : 00 pm. They clearly correspond to the trips to the offices and the return home after work, as it is a working day.

Users trips unbalance the Velib' system by making some stations problematic (almost empty or almost full). Based on the thresholds of station occupancy introduced before (10% and 90%), the current number of problematic stations is given in Figure 6. Despite the performed bike regulation using tracks, the number of problematic stations during the day remains high. The problematic stations are mainly composed of almost empty stations.

We propose in this section two incentive methods that encourage Velib' users to improve the homogeneity of the stations in terms of occupancy ratio by slightly changing their trips. In the trips dataset, let us denote by A the station where the trip begins and by B where the trip ends. The neighborhood of the station is defined by a distance less than 400 meters.

First scenario

The key idea is to change the trips as follows: For each trip, in terms of occupancy ratio,

• station A will be replaced by the busiest station in the neighborhood of A,

• station B will be replaced by the emptiest station in the neighborhood of B.

This idea can be easily implemented using a mobile application: The user begins by giving his source and destination stations, respectively A and B. Then the application will propose an alternative itinerary that enables him to save a given amount of time or money.

The departure station of this new route is the busiest station in the neighborhood of A and the arrival station is the emptiest station in the neighborhood of B. The new route is of course calculated on-line as it depends on the current state of the system. The new itinerary can be accepted by the user or not. It can be improved with other options for implementation. For example, the user could begin giving his source and destination positions, i.e. his departure and arrival locations, not Velib' stations. Then the application can propose him the best itinerary, choosing among the different stations around these positions. This proposition could be followed or not. The neighborhood is defined by a maximal distance of 400 meters. We checked in section 4 that a value of 400 meters is a suitable distance to have a reasonable number of neighboring stations. Changing the departure station of a trip engenders an extra walk of at most 400 meters to let the user reach the new optimal departure station. Notice that the actual distance is very likely to be greater than the calculated distance because it depends on the chosen path between the departure and the arrival station, whereas the calculated distance only uses the geographical positions of these stations. All the geographical distances used in this paper are Euclidean distances also called Line of Sight (LoS) distances. Given the urban planning in Paris, the real walking or riding distance is approximated by the LoS distance, for small distances. Maybe l 1 -distance is more suitable, and one could also derive exactly the stations at a fixed walking distance from a given station. The proposed methodology applies in both cases. To be realistic, the extra walk engendered by the proposed change has to be compared to the distance of the trip. For this purpose, we plotted in Figure 7 the distribution of trips lengths. Moreover, the trips having the same departure and arrival station are ignored (about 2.27% of the performed trips). The so obtained mean trip length is of 1917 meters, therefore a maximal extra walk of 400 meters can be reasonably proposed. The distribution of trips length can be explained by the fact that the first 30 minutes of the trip are free (45 minutes for the students). The mean duration of the trips is of 14.5 minutes leading to an average speed of 7.93 km/h. The proposed method is inspired by Velib' + which consists of offering users of Velib' an extra time (that can be cumulated) when they park their bike in a station having a high altitude. The main difference is that Velib'+ regulation is static: Velib'+ stations are well known and never change over time, whereas our preferred busiest and emptiest stations dynamically change. They vary during the time depending on their occupancy ratio and the occupancy ratio of their neighboring stations.

Second scenario

The second scenario consists of performing the same natural regulation proposed in the first scenario, only during the rush hours. The rush hours correspond to the trips to the offices in the morning and the return home at the end of the day. According to Figure 5, these peaks of activities are occurring in the following intervals: [7h, 10h], and [17h, 20h]. They represent 40% of the daily trips. The key idea behind this scenario is to focus on the demands of regular bikers: Changing only the trips of the regular bikers may significantly improve resource availability.

Experiments and results

Figure 8 presents the impact of the proposed incentive method using the first and the second scenario) on the number of problematic stations. Using the first scenario, The results show a clear decrease in the number of problematic stations throughout the day. The average number of problematic stations drops from 164 in real trips to only 27 by slightly modifying each trip. Starting from a relatively high number of problematic stations (almost 150), users are able to balance almost all these stations within three hours. With the second scenario, which limits the regulation to the rush hours, the number of problematic stations also drops significantly and reaches the performance of the first scenario (permanent regulation) within one hour of the regulation (at 8h and at 18h). Notice that for both scenarios, no new trips are either added or lost. The modification is done with exactly the same number of trips. The real trips are only locally modified. The obtained results confirm our intuition that the global availability of the resources in the Velib' system can be significantly improved by acting locally. This improvement would allow accepting new trips, where originally users are rejected due to a lack of bikes.

The performance of the proposed incentive method can also be measured by the number of spatial outliers in the Velib' system. They consist of stations with an occupancy ratio significantly different from the average occupancy ratio in their neighborhood. These outliers are depicted in Section 4 using Moran scatterplot. The comparison of the number of spatial outliers between the original and modified behaviors (using the first and the second scenario is given in Figure 9. With the improved user behavior, the number of spatial outliers drops significantly, which enhances the balance of the Velib' system.

In Figures 8 and9, all trips are modified according to the proposed method. It is not a realistic scenario as in real life, many users will not accept to change their departure or arrival station even if they are encouraged by a financial motivation or an extra time offered. To simulate a real-world situation, we plotted in Figure 10 the average number of problematic stations in the day under a variable collaboration rate of the users. One can see that, if only 20% of users accept to change their trips, the number of problematic stations will decrease by half. The decrease in the number of problematic stations is fast (faster than a linear decrease) which is an excellent result as we cannot expect that the majority of users will collaborate.

The number of problematic stations during the day is a good indicator to evaluate the quality of the service offered to Velib' users. However, it cannot entirely qualify service availability. For a given station, the service is considered as interrupted if there is no bike or no dock in this station. In this case, the station is said invalid or out of service. Note that this concerns just one resource: bikes or free docks. To have complete information, we plotted in Figure 11 the average duration of stations invalidity during each one-hour interval of the day, before and after the proposed improvement. One can notice that the mean duration of station invalidity has largely decreased, and likewise, the mean cumulative invalidity duration during the day has been widely improved (cf. Figure 12). According to Figure 12, at the end of the day, the mean cumulative invalidity duration of a Velib' station drops from 141 minutes to only 22 minutes using our first scenario and 68 minutes using the second scenario.

Conclusion

We studied in this paper the detection of anomalies in a spatial context. Our use case was the bike sharing system of Paris (Velib'). A Velib' station is considered as a spatial outlier if it is almost empty or almost full while the surrounding stations are globally balanced. This is due to the problem of heterogeneity in the Velib' system causing resources unavailability and user dissatisfaction.

To identify spatial outliers in this context, we used Moran scatterplot. In order to calculate the occupancy distance between two stations, we introduced a similarity weight that takes into account the geographical distance between the stations as well as the difference between their capacities. This degree of similarity is calculated using a robust distance metric called Gower's coefficient. Thereafter, we proposed and tested a new algorithm that locally improves the distribution of resources (bikes and docks) in the stations, and we experimentally validated its efficiency. The results showed that even applied only during the rush hours, the proposed algorithm improves the homogeneity of the Velib' system by reducing the number of outlier stations and the duration of unavailability of the stations during the day, which ultimately leads to the improvement of the user level of satisfaction.

  where W is a weight matrix that defines a local neighborhood around each location. The observations Z (x-axis) and W.Z (y-axis) are represented by their standardized values. Moran scatterplot contains four quadrants, corresponding to four types of spatial correlation. The upper-right and lower-left quadrants consist of the locations with positive spatial correlation: association between similar values. In the upper-right quadrant, the high values are surrounded by high neighbors values, while in the lower-left quadrant, the low values are surrounded by low neighbors values. The upper-left and lower-right quadrants incorporate the locations with negative spatial correlation: association between dissimilar values. The upper-left quadrant contains low values surrounded by high neighbors values, while the lower-right quadrant contains high values surrounded by low neighbors values. The objects located in these two quadrants are considered as spatial outliers and can be identified by the statistical test function:

Figure 1 :

 1 Figure 1: Distribution of the geographical distance to the nearest station in bike sharing system Velib', Paris, France (1226 stations, average distance 215.848 meters).

Figure 2 :

 2 Figure 2: Distribution of the number of neighbors, stations at distance less than a threshold distance d, with different values of d: 300, 400 and 500 meters in BSS Velib', Paris, France.

Figure 3 :

 3 Figure 3: Distribution of station capacity in BSS Velib', Paris, France.

Figure 4 :

 4 Figure 4: Improved Moran scatterplot based on occupancy data for Velib' bike sharing system, Paris, France, on Thursday 10/31/2013 10 : 00 am.

Figure 5 :

 5 Figure 5: Number of current trips as a function of time during one day time based on trip data for Velib' bike sharing system, Paris, France, on Thursday 10/31/2013.

Figure 6 :

 6 Figure 6: Number of problematic (empty and full), empty and full stations as a function of time during one day, based on occupancy data for Velib' bike sharing system, Paris, France on Thursday 10/31/2013.

Figure 7 :

 7 Figure 7: Distribution of the trip length based on trip data for Velib' bike sharing system, Paris, France, on Thursday 10/31/2013.

Figure 8 :

 8 Figure 8: Number of problematic stations as a function of time during one day, based on occupancy data for Velib' bike sharing system, Paris, France on Thursday 10/31/2013.

Figure 9 :

 9 Figure 9: Detected spatial outliers as a function of time during one day, based on occupancy data for Velib' bike sharing system, Paris, France on Thursday 10/31/2013.

Figure 10 :

 10 Figure 10: Average number of problematic stations as a function of the user collaboration during one day, based on occupancy data for Velib' bike sharing system, Paris, France on Thursday 10/31/2013.

Figure 11 :

 11 Figure 11: Mean duration of stations invalidity as a function of time during one day, based on occupancy data for Velib' bike sharing system, Paris, France on Thursday 10/31/2013, for the same scenarios as in Figure 10.

Figure 12 :

 12 Figure 12: Mean cumulative duration of station invalidity as a function of time during one day, based on both trip and occupancy data for Velib' bike sharing system, Paris, France on Thursday 10/31/2013, for the same scenarios as in Figure 10.

Table 1

 1 Number of detected outliers stations with the improved Moran scatterplot based on occupancy data for Velib' bike sharing system, Paris, France, on Thursday 10/31/2013 10 : 00 am

	Allowed distance Outliers Outliers with critical occupancy ratio
	300	297	53
	400	334	52
	500	339	54


Edoardo Idá
email: edoardo.ida2@unibo.it

Jean-Pierre Merlet
email: jean-pierre.merlet@inria.fr

Marco Carricato
email: marco.carricato@unibo.it

Automatic Self-Calibration of Suspended Under-Actuated Cable-Driven Parallel Robot using Incremental Measurements

Keywords: Cable-driven parallel robots, Underconstrained robots, Underactuated robots, Homing, Self-calibration

This paper focuses on the problem of the initial-pose estimation by means of proprioceptive sensors (self-calibration) of suspended under-actuated Cable-Driven Parallel Robots (CDPRs). For this class of manipulators, the initialpose estimation cannot be carried out by means of forward kinematics only, but mechanical equilibrium conditions must be considered as well. In addition, forward kinematics solution is based on cable-length measurements, but if the robot is equipped with incremental sensors cables' initial values are unknown. In this paper, the self-calibration problem is formulated as a non-linear least square optimization problem (NLLS), based on the direct geometricostatic problem, where only incremental measurements on cable lengths and on swivel pulley angles are required. In addition, a data acquisition algorithm and an initial value selection procedure for the NLLS are proposed, aiming at automatizing the self-calibration procedure. Simulations and experimental results on a 3-cable 6-degree-of-freedom robot are provided so as to prove the effectiveness of the proposed methodology.

Introduction

Cable-driven parallel robots (CDPRs) are a class of parallel manipulators which employs flexible cables instead of rigid links in order to control their end-effector (e-e) pose. A CDPR is completely-constrained when the pose of its e-e can be fully determined if actuators are locked and, thus, all cable lengths are known. On the contrary, a CDPR is under-constrained if the e-e preserves some degrees of freedom (DoFs) when the cable lengths are assigned. This may originate from the e-e being constrained by a number of cables smaller than its DoFs or from some cables being slack [START_REF] Carricato | Stability Analysis of Underconstrained Cable-Driven Parallel Robots[END_REF]. This paper focuses on suspended CDPRs, in which all cables are located above the e-e, so that their tension state only depends on the action of gravity and inertia forces. These robots are inherently under-constrained when they are underactuated, that is, the number of actuators is less than the number of variables needed to fully describe the manipulator. A growing number of studies are being conducted on under-actuated CDPRs [START_REF] Abbasnejad | Direct geometrico-static problem of underconstrained cable-driven parallel robots with n cables[END_REF][START_REF] Berti | Solving the direct geometrico-static problem of underconstrained cable-driven parallel robots by interval analysis[END_REF][START_REF] Carricato | Stability Analysis of Underconstrained Cable-Driven Parallel Robots[END_REF][START_REF] Hwang | Trajectory generation to suppress oscillations in under-constrained cable-driven parallel robots[END_REF][START_REF] Idá | Rest-to-rest trajectory planning for planar underactuated cable-driven parallel robots[END_REF][START_REF] Lin | Design and oscillation suppression control for cable-suspended robot[END_REF][START_REF] Zarei | Oscillation damping of nonlinear control systems based on the phase trajectory length concept: An experimental case study on a cable-driven parallel robot[END_REF].

A major issue in the practical use of under-actuated CDPRs is the estimation of the e-e initial-pose. When the machine is switched on in a generic start-up condition, the e-e pose is generally unknown, but its knowledge is fundamental for any subsequent operation. In order to directly measure the e-e pose, external measurement devices such as laser trackers [START_REF] Nubiola | Comparison of two calibration methods for a small industrial robot based on an optical CMM and a laser tracker[END_REF] or high-resolution cameras [START_REF] Daney | Interval method for calibration of parallel robots: Vision-based experiments[END_REF] can be employed. On the other hand, an indirect estimation of the pose can be performed by measuring some of the robot's internal joint variables, followed by the solution of the direct kinematic problem [START_REF] Merlet | Direct kinematics of cdpr with extra cable orientation sensors: The 2 and 3 cables case with perfect measurement and ideal or elastic cables, in Cable-Driven Parallel Robots[END_REF][START_REF] Pott | An algorithm for real-time forward kinematics of cable-driven parallel robots[END_REF]. When this approach is used for pose estimation in start-up conditions, the solution to this problem is sometimes referred to as self-calibration [START_REF] Borgstrom | Nims-pl: A cable-driven robot with self-calibration capabilities[END_REF] or internal-calibration [START_REF] Lau | Initial length and pose calibration for cable-driven parallel robots with relative length feedback, in Cable-Driven Parallel Robots[END_REF] of the e-e initial-pose.

In [START_REF] Borgstrom | Nims-pl: A cable-driven robot with self-calibration capabilities[END_REF], a self-calibration procedure for initial-pose estimation of a 2-DoF 4-cable over-constrained robot is proposed, which is based on cable tension and length increment measurements. The proposed method relies on the over-constrained nature of the robot. In [START_REF] Miermeister | Auto-calibration method for overconstrained cabledriven parallel robots[END_REF] a two-stage calibration procedure for generic over-constrained CDPRs is introduced, aiming at both optimizing robot static parameters and determining the initial-pose of the e-e. Ref. [START_REF] Lau | Initial length and pose calibration for cable-driven parallel robots with relative length feedback, in Cable-Driven Parallel Robots[END_REF] shows how to perform initial-pose estimation by means of a manual self-calibration procedure for over-constrained robots, only relying on cable length increment measurement.

This paper extends the method introduced in [START_REF] Lau | Initial length and pose calibration for cable-driven parallel robots with relative length feedback, in Cable-Driven Parallel Robots[END_REF] by proposing an automatic procedure to estimate the initial-pose of a generic suspended under-actuated CDPR, e.g. a 6-DoF CDPR actuated by n < 6 cables, that only relies on incremental measurements of cable length and orientation. An automatic data acquisition procedure is exploited in order to reliably and autonomously collect the information required for calibration purposes. The initial-pose estimation is formulated as a non-linear least square optimization problem (NLLS), whose initial guess is generated automatically according to the data acquisition algorithm. In Section 2, the geometrico-static model of an under-actuated CDPR is developed. In Sections 3 and 4, the NLLS optimization problem is formulated and the data acquisition algorithm employed for its solution is discussed. Finally, simulation and experiments are presented.

Geometrico-Static Modelling

Kinematics

A generic 6-DoF under-actuated CDPR consists in a mobile platform connected to the base by n < 6 cables, which are actuated by motorized winches (Fig. 1). Ox y z is an inertial frame, whereas P x y z is a mobile frame attached to the e-e. The platform pose is described by the position vector p P of P , and the rotation matrix R(φ, θ, χ) is parametrized by Tait-Bryant angles = [φ, θ, χ] T according to the x y z convention. The platform generalized coordinates are thus p = [p P , ] T ∈ R 6 .

Cables are modelled as massless and inextensible, and are guided by swivel pulleys [START_REF] Pott | Influence of pulley kinematics on cable-driven parallel robots[END_REF] to their attachment points A i , i = 1, . . . , n on the platform (Fig. 2a). G is the

p G = p P + p G = p P + R P p G , a i = p P + a i = p P + R P a i ( 1
)
The i -th swivel pulley has center C i , radius r i , and is mounted on a hinged support, whose swivel axis z i is tangent to the pulley (Fig. 2a). Vector d i defines the fixed position of point D i where the i -th cable enters the pulley's groove and which is on the z i -axis. Unit vectors i i , j i , k i , associated with an additional fixed reference frame D i x i y i z i attached in D i , describe the orientation of each swivel pulley in Ox y z. The i -th cable exit point from the pulley groove is denoted by B i and the vector ρ i = A i -B i is tangent to the pulley.

In static conditions, if friction in the pulley hinges is negligible, vectors z i and ρ i define the pulley plane, forming the swivel angle σ i with the coordinate plane x i z i .

If u i and w i are unit vectors orthogonal to k i and to each other:

w i = -sin σ i i i + cos σ i j i , u i = cos σ i i i + sin σ i j i (2)
the aforementioned static constraint can be analytically expressed as:

w i • (a i -d i ) = 0 (3)
which can be solved for σ i (p), if the pose is assigned. We define ψ i ∈ (-π, π) as the angle between u i and B i -C i (Fig. 2b). It is convenient to define two more unit vectors, orthogonal to w i and to each other:

n i = cos ψ i u i + sin ψ i k i , t i = sin ψ i u i -cos ψ i k i (4)
so that the i -th cable vector can be expressed as:

ρ i = a i -d i -r i (u i + n i ) = ρ i t i ( 5
)
where

ρ i = l i -r i (π -ψ i ) (6)
is the rectilinear cable length, and l i is the i -th total cable length, comprising the rectilinear part and the wrapped part B i D i . The tangency constraint between ρ i and the pulley can be expressed as:

n i • ρ i = 0 (7)
which can be solved for ψ i (p), if the pose (and thus the swivel angle) are known. Once ψ i (p) and σ i (p) are determined, l i (p) can be found from the geometrical constraint that the i -th cable imposes on the platform:

ρ T i ρ i -ρ i 2 = 0 (8)
where ρ i is evaluated as in Eq. ( 5) and ρ i as in Eq. ( 6).

Statics

For the purpose of this work, we will consider that the e-e is acted upon by a constant force only, e.g. gravity, which is parallel to the z axis and applied at point G. Thus, if the e-e mass is m, the gravitational acceleration is g , and p G is the skew-symmetric representation of the vector product, namely p G × = p G , the external wrench applied to the reference point

P is W = -mg [k; p G k].
If all cables are taut, the static equilibrium equations of the e-e is:

J T T -W = 0 (9)
where

T = [T 1 , • • • , T n ]
T , T i is the cable tension, and J ∈ R n×6 is the Jacobian matrix of the constraints in Eq. ( 8), whose i -th row is

J i = [t T i -t T i ã i ]. It is convenient to partition J in a n × n matrix J a and a n × (6 -n) matrix J u , namely J = [J a J u ],
and vector W as W = [W a ; W u ], so that Eq. ( 9) can be rewritten as:

J T a T -W a =0 (10)
J T u T -W u =0 (11)
The first part can be solved for T, thus yielding:

T = J -T a W a ( 12
)
and the result can be substituted in the second part, thus obtaining the static equilibrium constraint of the platform:

J T u J -T a W a -W u = 0 (13)

Forward Geometrico-Static Problem

The aim of the forward geometrico-static problem is to find the e-e pose, once the value of a suitable set of the robot's internal joint variables is known. In the case that only cable length measurements are available, it is possible to formulate such a problem by considering the n constraints in Eq. ( 8) and the 6n equations in Eq. ( 13) [START_REF] Abbasnejad | Direct geometrico-static problem of underconstrained cable-driven parallel robots with n cables[END_REF].

The resulting system of nonlinear equations is completely determined and, thus, can be numerically solved by using nonlinear solvers. Alternatively, if redundant measurements are available, one can (possibly) neglect the static constraint equations and replace them with additional kinematic constraint equations that are explicitly dependent on the measured variables [START_REF] Merlet | Direct kinematics of cdpr with extra cable orientation sensors: The 2 and 3 cables case with perfect measurement and ideal or elastic cables, in Cable-Driven Parallel Robots[END_REF]. This approach leads to a system of equations that, depending on the number and the type of redundant measurements, can be completely determined or overdetermined (the system is never underdetermined, because static equilibrium constraints can always be accounted for).

In this paper, we assume that both cable lengths and swivel pulley angles can be measured. For the sake of generality, we formulate the forward geometrico-static problem by considering both 2n kinematic constraints, derived from Eqs. ( 3) and ( 8), and the 6-n static equilibrium constraints of Eq. ( 13), thus leading to an overdetermined system of 6 + n equations in 6 unknowns (namely, the generalized pose coordinates in p). By letting:

F 1 (p) =    σ 1 (p) -σ * 1 . . . σ n (p) -σ * n    , F 2 (p) =    l 1 (p) -l * 1 . . . l n (p) -l * n    , F 3 (p) = J T u J -T a W a -W u (14)
where σ i (p) and l i (p) can be easily calculated from Eqs. ( 3) and ( 8), and σ * i and l * i are the corresponding measured values, the forward geometrico-static problem can be formulated as:

F(p) =   F 1 (p) F 2 (p) F 3 (p)   = 0 (15)
This formulation leads to a fairly simple analytical formulation of the first-order differentiation of Eq. ( 15), as described in Section 3.

Initial-Pose Estimation Problem

In this section, the initial-pose estimation problem will be formulated by extending the work presented in [START_REF] Lau | Initial length and pose calibration for cable-driven parallel robots with relative length feedback, in Cable-Driven Parallel Robots[END_REF]. If the under-actuated CDPR is equipped with incremental measurement devices on motors and swivel axes, i.e. incremental encoders, cable lengths and swivel angles at a generic pose p i can be measured relatively to the initial values σ 0 i and l 0 i at pose p 0 , namely:

σ * i = σ 0 i + ∆σ * i ( 16
)
l * i = l 0 i + ∆l * i ( 17
)
While ∆σ * i and ∆l * i are measures provided by the encoders, σ 0 i and l 0 i are generally unknown and are the objective of the self-calibration procedure. The direct geometricostatic problem in Eq. ( 15) can thus be expressed as:

F(σ 0 , l 0 , p) = 0 ( 18
)
where

σ 0 = [σ 0 1 , • • • , σ 0 n ] T and l 0 = [l 0 1 , • • • , l 0 n ] T
. This problem has 6 + n equations and 6 + 2n unknowns (σ 0 , l 0 , p), and it has generally an infinite number of solutions.

However, by assuming that λ different measurement sets are available, that is:

σ * i ,k = σ 0 i + ∆σ * i ,k l * i ,k = l 0 i + ∆l * i ,k k = 1, • • • , λ (19)
the following system of equations is obtained:

G(X) = G(σ 0 , l 0 , p 1 , • • • , p λ ) = [F(σ 0 , l 0 , p 1 ); • • • ; F(σ 0 , l 0 , p λ )] = 0 (20)
where X = [σ 0 ; l 0 ; p 1 ;

• • • ; p λ ] ∈ R 6λ+2n .
The system (20) has a total of λ(6 + n) equations and 6λ + 2n unknowns. Thus, if λ > 2, the initial-pose estimation problem is overdetermined and can be formulated as a non-linear least-square optimization:

X opt = arg min X G(X) 2 (21)
This problem can be solved by employing numerical techniques, such as the Levenberg-Marquardt algorithm. The efficient solution of Eq. ( 21) relies on a reasonable initial solution guess X g uess (see Section 4), and an analytical formulation of the Jacobian matrix in Eq.(20). While the former is fundamental for both the solution accuracy and the algorithm rapidity, the latter is critical only in terms of computational time.

The Jacobian ∂G/∂X can be expressed analitically as:

∂G ∂X =                    -I n×n 0 n×n ∂σ(p 1 )/∂p 0 n×6 • • • 0 n×6 0 n×n -I n×n ∂l(p 1 )/∂p 0 n×6 • • • 0 n×6 0 (6-n)×n 0 (6-n)×n ∂F 3 (p 1 )/∂p 0 (6-n)×6 • • • 0 (6-n)×6 -I n×n 0 n×n 0 n×6 ∂σ(p 2 )/∂p • • • 0 n×6 0 n×n -I n×n 0 n×6 ∂l(p 2 )/∂p • • • 0 n×6 0 (6-n)×n 0 (6-n)×n 0 (6-n)×6 ∂F 3 (p 2 )/∂p • • • 0 (6-n)×6 . . . . . . . . . . . . . . . . . . -I n×n 0 n×n 0 n×6 0 n×6 • • • ∂σ(p λ )/∂p 0 n×n -I n×n 0 n×6 0 n×6 • • • ∂l(p λ )/∂p 0 (6-n)×n 0 (6-n)×n 0 (6-n)×6 0 (6-n)×6 • • • ∂F 3 (p λ )/∂p                    (22)
where I n×n and 0 n×n are the n×n identity and zero matrices, 3),( 7) and ( 8) yields:

σ(p k ) = [σ 1 (p k ); • • • ; σ n (p k )] and l(p k ) = [l 1 (p k ); • • • ; l n (p k )]. Manipulating Eqs. (
∂σ i ∂p = w T i u T i (a i -d i ) - w T i ã i H u T i (a i -d i ) , ∂l i ∂p = t T i -t T i ã i H (23)
where, respectively denoting cos x and sin x by c x and s x :

H =   1 0 s θ 0 c φ -s φ c θ 0 s φ c φ c θ   (24)
Finally, by considering that:

∂J T i ∂p = ∂t i /∂p ã i (∂t i /∂p) -ti (∂a i /∂p) , ∂W ∂p = -mg 0 3×1 kr H (25)
and:

∂t i ∂p = w i sin ψ i ∂σ i ∂p + n i ∂ψ i ∂p , ∂a i ∂p = 0 3×3 -ã i H ( 26
)
∂ψ i ∂p = n T i ρ i - n T i ã i H ρ i ( 27
)
∂F 3 /∂p can be evaluated according to the partition defined in Section 2.2 as:

∂F 3 ∂p = ∂J T u ∂p -J T u J -T a ∂J T a ∂p J -T a W a + J T u J -T a ∂W a ∂p - ∂W u ∂p (28)

Data Acquisition Algorithm

It is beyond the scope of this work to determine an optimal data acquisition algorithm. However, a practical one, which enables autonomous and safe operation of the CDPR during calibration, is provided hereafter. For this aim, cable tensions or alternatively motor torques are assumed to be measurable or at least estimated, and actively controlled by a suitable feedback system.

In the instant the robot is switched on, its pose is generally unknown. It is possibly unsafe to start the self-calibration process in this configuration. Then, it is useful to pre-determine a safe start configuration, in which every cable is taut and sufficiently long, so that it can be coiled and uncoiled, and the e-e may attain different poses. By assigning a start cable tension vector T 0 , the static problem (9) can be solved for p as a non-linear system of six equations in six unknowns. Because of the non-linearity of the problem, a finite set of real solutions can be determined: this calculation may be done off-line and just once, i.e. during robot parameter calibration. Only stable solutions [START_REF] Carricato | Stability Analysis of Underconstrained Cable-Driven Parallel Robots[END_REF] among the possibly many available should be considered. Additionally, T 0 may be selected so that only one stable solution exist.

In the following, we will consider a start cable tension vector leading to a unique stable solution of problem [START_REF] Lin | Design and oscillation suppression control for cable-suspended robot[END_REF]. Accordingly, a (computed) start pose p 0,comp is unambiguously determined, as well as start cable lengths l 0 comp and swivel angles σ 0 comp . The real start pose p 0 attained by the CDPR can be fairly different from the ideal one p 0,comp , and its determination is the aim of the self-calibration procedure. A maximum cable tension vector T m should be set as well. The data acquisition algorithm objective is to ensure that every DoF of the e-e is varied during measurements, so that problem (21) is always well conditioned. The procedure workflow can be summarized as follows.

1. Start phase: command the CDPR so that cable tensions (or motor torques) quasistatically reach the assigned value T 0 . When T 0 is reached and static conditions are attained 3 , the j -th actuator is assigned an incremental cable tension (or motor torque) set-point, starting from j = 1. The change in a single actuator setpoint ensures that the pose of the end-effector is different at any iteration, thus being effective, as well as practical and easy to implement; 2. Tensioning phase: quasi-statically move the CDPR by assigning λ/(2n) positive increments of magnitude ∆T = 2n(T j ,m -T j ,0 )/λ to the tension set-point of the j -th actuator, namely T j ,k = T j ,k-1 + ∆T , where T j ,0 is the j -th component of T 0 . After each assignment k, the CDPR e-e could possibly oscillate during the transition. When static conditions are attained, record measurements ∆σ * i ,k and ∆l * i ,k , for i = 1, • • • , n; 3. Detensioning phase: assign λ/(2n) negative increments of magnitude ∆T to the tension set-point of the j -th actuator, namely T j ,k = T j ,k-1 -∆T . When static conditions are attained, record measurements ∆σ * i ,k and ∆l * i ,k , for i = 1, • • • , n. During the detensioning phase, the robot follows exactly the same cable tension (motor torques) set-points as in the tensioning phase: on a real machine, due to repeatability errors, this could lead to different cable lengths and swivel angles, which are possibly useful in the calibration procedure in order to minimize the repeatability error of the robot. λ/n measurement sets are thus obtained by varying a single actuator set-point; 3 Notice that cable tensions are only used to lead the platform to poses where the robot is stable and kinematic measures can be accurately performed. They have no role in the optimization problem solution, since they do not appear as variables in Eq. (21). They may be affected by appreciable errors without compromising the procedure, whereas platform stability plays a key role in the data acquisition process. The initial guess for the solution of problem ( 21) is computed as:

Table 1: Actuation unit properties

i d i [m] r i [m] P a i [m] x i y i z i 1 -2.
X g uess = σ 0 comp ; l 0 comp ; p 1,comp ; • • • ; p λ,comp (29)
where p k,comp , k = 1, . . . , λ, can be evaluated by solving the static problem [START_REF] Lin | Design and oscillation suppression control for cable-suspended robot[END_REF] with assigned tension T k . Finally, by employing

X g uess , ∆σ * k = [∆σ * 1,k • • • ∆σ * n,k ] T and ∆l * k = [∆l * 1,k • • • ∆l * n,k ] T for k = 1, • • • , λ,
it is possible to determine X opt as a solution of (21). Ideally, σ 0 opt and l 0 opt should converge to σ 0 and l 0 , respectively.

Simulations and Experimental Results

A simulation example, developed in MATLAB®, aims at testing the efficiency of the optimization routine as a function of λ and the different formulations of the Jacobian. Experimentation on a prototype shows the application of the method in a realword scenario. The geometrical properties of a 3-cable 6-DoF under-actuated CDPR used for both simulation and experimentation are summarized in Table 1, where e x = [1 0 0]] T , e y = [0 1 0] T and e z = [0 0 1] T are elements of the canonical basis of SO [START_REF] Borgstrom | Nims-pl: A cable-driven robot with self-calibration capabilities[END_REF]. The platform mass is m = 8 Kg, and P p G = 0 0 0.182 T m.

Simulation results

If the start tension vector T 0 = 40 40 40 T N is assigned, the only stable solution of the static problem ( 9) is p 0 = -1.363 0.963 -0.460 -0.056 -0.069 -0.557 T m, l 0 = 1.470 1.509 1.543 T m, and σ 0 = 0.738 0.607 -0.799 T rad. Since the simulation is based on a perfect data acquisition process, the above (theoretical) values of p 0 , l 0 and σ 0 are the expected output of the optimization algorithm.

According to the procedure proposed in Section 4, several simulations were performed with different values of λ, i.e. λ = 12 ÷ 372. The initial guess X g uess for the optimization solver was generated by considering a perturbed start tension T 0,per t = 46.9 23.9 32.6

T N. In this case, the only stable solution of the static problem (9) yields p 0,comp = -1.648 0.740 -0.732 -0.246 0.179 -0.523 T m, In the case of ideal measurements, a solution X opt to the problem (21) is always found that is not only a minimum, but a real numerical zero. Solutions were computed in Matlab by using the lsqnonlin solver. Figure 3 shows the computational time and the norm of the error e = p 0 -p 0,opt for a variable number λ of measurement sets, in the case a finite-difference estimation (FD) or an analytical formulation (J) of the Jacobian. It is clear that the analytical formulation is crucial for the algorithm timeperformance, whereas there is no noticeable difference in the numerical precision achieved by the two numerical formulations. Simulations also show that there are no advantages in acquiring an increasing number of data points for the initial-pose estimation. The norm of error e is consistent with the numerical tolerance.

Experimental results

The proposed data acquisition strategy and calibration method was tested on a prototype. Swivel angles were measured by 16-bit incremental encoders, mounted directly on the swivel axes of pulleys, whereas cable lengths were estimated by using 20-bit incremental encoders on each motor axis and a kinematic model of the winch. Swivel pulleys were manufactured by FDM technology, thus limited, but not negligible, errors in their geometry and elasticity exist. Cables were coiled on IPAnema winches [START_REF] Pott | IPAnema: A family of Cable-Driven Parallel Robots for Industrial Applications[END_REF]. Clearance and elasticity in the winch components, as well as cable elasticity itself, are possible sources of error in the estimation of cable lengths.

In order to measure the real pose of the platform during experiments, 8 VICON Motion Capture Systems cameras were employed to track the position of 4 markers mounted on the robot platform.The accuracy of the measure is ±0.2mm for each dimension of the marker (x, y, z), according to manufacturer specifications. In the end, the position of the reference point and the platform orientation were reconstructed from the recorded position of each marker. Because of the lack of force sensors in the robot set-up, motor torques were employed instead of cable tensions for the implementation of the algorithm presented in Section 4. The start tension vector and maximum cable tensions were set to T 0 = 40.0 40.0 40.0

T N and T m = 80.0 80.0 80.0 T N, respectively, and converted in motor torques according to static equilibrium of the cable transmissions. In the end, λ = 60 was chosen as a trade-off between accuracy of the initial-pose estimation and data-acquisition speed.

The results of five experiments are reported in Table 2, where p 0 is the real starting pose, as measured by the motion tracking system, p 0,opt is the estimated starting pose resulting from the solution of the problem (21), e p = p 0p 0,opt is the reference position error norm and e = 0 -0,opt is the norm of the error of the orientation parameters. Positions are expressed in millimeters and angles in degrees. The execution of the calibration procedure required, on average, 4 min for the data acquisition procedure, and 2.5 s for the initial-pose estimation by using the analytical formulation of the Jacobian (16 min by using a finite-difference estimation).

During experiments, it was observed that the orientation of the swivel pulley axes plays a crucial role for the conditioning of problem (21). Pulley orientations were set in order to achieve the best possible results with the robot architecture at hand, but are not optimal. Nonetheless, the results are satisfactory by considering the modelling simplification and the hardware at hand.

Conclusions

This paper presented an automatic self-calibration method for suspended underactuated CDPRs, extending the method proposed in [START_REF] Lau | Initial length and pose calibration for cable-driven parallel robots with relative length feedback, in Cable-Driven Parallel Robots[END_REF]. An analytical formulation of the problem was developed and an autonomous procedure to calibrate the initialpose of the CDPR by means of incremental measurements of cable lengths and swivel angles was shown. Experimental results on a prototype show that the application of the proposed method is promising, but it needs refining in order to achieve a higher accuracy. In particular, the optimal orientation of the swivel pulleys and the effect of uncertainties in the measurements will be investigated in the future.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Geometry of the cable-driven robot

T-e y e z -e x 4 .

 4 Increment phase: j = j + 1; if j ≤ n then go to point 2, otherwise the algorithm is finished because λ measurement sets have been recorded.

Fig. 3 :

 3 Fig. 3: Simulation results comparing FD and J formulations

i p 0 1 T

 01 [mm, • ] p 0,opt [mm, • ] e p [mm] e [ • ]1 -1314 900 -355 -5.0 -4.8 -32.0 T -1324 906 -361 -4.5 -4.9 -32.

Table 2 :

 2 Experiments results


Remerciements/Acknowledgements

Au Docteur Laura RICHERT Un immense merci pour m'avoir proposé ce projet passionnant, de m'avoir guidé, de ta disponibilité, de ta bienveillance, ainsi que de m'avoir fait confiance. C'est un honneur pour moi d'avoir réalisé ce travail sous ta direction. Cela a été un très grand plaisir de travailler ensemble ces dernières années et je suis extrêmement ravi que nous puissions continuer les années à venir sur de nouveaux projets.

Au Professeur Rodolphe THIEBAUT

Merci pour tes conseils avisés, de ton soutien depuis mon arrivée dans l'équipe et ton enthousiasme communicatif pour la recherche. Sois assuré de toute ma reconnaissance, c'est à la fois un honneur, un plaisir et une chance de pourvoir continuer à travailler avec toi.

A Béhazine COMBADIERE

Je vous remercie de l'honneur que vous me faites en venant à Bordeaux pour présider le jury de cette thèse. Soyez assurée de toute ma gratitude pour l'intérêt que vous portez à ce travail. C'est pour moi un grand plaisir que vous jugiez ce travail et que vous y apportiez votre oeil expert dans le champ de l'immunité vaccinale.

To Professor Marie REILLY

I am very grateful that you accepted to judge my thesis as one of the main examiners. Your great experience in epidemiology and biostatistic will certainly provide me elements to improve my thinking on the subject. Thank you for coming to Bordeaux for the defense.

To Professor Niel HENS

I would like to express my great gratitude and for accepting to judge my thesis as one of the main examiners. Your great experience in mathematical and statistical methods in infectious disease epidemiology will be very relevant to get your thoughts on this work.

To Professor Sarah WALKER

It is a great honor that you accepted to judge my thesis and thank you very much for travelling to Bordeaux for the defense. This will be a real pleasure to meet you again, this is the opportunity for me to thank you for the fruitful collaboration on CHAPAS 3 during my time in London. It is an honor for me to benefit once again from your great knowledge in infectious diseases and statistics. I am looking forward to exchanging with you on this subject.

Je remercie les personnes qui ont contribués à ce travail :

• Equipe SISTM : Boris Hejblum pour son aide précieuse sur les aspects statistiques et sa disponibilité, Lise Mandigny, Chloé Pasin

• Immunitoring plateforme du VRI : Aurélie Wiedemann, Christine Lacabaratz, au Pr Yves Levy, pour leur expertise immunologique qui a enrichi les réflexions dans le cadre de cette thèse, pour leur accueil à l'hôpital Mondor et la découverte des ICS

• L'ensemble de l'équipe PREVAC, de Bordeaux, Paris, en Guinée et d'ailleurs, qui aura été un formidable projet fil rouge de ces trois dernières années. Merci à Moses Badio pour sa collaboration dans la rédaction du protocole.

• Aux participants des essais HVTN 068, LIGHT et VRI01, et PREVAC, ainsi qu'à leurs investigateurs sans lesquels les exemples d'application dans cette thèse n'auraient pas été possibles.

• Merci également à toutes les autres personnes qui m'ont aidées de près ou de loin pour réaliser ce projet.

Merci à toutes les personnes qui ont contribué à ma formation, et mon épanouissement dans mon début d'activité professionnelle (et non citées précédemment) :

• Au Professeur Geneviève Chêne, pour tout ce que vous m'avez apporté, de votre soutien tout au long de mon parcours professionnel et votre confiance, et en vous souhaitant tous mes voeux de réussite pour vos nouvelles fonctions à la tête de Santé Publique France,

• A toutes les personnes qui ont contribuées à ma formation au cours de mon internat, du CPIAS Nouvelle Aquitaine, Simone Mathoulin-Pelissier et l'équipe de l'UREC de l'Institut Bergonié, Renaud Bessellère de l'HIA Robert Picqué, Robin Callard, Nigel Klein et l'équipe de l'Institute of Child Health de London, et Christophe Tzourio et l'équipe I-Share/Espace santé étudiant,

• Aux enseignants de l'ISPED, • A l'ensemble des nombreux collègues de l'équipe SISTM et l'équipe de l'USMR, du CMG-EC et de la plateforme EUCLID/F-CRIN avec qui je suis ravis de continuer l'aventure.

A tous les internes de santé publique de Bordeaux, DOM-TOM et d'ailleurs croisés au cours de mon internat et les thésards de l'ISPED pour notamment la bonne ambiance qui a régnée pendant toutes ces années que nous avons partagées. Aux collègues et amis croisés lors de mes différentes activités associatives, à l'AIHB, l'AISPB, Jeunes Médecins, et Aquithealth/WDMH pour le travail commun, le partage d'expériences professionnelles, et leur amitié. Mes collègues et amis AHU Eric, Sébastien, Florence. To Diana for sharing my office and for the support.

A tous mes amis, de Poitiers, Bordeaux et d'ailleurs, notamment Nath, Jean-Marc, Guilhem, Dimso, Benoit, Remi, Lassina, Duc, Kevin, Pierre-Antoine, David, et Thomas.

A ma famille, pour leur soutien bénéfique tout au long de ces années.

Résumé

Les essais cliniques de vaccins posent des défis méthodologiques particuliers, principalement liés à la spécificité du développement des vaccins, les schémas des essais cliniques, de l'absence d'un corrélat de protection validé et de la complexité des techniques immunologiques évaluant l'immunogénicité des candidats vaccins. Celles-ci nécessitent des recherches méthodologiques pour définir les méthodes les plus appropriées. Cette thèse porte sur la recherche méthodologique visant à optimiser les méthodes utilisées dans le développement clinique des vaccins, notamment pour proposer et développer des méthodes statistiques de modélisation de l'immunogénicité, en prenant comme exemple les essais cliniques des vaccins contre le VIH et le virus Ebola. Nous avons d'abord étudié la dynamique des réponses immunitaires après la vaccination et montré qu'il faudrait tenir compte des temps de mesure précoce dans les futurs essais cliniques pour mieux comprendre le rôle des lymphocytes T CD4 auxiliaires et évaluer leur rôle prédictif dans la réponse immunitaire aux vaccins. Ensuite, nous avons développé une nouvelle approche de modélisation bivariée pour l'analyse de la réponse cellulaire mesurée par la technique de cytométrie en flux de marquage des cytokines intracellulaires (ICS). Cette nouvelle méthode a montré de très bonnes performances statistiques et devrait devenir la nouvelle méthode statistique standard pour les analyses ICS dans les essais vaccinaux. Ce travail aura un impact direct sur l'évaluation de la réponse du ICS dans les essais cliniques de vaccins. En ce qui concerne la réponse humorale, nous avons montré qu'il subsiste des incertitudes importantes sur les déterminants de la réponse anticorps après une vaccination préventive contre le virus Ebola. Cela met l'accent sur l'intérêt d'harmoniser les méthodes de mesure et les schémas d'étude. De plus, il indique le besoin de mettre en place des essais cliniques randomisés à bras multiples sur le vaccin contre le virus Ebola pour une comparaison précise de l'immunogénicité entre les différentes stratégies vaccinales. Nous avons présenté enfin la méthodologie d'un essai randomisé de phase 2 contre le virus Ebola évaluant trois stratégies vaccinales dans quatre pays d'Afrique de l'Ouest, et plus particulièrement pour finir une réflexion méthodologique et éthique sur la question de l'inclusion des personnels de l'essai dans un essai clinique vaccinale contre le virus Ebola en période non-épidémique. Les méthodes développées dans le cadre de cette thèse contribueront à améliorer la conception et l'analyse des futurs essais vaccinaux, et pourraient également être transposées plus largement à d'autres domaines de recherche. Mots clés : vaccin ; immunogénicité ; essai clinique ; modélisation ; VIH ; Ebola Summary Specific methodological challenges exist in vaccine clinical trials, due principally to specificities of vaccine development, clinical trial design, absence of validate correlate of protection, and complexities of new immunological assays for evaluating immunogenicity of vaccine candidates. These require methodological research to define the most appropriate methods. This thesis focuses on methodological research to optimize methods used in the clinical development of vaccines, especially to propose and develop statistical methods to model immunogenicity, using HIV and Ebola vaccine clinical trials as an example. We first investigated the dynamics of the immune responses post-vaccination and showed that early sampling time points should be considered in future clinical trials to better understand the role of the early CD4 helper T cells and to evaluate their predictive role in the immune response to vaccines. Then, we developed a new bivariate modelling approach for the analysis of the cellular immune response (assessed by intracellular cytokine staining, ICS) that showed good statistical performances and should become the new statistical standard method for ICS analyses in vaccine trials. This work will have a direct impact on the assessment on the ICS response in vaccine clinical trials. Regarding the humoral response, we showed that there are still significant uncertainties in the determinants of the antibody response after preventive vaccination against Ebola virus disease. This emphasizes the interest of harmonizing measurement methods and study designs. Furthermore, it indicates the need of randomized multi arm Ebola vaccine trials for accurate comparison of immunogenicity between different vaccine strategies. Finally, we presented the methodology of an international randomized phase 2 trial against Ebola, and in particular a methodological and ethical reflection related to the enrollment of study personnel in Ebola vaccine trial in a nonepidemic context. Methods developed in this thesis will contribute to improve the design and analysis of future vaccine trials, and also could be transposable more widely to other research domains. Key words: statistical modeling; immunogenicity; vaccine; clinical trials; HIV; Ebola

Unité de recherché / Research unit

Bordeaux Population Health Research Center -UMR-S U1219 BPH

  ISPED, Institut de Santé Publique, d'Épidémiologie et de Développement Université de Bordeaux, case 11 146, rue Léo Saignat, 33076 Bordeaux Cedex, France

Scientific production related to the subject of the thesis

Communications in conferences with a peer-review committee

Awards related to the subject to the thesis

Award by the Scientific Committee of the 13rd « Conférence Francophone d'Epidémiologie Clinique » (EPICLIN), Toulouse, France, 15-17 May 2019 for the oral communication « Evaluation de l'immunogénicité cellulaire dans les essais cliniques vaccinaux : un modèle bivariée pour mieux prendre en compte la réponse non spécifique ».

Table of contents


Alexandre Routier

Arnaud Marcoux

Mauricio Diaz Melo

Jérémy Guillon

Jorge Samper-González

Junhao Wen

Simona Bottani

Alexis Guyot

Elina Thibeau-Sutre

Marc Teichmann

Marie-Odile Habert

Stanley Durrleman

Ninon Burgos

Olivier Colliot

New advances in the Clinica software platform for clinical neuroimaging studies

Introduction

Carrying out clinical neuroscience studies usually involves many data analysis steps including image preprocessing, extraction of image-derived measurements and statistical analysis. Most of these steps also need to be performed when developing machine learning methods for neuroimaging. Clinica aims to automate the processing and statistical analysis of neuroimaging data and ease the development of machine learning approaches. We present the new functionalities and improvements made to Clinica, including framework upgrades, continuous integration and new pipelines (structural connectome, surface-based PET analysis, fMRI preprocessing, machine learning).

Methods

Clinica is based on the Nipype library [START_REF] Gorgolewski | Nipype: A Flexible, Lightweight and Extensible Neuroimaging Data Processing Framework in Python[END_REF] and relies on tools from the neuroimaging community to build pipelines (FreeSurfer [START_REF] Fischl | FreeSurfer[END_REF], FSL (Jenkinson et al., 2012), SPM (Frackowiak et al., 1997), ANTs [START_REF] Avants | Spatial and Anatomical Regularization of SVM: A General Framework for Neuroimaging Data[END_REF], MRtrix3 [START_REF] Tournier | MRtrix: Diffusion tractography in crossing fiber regions[END_REF], and PETPVC [START_REF] Thomas | PETPVC: a toolbox for performing partial volume correction techniques in positron emission tomography[END_REF]) or build machine learning modules (scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]). The only assumption on the raw data is that they follow the BIDS format [START_REF] Gorgolewski | The brain imaging data structure, a format for organizing and describing outputs of neuroimaging experiments[END_REF]. New functionalities of Clinica can be divided into three main parts.

First, new pipelines were developed in this release of Clinica:

fmri-preprocessing: This pipeline corrects fMRI datasets for head motion, slice-timing, susceptibility distortion and mainly relies on SPM. It also performs coregistration with T1, spatial normalization into MNI space and spatial smoothing.

dwi-connectome: This pipeline computes a structural connectome, i.e. a weighted graph encoding anatomical connections between a set of brain regions. It relies on the constrained spherical deconvolution diffusion model implemented in MRtrix [START_REF] Tournier | MRtrix: Diffusion tractography in crossing fiber regions[END_REF], performs probabilistic tractography and computes a connectome using the Desikan & Destrieux atlases from FreeSurfer.

pet-surface (Marcoux et al., 2018): This pipeline enables the analysis of PET data on the cortical surface. It includes co-registration of PET and T1 MR images, intensity normalization, partial volume correction, robust projection of the PET signal onto the cortical surface extracted from the coregistered T1 image, and spatial normalization to a template.

machinelearning-prepare-spatial-svm [START_REF] Avants | Spatial and Anatomical Regularization of SVM: A General Framework for Neuroimaging Data[END_REF]: This pipeline prepares T1 MRI and PET data for performing classification with an SVM that accounts for the spatial and anatomical structure of neuroimaging data, leading to more regular and anatomically interpretable classification results.

The pipelines providing image-derived measurements are listed in Fig 1 . Secondly, the core of Clinica was consolidated and upgraded to use Python 3 and Nipype 1. Thirdly, we developed a robust continuous integration framework, addressing some specific challenges raised by medical image processing. For instance, pipelines including stochastic models make it difficult to reproduce any experiment with machine precision. To solve this problem, we implemented comparison criteria adapted for each type of data and processing pipelines (correlation coefficient, metric to control the quantity of voxels exceeding certain relative difference thresholds). This allows developers to monitor the consistency of their results with respect to their developments or after a software package upgrade.

Results

The new pipelines expand the scope of modalities covered by Clinica and offer new functionalities. Fig 2 illustrates how pipelines can be chained to carry out a study, from the conversion of the raw imaging data to the statistical analysis.

Clinica has already been used for clinical research studies [START_REF] Routier | Structural, Microstructural, and Metabolic Alterations in Primary Progressive Aphasia Variants[END_REF][START_REF] Wen | Neurite density is reduced in the presymptomatic phase of C9orf72 disease[END_REF]) and to develop machine learning frameworks for Alzheimer's disease diagnosis [START_REF] Samper-Gonz Ĺez | Reproducible evaluation of classification methods in Alzheimer's disease: Framework and application to MRI and PET data[END_REF][START_REF] Toro | Supervoxels-Based Histon as a New Alzheimer's Disease Imaging Biomarker[END_REF]. This is in line with the target audience of Clinica, namely neuroscientists and clinicians conducting neuroimaging studies, and researchers in machine learning.

Conclusions

  Clinica (www.clinica.run) is an open source software platform for reproducible clinical neuroscience studies. It aims to make clinical research studies easier and pursues the community effort of reproducibility.

Figures

  Figures