Maurice Aron
email: maurice.aron@ifsttar.fr.

Régine Seidowsky

Simon Cohen

Ex-ante assessment of a speed limit reducing operation -A data-driven approach

Keywords: Injury accidents, prediction, ex-ante assessment, speed limit, traffic conditions, power model

Reducing speed limits often leads to a decrease in the injury and fatality accident counts. This depends on the type of accident, on speed limits, on compliance and on traffic conditions, before and after implementation of the speed limit. In France, the speed limit is 130 km/h and 110 km/h on interurban and urban motorways respectively. The safety assessment of a speed limit reducing operation before its implementation is the subject of this paper. Such an assessment was required before the decision to reduce speed limits in an urban motorway network in the north of France. The safety assessment process consists of four steps: first, accident analysis for the prior period, in order to count the speed-related or density-related accidents; second, analysis of the empirical six-minute average speed and density before and after the introduction of the speed limit reduction; third, predicting new traffic conditions; fourth, using relationships linking the accident count to changes in speed and density. Results predict a decrease in the accident count, which is sensitive to various possible assumptions.

Résumé

Les décisions de réduction de la limite de vitesse ont souvent conduit à une baisse du nombre des accidents corporels, le niveau de la baisse dépendant des vitesses, limites et pratiquées, avant et après l'introduction de la nouvelle mesure, des conditions de trafic, des types des accidents. Cette contribution vise à présenter l'évaluation a priori de l'impact sur le nombre des accidents corporels, d'une réduction de la limite de vitesse de 20 km/h sur huit autoroutes autour de Lille au nord de la France. Quatre étapes sont nécessaires pour établir une prévision du nombre des accidents : l'analyse des accidents survenus avant l'application de la mesure, l'analyse des conditions de trafic (vitesse, concentration) avant, la prévision des conditions de trafic futures, la spécification et l'utilisation des liens accidents-vitesse et accidentsconcentration. Les résultats prédisent une réduction du nombre des accidents, réduction qui est sensible aux différentes hypothèses que l'on peut formuler. Mots-clé: Accidents corporels ; prévision ; évaluation a priori ; limites de vitesse ; conditions de circulation.

Introduction.

Driving at an inappropriate speed contributes to the occurrence of road accidents. A number of authors have used traffic conditions and accident databases to analyse the relationship between accidents and speed [START_REF] Golob | Freeway safety as a function of the traffic flow[END_REF], [START_REF] Abdel-Aty | Identifying crash propensity using specific traffic speed conditions[END_REF]; accidents and Volume/Capacity ratio [START_REF] Lord | Modeling Crash-Flow-Density and Crash-Flow-V/C Ratio Relationships for Rural and Urban Freeway Segments[END_REF]; and accidents and traffic states [START_REF] Yeo | Impact of traffic states on freeway crash involvement rates[END_REF]. The relationship between speed limits and accidents has been analysed from national accident databases before and after the setting of new speed limits [START_REF] Nilsson | Traffic Safety Dimensions and the Power Model to Describe the Effect of Speed on Safety[END_REF], (Elvik, 2009[START_REF] Elvik | Are-parameterisation of the Power Model of the relationship between the speed of traffic and the number of accidents and accident victims[END_REF], [START_REF] Pauw | Safety effects of reducing the speed limit from 90 km/h to 70 km/h Accident Analysis & Prevention[END_REF], who suggested and used the power model, relating the after/before ratio of accident counts to the speed limit ratio.

The safety impact of a speed limit reduction (prior to its implementation) is assessed in this paper. The safety assessment process consists of four steps that are introduced in Section 2. Such an assessment was required before making the decision to reduce the speed limits in part of the ALLEGRO urban motorway network in the North of France -the site, traffic and accidents data are presented in Section 3. The method requires some relationships between accidents and traffic conditions: their calibration is given in Section 4, based on data coming from another motorway set, the Marius network in the South of France. Section 5 gives the results of the four steps for the ALLEGRO network, and Section 6 outlines the validation of this approach and the links with a future ex-post assessment, which has not yet been carried out. Concluding remarks are given at the end.

Method for predicting the injury or fatality accident count.

Making a prediction is a twofold process:

• Identifying what things should not change. There is no point in extrapolating them since either they remain constant or , they are excluded (subtracted) from the analysis and then added to the predictions. • Identifying things that will change, formulating a model for them, calibrating this model and then using it.

In both cases hypotheses are made, which should be either validated or commonly accepted in the community. The prediction process consists of four steps: 2.1. Accident analysis [START_REF] Golob | Freeway safety as a function of the traffic flow[END_REF] have shown strong links between accident counts (by type of accident) and traffic conditions. The first step of our method consists in simplifying his approach by classifying the accidents into just three types and reducing their links to traffic conditions to two dimensions (speed and density): Type 1. Some accidents are mainly due to the driver suffering from fatigue or falling asleep, to alcohol, or to a vehicle failure. Let us consider here that the number of such accidents is independent of any speed limit. Type 2. Other accidents are linked to a high speed (run off the road, rollovers); they occur in a low traffic flow and concern a single vehicle. Let us consider that changes in speed affect only those single-vehicle accidents. Type 3. The remaining accidents concern several vehicles (side collisions, rear-end collisions, multi-vehicle collisions) and are linked to a high relative speed, a high traffic density or a lane change. Let us consider that these accidents, concerning two or more vehicles, are affected only by the changes in traffic density.

Empirical speed analysis

Empirical speeds are analysed in two components:

• Fundamental Diagrams (FD), which link speed to traffic flow; FDs vary according to the speed limits.

• A deviation round the FD, because the FD is valid only on average and not at every single moment. These deviations are different at each six-minute period. Let us consider that they do not depend on the speed limit.

Traffic conditions prediction

It is the traffic authorities which decide speed limits (from safety, environmental or traffic efficiency points of view), but compliance is not perfect. Our task is thus twofold: to predict what will be and what ideally should be.

Predicting what will be

Let V FD (q) be the FD giving the speed corresponding to the traffic flow q. The empirical speed is split here into two components, the first given by the FD, the second being the deviation from the FD. For every six-minute period i, let q i , v i ,k i be the empirical traffic flow, speed and density:

( ) ( ) ( ) FD FD i ii i speed _ limit speed _ limit v vq v vq = +- (1
) The relationship q i =k i v i implies that an average six-minute speed decrease leads to an increase in density k i.. For the non-congested branch of the FD, provided the new density (after the speed limit implementation) remains lower than the critical density, the traffic demand q i will go through. If it is not the case, the period passes into congestion, with a queue that needs to be managed. However, in the simulation study carried out on the same subject [START_REF] Cohen | Evaluation a priori de l'impact des modifications des vitesses limites sur le réseau autoroutier de l'agglomération lilloise[END_REF] has shown that the increase in congestion was negligible. In addition, capacity was virtually maintained (see the FDs Fig. 4), except for the slow lane, when the new speed limit is 90 km/h. In congestion, when the demand is greater than the supply the FD congested branch is not impacted by the reduction of the speed limit (Fig. 4). Let us therefore assume that after the implementation of the new speed limit the traffic flow and the traffic state remain constant by six-minute period. The new speed is obtained by adding to the empirical one the change in FD: Thus the traffic density is deduced from the equation Traffic flow = Speed.Density. The variability of the speed and density is conserved, since empirical six-minute deviations from the FD are taken up again.

( ) ( ) FD FD ii new _ speed _ limit old _ speed _ limit v qv q - (2)

Predicting what should be

In 2.3.1, the new empirical FD takes into account actual compliance with the new speed limit where already implemented. Computing the ideal result of the speed limit reduction (with full compliance) requires another model, where the empirical speed (when higher than the new speed limit) is replaced by the limit (or by a function of the limit, in the case of partial compliance). Individual speeds should be used, but due to unavailability on that site they have been replaced by six-minute average speeds. Compliance is said to be partial either where the rate τ of compliant drivers or periods is <100%, or where the speed decrease is only ρ.(V old_lim -V new_lim ), (0<ρ<1). The empirical speed v, when greater than the new speed limit, is changed to v' or v": v'= (1-τ) v + τ .V new _ lim , 0%≤τ≤100% -see results in Tables 3456for τ =100% or 50%

v"= Max{V new_lim , v -ρ (V old_lim -V new_lim )} see results in Tables 3456for ρ =50% (4)

Predicted accidents related to speed and traffic density

Two models are considered here, for modelling the impacts of speed and density on accident count.

Relationships between accident risk per vehicle-kilometre and empirical speed

Since national traffic and accident data before-after a decision to reduce speed limit was not available, the speed power model is calibrated from type-2 accidents and speed data from a motorway network; this model therefore has a different meaning and final use from [START_REF] Nilsson | Traffic Safety Dimensions and the Power Model to Describe the Effect of Speed on Safety[END_REF]. The accident count is based on the predicted speeds compared with the empirical speeds before, which differ from the speed limits.

Let nb_acc(v) be the injury accident count at speed v and v b ; the "before" speed the power model is:

nb_acc(v)/nb_acc(v b )= (v/v b ) α or nb_acc(v)=βv α with β=nb_acc(v b )/(v b ) α (5)
The risk per vehicle-kilometre appears when dividing nb_acc(v) by the number of vehicle-kilometres:

Risk(v) nb_ac c(v)/nb_veh_km(v) v /nb_veh_km(v) ' v α α = = β⋅ = β ⋅ (6) with ( ) b b ' nb_ac c(v )/ v .nb_veh_km(v) α β= (7)
The (theoretical) number of accidents Nb_acc is the sum, over the speed classes, of the risks per vehiclekilometre, multiplied by the number of vehicle-kilometres of the class nb_veh_km(v); so β' is solution of:

vv v Nb_ac c nb_ac c(v) nb_veh_km(v) ' v ' nb_veh_km(v) v αα = = ⋅β⋅ =β⋅ ⋅ ∑∑ ∑ (8)
Empirical speed, flow and state (congested or not) BEFORE FD speed +deviation from FD speed BEFORE New FD speed + same deviation from FD AFTER Predicted speed (New FD+deviation) AFTER Same flow, same state (congested or not) R. Elvik (2009) gave the following values for α: for all injury accidents (1.6), for slight injury accidents (1.1), for serious injury accidents (2.6), for fatality accidents (4.1). Note that [START_REF] Elvik | Are-parameterisation of the Power Model of the relationship between the speed of traffic and the number of accidents and accident victims[END_REF], following [START_REF] Hauer | Speed and Safety[END_REF], found an even better fit with an exponential model, the number of accidents being proportional to e 0.034v .

Relationships between accident risk per vehicle-kilometre and traffic density or sensor occupancy

Occupancy and traffic density are related: assuming a constant vehicle length of 5 metres, the maximum density per lane is 200 vehicles/km, whereas the maximum occupancy is 100%, so the density is twice the occupancy (in percentage terms). (Ait-Belkacem, 2012) calibrated a power model linking the third-type accident count to occupancy (Table 1) from a French motorway network. However, when predicting the new traffic conditions on the ALLEGRO network, we remarked that while the speed limit reduction affects the proportion of high-density vehicle-kilometres (we say that a density is "high" when greater than the threshold K h = 24 vehicles per lane), it does not much affect the traffic density distribution (conditionally to be in the low/high class), which simplifies the calculation. Indeed, let us note: nb_veh_km(k) the number of vehicle-kilometres for density class k, Nb_Veh_Km(K h ) (with a capital N) the cumulated number of vehicle-kilometres for k > K h , Risk(k) the risk per vehicle-kilometre for density k, f k the conditional density (conditionally to

k > K h ): f k = nb_veh_km(k)/ Nb_Veh_Km(K h ). Nb_Acc(K h )=Σ k|k>Kh nb_acc(k)= Σ k|k>Kh nb_veh_km(k).Risk(k)= Nb_Veh_Km(K h ). Σ k|k>Kh f k .Risk(k) (9) Σ k|k>Kh f k .Risk(k) remains constant (equal to Nb_acc / Nb_Veh_Km(K h ),
provided the conditional density f k does not change. Thus, whatever the risk Risk(k), the number of high-density accidents is proportional to the number of high-density vehicle-kilometres. Two auxiliary problems remain: for the determination of K h and, when traffic data corresponding to accidents is missing, for the assignment of each accident into one of the two cases.

The predicted accident count Nb_acc is obtained through the following scheme: The speed limit will decrease from 130 km/h to 110 km/h in three motorway areas: A (A25 motorway), C (A23), and D (A27); it will decrease from 110 km/h to 90 km/h in areas E (A25), F (A1), G (weaving section A22-A27-A23), I (A22); additional sections PC1 and PC2 (A1), B (national road 41) and H (national road 356) are not considered here. 345 injury accidents, described in the French accident database BAAC (Bulletins d'Analyse des Accidents Corporels) occurred in 2003-2010 on this 23 km-long part of the Allegro Network. Traffic data from 15 stations was available only for 2009-2010, covering lane flow, speed and occupancy every six minutes. However it was not possible to link an accident to the traffic conditions before, because there were too few traffic stations

Calibration of the average speed-accident and density-accident relationships

Since accidents often result from just one person driving at an inappropriate speed, the link between the average speed of a set of drivers and accidents is more tenuous, which might jeopardize the approach. As the initial lane of the accident is not available in the BAAC, the speed-accident or density-accident relationships have been calibrated independently by lane, as if each accident were independently related to the traffic conditions on each lane. The distributions of the 6-minute average speed and occupancy have been established for the full year; in addition, for every accident the 6-minute average speed and occupancy recorded 18 minutes before the accident, at the nearest upstream or downstream traffic station, have been considered, in order to ensure that this data has not been impacted by the accidents (AIT-BELKACEM, 2012). The following results have been obtained: Density power model (Daylight, no rain, 144 multi-vehicle accidents)* 1.4 (1.0;1.9) 1.4 (0.9; 1.8) 0.9(0.4; 1.3) % of multi-vehicle accidents by density ≥ 24 vehicles/km/lane 70% of accidents concerning at least 2 vehicles (*) Unused, due to the assumption that the density distribution (restricted to the class of low densities) is constant; it is the same for high densities.

The arithmetic average over the 3 lanes of the speed exponents α is α=2.3, greater than the values 2 in [START_REF] Nilsson | Traffic Safety Dimensions and the Power Model to Describe the Effect of Speed on Safety[END_REF] or 1.6 in (Elvik, 2009). This increase is not unexpected because we discarded from the calibration set accidents involving several vehicles, which tends to sharpen the speed-accident relationship.

5. Ex-ante assessment results of the speed decrease in the ALLEGRO motorway network.

Fundamental Diagrams on the ALLEGRO network

One way to obtain the relationship between speed and density is to classify the traffic flow, then to divide each class into two sub-classes (one congested, the other not) and to assign a sub-class to each six-minute period according to its traffic flow. The FD speed is the average of speeds by sub-class. Congestion is selected when the speed is lower (or density greater) than the critical speed (or density). These critical speeds and density can be:

• derived from the calibration of an analytic curve which globally models the FD;

• or estimated from the data subset corresponding to the maximum observed flow;

• or the values which correspond to the minimum of the sum of square deviations between the empirical data and the FD resulting from these critical values. Here this method has been employed, while forcing the same critical value on the three lanes (35 vehicles/lane/kilometre).

Unlike usual traffic simulation, no global form of a speed-density relationship is required. The main constrainta decrease in speed when density increases -is not always satisfied with empirical data, and it is necessary to process monotonically the obtained average speeds by class. We have developed a heuristic which does so.

On some parts of the ALLEGRO network, the 90 km/h, 110 km/h and 130 km/h speed limits are already used, which makes it possible to build the three FDs (see Fig. 4).

Overestimation of the critical speed induces overestimation of the number of congested periods and the speed levels of both branches of the FD. This is because the average speed for a given flow class corresponds to the weighted average of the mean non-congested and congested speeds, the weights being the proportions of vehicles travelling during non-congested and congested periods. Overestimation of the critical speed implies that some non-congested periods, with a rather low speed, are wrongly classed in congestion; thus the non-congested average speed increases, because the lowest speed values are not integrated in this average; the congested average speed increases too, because of the addition of these wrongly classified speeds, which are rather high for the congested case. Of course, these two overestimations do not lead to any overestimation of the "average of the two averages", which does not vary, because of the increase in the weight of the lowest average. In the central and fast lanes, capacity is the same whatever the speed limit; in the slow lane, capacity is the same at 110 km/h and 130 km/h, but decreases at 90 km/h. Since the sections used for determining the FDs were different, the FD modification cannot be completely attributed to the change in speed limits, and part of it is probably due to local geometric or traffic characteristics which, in turn, are not independent of the speed limit decisions.

Splitting the accidents into three types and predictions

By analogy with the Marius network where detailed traffic data was available, 70% of type-3 accidents concerning at least 2 vehicles have been assigned in a low density period, and 30% in high density periods. Empirical speeds are rather low compared with the actual speed limits, but the predicted FDs are lower even by medium density; their use leads to a big decrease in speed, bringing high densities instead of medium densities and thus a high increase in high-density type-3 accidents. Conversely, when the new speed is a function of the speed limit and compliance (in the last 3 columns), the decrease in speed (and in the number of type-2 accidents) is lower and effective only at high speeds and low densities; the passage from medium to high density is more rare and the increase in high-density type-3 accidents is very low. Type-3 risk for high density is four times greater than for low density. For 8 traffic flow regimes of the fundamental diagram, [START_REF] Golob | Freeway safety as a function of the traffic flow[END_REF] gave the accident rate (property damage/injury/fatality) per million vehiclemiles):

• On the non-congested branch of the FD, the accident rate passes from 1.28 per million vehicle-miles (light flow) to 1.49 (mixed free flow), then 1.03 (heavy, variable free flow) and lastly 0.55 (approaching capacity). • On the congested branch of the FD, the accident rate passes from 1.24 (heavy flow at moderate speed) to 2.97, then 3.21, and lastly 5.99 (heavily congested flow).

Since in [START_REF] Golob | Freeway safety as a function of the traffic flow[END_REF] the units differ from ours, and injury/fatality accidents are only in 28% of his accident database (which includes property damage), multiplying these accident rates by 100*0.28/1.852 is required before any matching with our results. Having said this, what is important is the range of the accident rate according to the traffic flow regime (which is very high in both researches) rather than the numerical values. The accident count decreases from 345 to 338.4 (-1.9%); in the case of full compliance with the new speed limits it would reach 10.5 accidents (3% of the accidents). This decrease might be underestimated, due to conservative approximations:

• Use of average speeds, which decrease the link between speed and number of accidents.

• The factor (V/V b ) 2,3 has been applied here on the average speed and an "average lane" over the two years 2009-2010, and not for each 6-minute period (weighted per vehicle-kilometre)and lane. Given the convexity of the function (x) 2,3 , the application of this function on the average is greater than the average of the applications on every point; thus (V/V b ) 2,3 is overestimated (although remaining less than 1), leading to an overestimated speed and an underestimation of the decrease in accidents. • No decrease in the risk of accidents concerning several vehicles at low density, although relative speed should decrease; no decrease in type-1 accidents (due to fatigue), although it should be easier to control a slower vehicle. • The night speed exponent (between 4.1 and 6.1, see Table 1.) was discarded because it was insignificant, and because type-1 accidents were not excluded from the calibration set. • The rain exponent was discarded because it was insignificant; the predicted decrease in accidents should be higher.

Three problems have not been considered here, due to a lack of knowledge for modelling them:

• A possible change in lane assignment and thus in density in certain lanes;

• An increase in motorcycles between the lanes.

• Inside a lane, an increase in relative speed if some drivers comply and others not.

Results are highly sensitive to the speed model. This becomes more evident when the two parts of the network are analysed independently. In "interurban areas" A, C, D (where the speed limit decreases from 130 km/h to 110 km/h) and "urban areas" E, F, G, I (where the speed limit decreases from 110 km/h to 90 km/h), variations in speed (and therefore in the number of accidents) are quite different, as shown in Tables 56. Because of the non-linearity of the power model, the decrease in type-2 accidents is not "additive", i.e. equal to the sum of the decreases obtained when the power model is applied independently on the two parts (urban and interurban) of the network. A better application would have been to apply the power model by speed class.

Variations in the number of type-3 accidents on the urban and the interurban parts of the network are not additive either. Let: -N u , N i , N=N u +N i be the vehicle-kilometres travelled in the urban or interurban parts, or the whole network, n u , n i , n=n u +n i be the vehicle-kilometres passing from low to high density due to the new speed limit, -δ u , δ i , δ be the differences between the high density and the low density risks per vehicle-kilometre.

The variations in the number of type-3 accidents are respectively n u .δ u , n i .δ i , n.δ. However (n u +n i ) δ ≠ n u .δ u +n i .δ i , since δ is an average between δ i and δ u , weighted by N u and N i and not by n u and n i (Tables 56).

With the FD model, most decreases in accident count are obtained in the interurban part (Table 5), where speed decreases from 102.4 to 89.7 km/h; in the urban part, when the speed limit passes from 110 km/h to 90 km/h (Table 6) the speed and accident decreases are very low. It is quite the opposite with the full compliance model (speed decreases from 89.1 to 83.1 km/h on the urban part)! There are two possible explanations:

• An underestimation of the speed in the FD established for a 110 km/h speed limit would imply an overestimation of the decrease in accidents (Table 5) when passing from 130 km/h to 110 km/h, and an underestimation of the decrease in accidents when passing from 110 km/h to 90 km/h (Table 6). Indeed, the speed in the central lane, in the 110km/h FD (Fig. 2b), is always less than 100 km/h, which means that following this FD leads to a lower speed than compliance in the 110 km/h speed limit. In the fast lane (Fig. 2c) the FD speed is always less than 110km/h. The FD speed was not really underestimated, but using for the interurban part of the network an FD coming from the urban part of the network, where the network configuration (with more ramps) and driver behaviour are different, is questionable. Furthermore, is it correct to say that, regarding the urban part of the network, a speed lower than 100 km/h or 110 km/h has incited traffic operators to decrease the speed limit? • Low compliance with the 90 km/h speed limit in France, which limits the decrease in accidents (Table 6).

The fourfold validation of the approach.

The approach is based on two traffic models and on two models linking traffic conditions to accidents. Some assumptions intervene either in the models themselves, or in their calibration, or in their use. Validation will therefore concern the models, their calibration, their use and the results.

Validation of the models

Two traffic models are used to predict the traffic conditions, and two models linking accidents to traffic:

• The first traffic model (prediction of what will be), based on FDs, does not manage a change in the queues from one period to the next. Changes in lane assignment are also ignored. If, when considering the ex-post traffic conditions, these assumptions are not validated, the use of a simulation model is required. • The second one (prediction of what should be) replaces the actual speed (when higher than the new speed limit) by a function of this limit, based on the compliance rate. This must be applied on individual (and not average) speeds. The impact of aggregation will be computed in a network where individual speeds are available (like the Marius network) in order to validate or invalidate this model. • The form of the relationship between speed and accident (the power model) is commonly accepted by the scientific community, although confidence in the exponential model is progressing. • The accident/ density relationship is simplified, thanks to the assumption that there will be no change in the density distribution, conditional to density < 24 vehicles/lane (respectively density > 24 vehicles/lane). If this is not validated, the calibration and use of another model is required.

Validation of the calibrations

New calibrations of the relationships between traffic conditions and accidents (for instance with the ALLEGRO traffic and accident database, on the before period) will validate/invalidate the relationships used: 1. For the first accident type, validating "no relationship with traffic condition" requires, first of all, the calibration of a model such as a logistic model linking the risk of accident to speed and/or density. In the second step, the part of the accident variance explained by this model has to be analysed. The logistic model can be withdrawn (e.g. the assumption that there is no relationship for type-1 accidents is validated) if this part of the variance is low (10% or 20% of the variance, say) or if the parameters α' obtained are close to zero. 2. For the second accident type, a new calibration α' of the exponent of the speed power model is necessary. 3. For the third accident type, a new calibration of threshold separating low and high densities is necessary.

Let traffic operators define, by type of accident, a maximum bias (positive or negative) on the predicted accident count. It is possible to convert these biases into positive and negative deviations round parameters α of the relationships linking accident counts with the traffic conditions, and thus to define an "acceptable" interval round α; the relationships used will be validated if α' falls within this interval.

Verification/Validation of the use of the models

• It is not recommended, in a speed limit reduction framework, to use traffic conditions/accident relationships that are calibrated without any specific speed limit decrease. But road safety measures cannot wait. • Using the 110 km/h FD which has been calibrated on an urban network for an interurban network is questionable; this urban FD should be compared with the ex-post interurban one.

Comparisons with the final results

• An ex-post analysis of the traffic conditions will improve the understanding of the speed limit reduction, and will validate/invalidate the predicted speed and occupancy distributions. The distance between the predicted and the new empirical distributions can be measured by the Kullback-Leibler divergence [START_REF] Kullback | Letter to the Editor: The Kullback-Leibler distance[END_REF], which is the basis of the Akaike Information Criterion, or tested using the Kolmogorov-Smirnov statistics. A discrepancy between the two distributions corresponding to a shift of 2 or 3 km/h in speed seems acceptable. • Unfortunately, due to the low number of accidents, it is not expected that the comparison of the ex-post accident count with the ex-ante prediction will provide statistically significant elements for validation.

Conclusions.

A speed limit decrease generally induces a decrease in the road injury and fatality accident count. This paper presents a method for the ex-ante assessment of such a decrease, and its application on a French motorway network. Some results on the calibration and the use of relationships linking traffic density or speed to accident counts are given. The envisaged decrease in the speed limit by 20 km/h should lead to a 1.9% decrease in the accident count with the current compliance, and to a 3% decrease with full compliance. These figures are a trade-off between a null trend for certain accidents (due to driver or vehicle failure), a decrease in accidents linked to speed, and an increase in accidents linked to a high density. However, this prediction is very sensitive to the traffic and safety models, to their calibration and to the conditions of their use.

Breaking down the prediction into the urban and interurban parts of the network shows that most decreases are located on the interurban part of the motorway network (4.7 accidents) where the speed limit decreases from 130 km/h to 110 km/h and the average speed decreases by 13.7 km/h. On the urban part, where the speed limit decreases from 110 km/h to 90 km/h, the predicted speed and accident decreases are much lower (0.7 accidents).

Fig. 1 .

 1 Fig. 1. Traffic conditions analysis and prediction

Fig. 2 .

 2 Fig. 2. Predicting accidents by type 3. The part of the ALLEGRO motorway network concerned by speed limit reduction.

Fig. 3 .

 3 Fig.3. Speed limit reduction scheme near Lille

Fig. 4 .

 4 Fig.4. Fundamental Diagram in the ALLEGRO motorway network according to the speed limit for daylight, no rain periods.

  Since traffic data covers two years and accident data covers eight years, risks per vehicle-kilometre are computed with the assumption that the traffic in 2003-2008 was the same as in 2009-2010.

Table 1 .

 1 Calibration of the exponents of the power model -Marius network.

	Exponent α and (5% confidence interval) of the power model	Slow lane	Central lane	Fast lane
	Speed (Daylight, no rain, based on 45 single-vehicle accidents)	0.2(-1.7;2.7)	0.7 (-1.5;2.9)	6.0(2.4;9.7)
	(Night, no rain, based on 17 single-vehicle accidents	4.1(0.6; 8.8)	6.1 (-5.1;17.3)	6.1(2.1;10.0)

Table 2 .

 2 Observed accidents and predicted accidents by accident type according to speed limit compliance assumptions

	Accidents, areas A to I	Observed before	Predicted by FDs	Compliance τ=100%	Compliance rate τ=50%	ρ=50% (-10km/h)
	2003-2010 Type 1	46.0	46.0	46.0	46.0	46.0
	2003-2010 Type 2 (due to speed)	99.0	85.9	88.5	93.5	90.8
	2003-2010 Type 3 | density<24) (*)	140.0	137.8	139.8	139.9	139.8
	2003-2010 Type-3 | density>24 (**)	60.0	68.7	60.9	60.3	60.7
	(*) 70% of type-3 accidents (**) 30% of type-3 accidents

Table 3 .

 3 Observed and predicted speed and density, according to speed limit compliance assumptions

	Global predictions speed, density, accidents	Observed FDs	τ=100%	rate τ=50%	ρ=50%
	Speed (km/h)	93.60	88.04	89.33	91.47	90.30
	Decrease for the 99 Type-2 accidents	0.00	-13.10	-10.50	-5.50	-8.20
	2009-2010 millions vehicle-km |density<24	652.30	641.90	651.20	651.90	651.40
	2009-2010 millions vehicle-km |density>24	71.50	82.00	72.60	71.90	72.40
	% vehicle-km| density>24	9.88%	11.32%	10.03%	9.94%	10.00%
	% increase in high density vehicle-km	0.00%	14.60%	1.50%	0.60%	1.20%
	2003-2010 Accidents (total)	345.00	338.40	335.20	339.70	337.30

Table 4 .

 4 Empirical and predicted risks according to speed limit compliance assumptions

	Risks per 100 million vehicle-kilometres	Observed before	Predicted by FDs	Compliance τ=100%	Compliance rate τ=50%	ρ=50% (-10km/h)
	Type-1 risk	1.59	1.59	1.59	1.59	1.59
	Type-2 risk	3.42	2.97	3.07	3.24	3.15
	Type-3 risk	6.91	7.13	6.93	6.91	6.93
	Type-3 risk | density<24/lane/km)	5.37	5.37	5.37	5.37	5.37
	Type-3 risk | density≥24 vehicle /lane/km)	20.98	20.95	20.97	20.97	20.96
	Total risk	11.92	11.69	11.59	11.75	11.66

Table 5 .

 5 Observed accidents and predicted accidents, interurban part

	130 km/h  110 km/h	Observed before	Predicted by FDs	Compliance τ=100%	Compliance rate τ=50%	ρ=50% (-10km/h)
	Average speed	102.4	89.7	101.4	101.9	101.4
	Type-2 accidents	44.0	32.5	42.9	43.4	43.0
	Millions vehicle-km, density <24	235.3	228.4	235.3	235.3	235.3
	Millions vehicle-km, density >24	11.8	18.7	11.8	11.8	11.8
	Type-3 accidents, low density	30.8	29.9	30.8	30.8	30.8
	Type-3accidents, high density	13.2	20.9	13.2	13.2	13.2
	Type-2+ Type-3 accident decrease		-4.7	-1.1	-0.6	-1.0
	Table 6. Observed accidents and predicted accidents, urban part
	110 km/h  90 km/h	Observed before	Predicted by FDs	Compliance τ=100%	Compliance rate τ=50%	ρ=50% (-10km/h)
	Average speed	89.1	87.2	83.1	86.1	84.5
	Type-2 accidents	55.0	52.4	46.1	50.4	48.2
	Millions vehicle-km, density <24	417.0	413.5	416.0	416.6	416.2
	Millions vehicle-km, density >24	59.7	63.3	60.8	60.2	60.6
	Type-3, accidents low density	109.2	108.3	108.9	109.1	109.0
	Type-3, accidents high density	46.8	49.6	47.6	47.1	47.5
	Type-2+ Type-3 accident decrease		-0.7	-8.4	-4.4	-6.3

Acknowledgements.

We thank the French traffic operators DIR-Nord and CETE-Méditerranée, who kindly gave us access to their traffic data. Thanks also to our students Y. Saas and I. Ait-Belkacem, who helped us to do the work.


Loïc Peter

Diana Mateus

Pierre Chatelain

Noemi Schworm

Stefan Stangl

Gabriele Multhoff

Nassir Navab

Leveraging Random Forests for Interactive Exploration of Large Histological Images

The large size of histological images combined with their very challenging appearance are two main difficulties which considerably complicate their analysis. In this paper, we introduce an interactive strategy leveraging the output of a supervised random forest classifier to guide a user through such large visual data. Starting from a forest-based pixelwise estimate, subregions of the images at hand are automatically ranked and sequentially displayed according to their expected interest. After each region suggestion, the user selects among several options a rough estimate of the true amount of foreground pixels in this region. From these one-click inputs, the region scoring function is updated in real time using an online gradient descent procedure, which corrects on-the-fly the shortcomings of the initial model and adapts future suggestions accordingly. Experimental validation is conducted for extramedullary hematopoesis localization and demonstrates the practical feasibility of the procedure as well as the benefit of the online adaptation strategy.

Introduction

Analyzing histological images is usually an extremely challenging task. Due to the complex appearance of objects of interest, the accuracy of fully-automatic techniques is generally insufficient for clinical use, e.g. in the case of mitosis detection [START_REF] Veta | Breast cancer histopathology image analysis: A review[END_REF] for which a variety of automatic approaches has been recently compared quantitatively [2]. In fact, the accurate identification of patterns within such images is often only achievable by well-trained human experts and remains prone to inter-experts disagreements in some cases [START_REF] Crowley | Development of visual diagnostic expertise in pathology -an information-processing study[END_REF]. Moreover, histological data are very tedious to process for a human because of their large dimension, which commonly reaches tens of thousands of pixels along each direction. When aiming at finding rare objects within images of this size, a manual search for these instances requires a painstaking exploration of the whole content, and a huge amount of time is spent scrolling through uninteresting background areas.

To overcome this, we propose to leverage fully-automatic pixelwise classification techniques to recover candidate areas of interest (i.e. where positive instances are found), leaving the task of interpreting the visual content within these regions to an expert user. To this end, assuming that an available automatic method provides a label confidence at each pixel, an intuitive strategy would be to rank subregions of the images in decreasing order of confidence and display them accordingly one after the other. Having thereby patterns of interest shown early in the process, an exhaustive search would not be necessary anymore. However, such a scenario relies entirely on the accuracy of the automatic detector and can fail if the latter suffers from difficulties on the images at hand. In the case of detectors based on supervised learning procedures, difficulties may arise due to the common mismatches between the training and test data, caused by e.g. differences in noise or illumination. In this context, it would be beneficial to allow the user to interact after each suggestion (in particular after wrong ones) so that the procedure can be adapted to the characteristics of the test data and eventually avoid mistakes for the upcoming region suggestions. Online transfer learning (OTL) frameworks [START_REF] Zhao | OTL: A framework of online transfer learning[END_REF][START_REF] Tommasi | Leveraging over prior knowledge for online learning of visual categories[END_REF] address this problem by combining (i) a kernel-based classifier pre-trained on some source data, and (ii) an online classifier continuously trained on arriving instances from the target data, given by the user input in our case. On the closely related problem of domain adaptation, an unsupervised algorithm based on Gaussian process regression was introduced to adapt the decision boundary of a face detector to a test image [START_REF] Jain | Online domain adaptation of a pre-trained cascade of classifiers[END_REF]. In the case of medical data, random forest classifiers [START_REF] Breiman | Random forests[END_REF][START_REF] Criminisi | Decision Forests for Computer Vision and Medical Image Analysis[END_REF] demonstrated high accuracy and tractability and are especially used for their ability to handle high amounts of training data and descriptors. Although an online version of random forests exists in the literature [START_REF] Saffari | On-line random forests[END_REF], replacing the kernel classifiers of the OTL approach by forests is not straightforward: only few samples are available to train the target classifier, and a forest-based OTL scheme would handle the source classifier as a whole, without exploiting the diversity and possibly heterogeneous relevance of its local models on the target data.

In this work, we introduce a regression-based method able to leverage the knowledge carried by a pre-trained random forest classifier to guide the user through the data by suggesting candidate regions of interest. Our approach exploits the multiple partitionings of the feature space defined by the trained forest to perform local updates of the model from one-click inputs provided by the user.

Thereby, characteristics of the test data can be captured and future suggestions are adapted accordingly. These updates are performed in real time between two iterations and do not require the original training data, which makes our method particularly tractable for clinical use. Experimental evaluation is conducted on a high-resolution dataset for hematopoesis identification within mouse liver slices, as well as on synthetic deformations of it, and demonstrates the general feasibility of our approach and the benefit of the adaption scheme.

Methods

We first summarize the general scenario of our approach. Let us consider a given set of images composed of pixels. Each pixel p has a true label y(p) ∈ {0, 1} and we aim at detecting the positive instances within our images, i.e. the pixels p such that y(p) = 1. We make the assumption that a random forest classifier has been trained beforehand on some available labeled data and hence provides a probability P RF (y(p) = 1) ∈ [0, 1] for each pixel p. Here, we implicitly assume that this forest has been originally designed for a binary classification task, but any multi-class forest could also be used by grouping as positive the labels considered as interesting for the given application. For tractability, the original training data is not made available during the whole process and is thus only encoded through this forest. In our scenario, a user runs this classifier which is based on prior data on a new set of images that are correlated but potentially slightly different, e.g. because of different imaging conditions. From the resulting pixelwise output, an area of interest, expected to contain positive labels, is then displayed to the user, who provides in return an estimate of the proportion of positive pixels actually observed. In doing so, the bias of the initial model can be progressively assessed and compensated to ultimately increase the relevance of the upcoming region suggestions.

Under these conditions, we can formalize our scenario mathematically. The test images are partitioned into a set of rectangular regions R = {R 1 , . . . , R n } whose sizes are tractable for a human user. At iteration k, i.e. after the user saw k -1 regions, we proceed as follows:

1. The region R ∈ R, maximizing a scoring function φ k : R → R stating the expected relevance of each region, is displayed to the user. 2. The user interacts to provide information about the content observed in R. 3. Using this feedback, φ k is updated, resulting in a new scoring function φ k+1 . 4. R is removed from the pool of candidate regions R.

The initial scoring function φ 1 is based on the initial random forest classifier only. We summarize the random forest model in Sec. 2.1. In Sec. 2.2, we expose our choice of scoring functions φ k and how they can be parametrized. Finally, the update strategy of φ k from the user inputs is described in Sec. 2.3.

Random Forest Model

Our approach uses a trained random forest classifier as starting point, which provides for each pixel p a probability of being a positive instance. A random forest is a collection of T decorrelated binary decision trees. Each tree (indexed by t ∈ {1, . . . , T }) is a hierarchical collection of decision rules based on visual features, leading to a partition of the domain into l(t) leaves L t,1 , . . . , L t,l(t) . Each leaf L t,i contains a probability model π t,i ∈ [0, 1]. For each tree t, every newly observed pixel p reaches exactly one leaf, and we will denote π t (p) the corresponding probabilistic model. Averaging over trees, the final pixelwise probability provided by the random forest for each pixel p is

P RF (y(p) = 1) = 1 T T t=1 π t (p) . (1)
We follow the standard procedure [START_REF] Breiman | Random forests[END_REF][START_REF] Criminisi | Decision Forests for Computer Vision and Medical Image Analysis[END_REF] to train such a classifier, i.e. to infer the decision rules and the leaf models. Starting from a root node, the decision rules are recursively chosen within a set of randomly drawn splitting functions to maximize an information gain criterion. The proportion of positive training samples contained in a leaf defines its probabilistic model.

Region Scoring Function

In this subsection, we describe our region scoring model. Given a pixelwise estimate P k (y(p) = 1) at iteration k, we define the relevance of a region R ∈ R as the expectation of the proportion of positive pixels in R. Modeling the pixel labels as independent Bernoulli distributions of success probability P k (y(p) = 1), this expectation and hence the scoring function φ k can be written as

φ k (R) = 1 |R| p∈R P k (y(p) = 1) , (2)
where |R| denotes the number of pixels in the region R. In this work, we would like to define the pixelwise probabilistic model P k as a combination of (i) the prior knowledge acquired on training data that is encoded by the initial random forest classifier, and (ii) the k -1 inputs provided by the user based on the test data. Therefore, we propose to obtain P k (y(p) = 1) from the forest estimate of Eq. 1 by adding to each prior leaf model π t,i a signed offset ǫ k t,i ∈ R whose value will be progressively adapted from the user inputs. Let us denote ǫ k the vector (ǫ k t,i ) t,i and rewrite the scoring function φ k (R) as φ(R|ǫ k ) so that this parametrization clearly appears. Introducing this additive model in Eq. 2 and after some algebraic manipulations, the scoring function can be rewritten as

φ(R|ǫ k ) = φ RF (R) + φ new (R|ǫ k ) , (3)
where φ RF (R) is obtained by using the forest probabilistic output (Eq. 1) as pixelwise model in Eq. 2, and

φ new (R|ǫ k ) = 1 |R| 1 T T t=1 l(t) i=1 ǫ k t,i h t,i (R) = r, ǫ k . (4)
In Eq. 4, h t,i (R) denotes the number of pixels p in R which reached the leaf L t,i , and r = ( 1 |R| 1 T h t,i (R)) t,i . Before the user starts providing inputs, the ǫ 1 t,i are set to 0 such that φ 1 = φ RF . As the user gives information after each suggestion, the acquired knowledge is encoded through the leaf-dependent ǫ k t,i that allow the scoring function φ k to vary around the initial model φ RF . The next subsection describes how these ǫ k t,i are inferred from the collected user feedbacks.

Online Model Update

After each suggestion of a region R j , a feedback P user (R j ) ∈ [0, 1] stating the true proportion of positive pixels in R j is requested from the user. We assume for the sake of formalism that the exact proportion is given, but we will demonstrate in our experiments that providing a discretized approximation of this input does not decrease the performance (Tab. 1). This reduces drastically the necessary amount of interaction, since giving a discrete estimate of this proportion requires only one click from the user, whereas computing the exact value would involve a pixelwise region labeling. Note that this feedback does not provide any information at pixel level and hence excludes a direct update of the leaf statistics [START_REF] Saffari | On-line random forests[END_REF]. Assuming that k regions R 1 , . . . , R k have been suggested and their respective user inputs P user (R 1 ), . . . , P user (R k ) collected (k ≥ 1), we would like to find a new scoring function φ k+1 , or equivalently a coefficient vector ǫ k+1 , which matches with these k user inputs. Considering a squared loss function, we solve the following minimization problem

ǫ k+1 = argmin ǫ k j=1 (φ(R j |ǫ) -P user (R j )) 2 + λ ǫ 2 , (5)
which can be rewritten from Eq. 3 and Eq. 4 as a linear ridge regression problem

ǫ k+1 = argmin ǫ k j=1 ( r j , ǫ -δ j ) 2 + λ ǫ 2 , ( 6
)
where δ j = P user (R j ) -φ RF (R j ). The regularization term aims at keeping the deviation from the initial random forest classifier small and is driven by a hyperparameter λ. The minimization problem defined by Eq. 6 is equivalent to a linear ridge regression problem with k training instances of feature vectors r 1 , . . . , r k and labels δ 1 , . . . , δ k , i.e. a regression problem at the region level. As a consequence, Eq. 6 could be naturally extended to a kernel ridge regression problem if one desires to remove the linearity assumption. The solution ǫ k+1 of Eq. 6 can be found in closed form [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF]. However, computing this exact value at each iteration involves the computationally expensive inversion of a large matrix, which could force a user to wait between two suggestions of regions. Instead, we adopt an online gradient descent scheme [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF] leading to the following incremental update rule:

ǫ k+1 = ǫ k - 1 λ ( r k , ǫ k -δ k ) r k . (7)
Thereby, the k th user input can be incorporated in real time in our model.

Experiments

Dataset: Extramedullary hematopoesis, i.e. hematopoesis located outside of the bone marrow, is a rare event in human adults. Yet, in some animals like mice, its presence in the liver can be a reaction to an immense inflammation or to a tumor stimulating extensively the immune system of the organism [START_REF] Tao | Imagable 4t1 model for the study of late stage breast cancer[END_REF]. Localizing these relatively small lesions as well as estimating their number can be useful to compare the influence of different therapeutic strategies. 34 high-resolution images entirely labeled at pixel level were extracted from 8 mouse liver slices of resolution 0.498 µm/pixel. This restriction to subimages was necessary to obtain accurate labels from a pathologist, since a slice is typically too big (34200×21052 pixels) to be labeled entirely. Since images taken from a same slice share common visual properties, we kept them together for training or testing. Images were downsized by a factor 3 to speed up the training and testing procedures.

Experimental Setup: We performed our experiments at two levels of supervision: (i) a standard leave-one-out strategy (strong supervision), where the initial random forest is trained on 7 slices (i.e. approximately 30 images) and tested on the 8 th one, and (ii) the inverse setting where one trains only on one slice and tests on the 7 others (weak supervision). While the second scenario is plausible in practical situations due to the difficulty of obtaining pixelwise labels, it also suggests another application of our approach as a way to facilitate collecting hard-to-find positive samples to enrich an existing weak training set. Training samples were collected on a regular grid of step 30, resulting in approximately 500 k instances in the leave-one-out setup. Each pixel was described by its visual content at offset locations. We used the CIELab color space that we enriched with the output of a bank of filters [START_REF] Winn | Object categorization by learned universal visual dictionary[END_REF]. 40 trees were grown without depth limit but with at least 10 training samples per leaf. The images were partitioned into regions of size 123 × 123 pixels, which fits to the size of usual hematopoesis patterns. If a hematopoesis pattern is not fully included in a displayed region, the field of view can be extended with a simple click revealing the next adjacent region in the indicated direction. We investigated ways to handle pools of regions of heterogeneous sizes, but this raised issues in terms of scoring function (finding a reliable rule that does not always favor the smallest or the largest regions is not straightforward) and evaluation. The retained solution requires almost no effort from the user and eventually seemed to be the most practical one.

Evaluation: Our approach aims at performing an interactive online adaptation to new imaging conditions, and more generally to mismatches between training and test data. The aforementioned weakly supervised scenario favors naturally such a bias, since training on only one slice is unlikely to generalize well. To go further into this direction, we propose, as additional experiments, to keep the initial training sets and apply three types of synthetic modifications on the test images to simulate classical difficulties occurring in histological imaging: (i) global illumination change, (ii) local illumination change, and (iii) local blurring. These changes were randomly generated on each image. To evaluate our method, we studied the average proportion of positive pixels recovered over time (Fig. 2) and measured the area under these curves (Tab. 1). The results demonstrate that our update framework allows a faster retrieval of interesting areas. As mentioned in Sec. 2.3, we also show that discretizing the user input into 7 predefined bins (i.e. enforcing P user (R j ) ∈ {0, 0.1, 0.3, . . . , 0.9, 1}) does not reduce the performance and thus lightens the user interaction. Finally, in the leave-one-out setting on real data, we assessed the influence of the parameter λ (Eq. 6). From Eq. 7, low (resp. high) values of λ are expected to lead to numerical instabilities (resp. negligible updates). These asymptotic behaviors were experimentally observed for λ ≤ 0.05 and λ ≥ 100, respectively. Between these two extremes, a wide range of values of λ led to improvements, and we chose λ = 1 for simplicity.

Conclusion

We proposed an interactive method able to suggest sequentially regions of interest within large data from the output of a pre-trained random forest classifier at pixel level. After each suggestion, one-click user inputs are collected to capture the specificities of the images at hand and update the region scoring function accordingly and in real time via an online gradient descent scheme. We evaluated our approach in the context of hematopoesis localization, with and without synthetic deformations simulating common sources of variability in histological imaging. As future work, we would like to extend our method to mitosis detection tasks [2] and generalize it to 2D video sequences and 3D volumes.

Fig. 1 .

 1 Fig. 1. Summary of our approach

Fig. 2 .

 2 Fig. 2. Results The curves display the average proportion of recovered positive pixels as a function of the proportion of the slice shown to the user. These experiments are driven under different synthetic conditions, examples of which are shown on the left. From top to bottom: real image, global illumination change, local illumination change, local blurring. The discretized version of the user input was used in these experiments.

Table 1 .

 1 Area under the curves in our different settings.

			Weak supervision
		Real data	Global change	Local change	Local blurring
		Min Median Max Min Median Max Min Median Max Min Median Max
	No update	82.5 90.9 95.9 63.8 77.0 95.7 56.2 76.4 81.8 78.7 83.7 93.8
	Exact input	87.6 92.2 96.1 74.4 87.8 95.7 62.9 82.8 91.6 78.0 86.1 94.1
	Discretized input 88.5 92.3 96.1 74.6 87.7 95.9 64.0 83.9 91.6 77.9 86.3 94.1
			Strong supervision
		Real data	Global change	Local change	Local blurring
		Min Median Max Min Median Max Min Median Max Min Median Max
	No update	94.1 96.6 97.7 71.1 93.1 97.1 78.5 85.4 92.7 91.2 95.4 96.0
	Exact input	94.7 96.7 98.1 89.1 93.6 97.2 83.2 91.9 95.0 92.4 95.8 97.1
	Discretized input 94.8 96.7 98.3 88.9 93.7 97.4 82.8 91.9 95.0 92.4 95.6 97.0

Acknowledgments This work was partially supported by the DFG-funded Collaborative Research Centre 824: "Imaging for Selection, Monitoring and Individualization of Cancer Therapies".


Didier Devaurs

Thierry Siméon

Juan Cortés

A Multi-Tree Extension of the Transition-based RRT: Application to Ordering-and-Pathfinding Problems in Continuous Cost Spaces

The Transition-based RRT (T-RRT) is a variant of RRT developed for path planning on a continuous cost space, i.e. a configuration space featuring a continuous cost function. It has been used to solve complex, high-dimensional problems in robotics and structural biology. In this paper, we propose a multiple-tree variant of T-RRT, named Multi-T-RRT. It is especially useful to solve ordering-and-pathfinding problems, i.e. to compute a path going through several unordered waypoints. Using the Multi-T-RRT, such problems can be solved from a purely geometrical perspective, without having to use a symbolic task planner. We evaluate the Multi-T-RRT on several path planning problems and compare it to other path planners. Finally, we apply the Multi-T-RRT to a concrete industrial inspection problem involving an aerial robot.

I. INTRODUCTION

Sampling-based path planning has traditionally aimed at finding collision-free paths to solve complex planning problems in high-dimensional spaces [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF], [START_REF] Lavalle | Planning Algorithms[END_REF]. However, beyond feasible solutions, in many applications it is important to compute high-quality paths with respect to a given cost criterion. When a cost function is defined on the configuration space, we call the latter a cost space.

Several approaches to sampling-based cost-space path planning have been proposed based on the Rapidly-exploring Random Tree (RRT) algorithm [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF], such as RRT* [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF] or the Transition-based RRT (T-RRT) [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF]. T-RRT combines the exploratory strength of RRT with a transition test favoring lowcost regions. It has been successfully applied to diverse robot path-planning problems [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF]- [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF] (some involving humanrobot interactions [START_REF] Mainprice | Planning human-aware motions using a sampling-based costmap planner[END_REF]) and structural biology problems [START_REF] Iehl | Costmap planning in high dimensional configuration spaces[END_REF], [START_REF] Jaillet | Randomized tree construction algorithm to explore energy landscapes[END_REF]. Contrary to RRT*, T-RRT does not offer asymptoticoptimality guarantees, but, in high-dimensional spaces, it may converge faster than RRT* [START_REF] Iehl | Costmap planning in high dimensional configuration spaces[END_REF], [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF].

In this paper, we propose a multi-tree variant of T-RRT, named Multi-T-RRT. Since there exist numerous multi-tree path planners that involve RRT [START_REF] Esposito | Adaptive RRTs for validating hybrid robotic control systems[END_REF]- [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF], we have evaluated existing techniques and selected the most effective ones. The Multi-T-RRT is particularly useful when looking for a path going through a given set of unordered waypoints. Such ordering-and-pathfinding problems involve two aspects: a low-level path planning problem that consists of connecting pairs of waypoints, and a high-level ordering problem that consists of finding an efficient way to visit all the waypoints (which is a simple kind of task planning problem).

All authors are with CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France and Univ de Toulouse, LAAS, F-31400 Toulouse, France (e-mails: devaurs@laas.fr, nic@laas.fr, jcortes@laas.fr) This work has been partially supported by the European Community under Contract ICT 287617 "ARCAS". Hybrid approaches to solve task-and-path planning problems are often based on decoupling the two aspects: a symbolic task planner computes a high-level plan (possibly based on geometrical data) that is refined by a path planner computing precise low-level paths [23]- [START_REF] Kiesel | Integrating vehicle routing and motion planning[END_REF]. In some cases, when tasks are simple enough, the overall problem possesses a purely geometrical formulation, and no symbolic task planner is needed [START_REF] Siméon | Manipulation planning with probabilistic roadmaps[END_REF]. In this work, we also follow a purely geometrical approach: the geometric path planner (i.e. Multi-T-RRT) yields high-quality high-level solutions based on the costs of the low-level paths it computes between waypoints. To achieve that, we have enhanced the Multi-T-RRT with a useful-cycle addition mechanism enabling it to continually improve the solution path in an anytime fashion (which is illustrated by the accompanying video).

After a brief review of T-RRT (Section II), we present the Multi-T-RRT, based on the strategies we have selected to expand and connect trees (Section III). Then, we report some evaluation results and compare the Multi-T-RRT to planners involving the Bidirectional T-RRT [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF] (Section IV). Finally, we apply the Multi-T-RRT to ordering-and-pathfinding problems, including a concrete industrial inspection problem involving an aerial robot (Fig. 1, Section V).

1 if cj ≤ ci then return True 2 if exp(-(cj -ci) / T ) > 0.5 then 3 T ← T / 2 (c j -c i ) / costRange(T ) ; return True 4 else 5 T ← T • 2 T rate ; return False

II. TRANSITION-BASED RRT (T-RRT)

Starting from an initial configuration q init , RRT iteratively builds a tree T on the configuration space C [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF]. At each iteration, a configuration q rand is randomly sampled in C, and an extension toward q rand is attempted, starting from its nearest neighbor, q near , in T . If the extension succeeds, a new configuration q new is added to T , and connected by an edge to q near . The criteria on when to stop the exploration can be reaching a given target configuration q goal , a given number of nodes in the tree, a given number of iterations, or a given running time.

T-RRT (shown in Algorithm 1) is a variant of RRT used to explore cost spaces [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF], [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF]. It extends RRT by integrating a transition test favoring the exploration of low-cost areas of the space. The transitionTest presented in Algorithm 2 is used to accept or reject the move from q near to q new based on their respective costs. A downhill move (c j ≤ c i ) is always accepted. An uphill move is accepted or rejected based on the probability exp(-(c j -c i ) / T ) that decreases exponentially with the cost variation c j -c i . In that case, the level of selectivity of the transition test is controlled by the temperature T , which is an adaptive parameter of the algorithm. Low temperatures limit the expansion to gentle slopes, and high temperatures enable it to climb steep slopes. After each accepted uphill move, T is decreased to avoid over-exploring high-cost regions: it is divided by 2 (cj -ci) / costRange(T ) , where costRange(T ) is the cost difference between the highest-cost and the lowest-cost configurations in the tree. After each rejected uphill move, T is increased to facilitate the exploration and to avoid being trapped in a local minimum: it is multiplied by 2 Trate , where T rate ∈ (0, 1] is the temperature increase rate. Following [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF], we set T rate to 0.1 and initialize T to 10 -6 .

III. MULTI-TREE VARIANTS OF T-RRT

This section introduces our multi-tree variant of T-RRT, named Multi-T-RRT. To develop it, we have surveyed several techniques proposed in similar work on multi-tree approaches to sampling-based path planning [START_REF] Esposito | Adaptive RRTs for validating hybrid robotic control systems[END_REF]- [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. Some approaches aim at solving single-query problems, the way RRT usually works, but involve the construction of several RRTs to reach a solution [START_REF] Esposito | Adaptive RRTs for validating hybrid robotic control systems[END_REF]- [START_REF] Flavigné | Improving motion planning in weakly connected configuration spaces[END_REF]. Others are multiplequery approaches similar to the Probabilistic Road-Map (PRM), where RRT is used as a local planner [START_REF] Morales | Improving the connectivity of PRM roadmaps[END_REF]- [START_REF] Plaku | Samplingbased roadmap of trees for parallel motion planning[END_REF]. Others focus on dynamic environments and build several RRTs at different points in time [START_REF] Gayle | Lazy Reconfiguration Forest (LRF) -an approach for motion planning with multiple tasks in dynamic environments[END_REF], [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. The version of the Multi-T-RRT we present here is a single-query planner building several T-RRTs to find a path. We do not deal with multiple queries or dynamic environments.

Growing multiple trees on the configuration space can be done in various ways. The aim can be to have several RRTs rooted in different regions of the space to ensure a broader exploration [START_REF] Clifton | Evaluating performance of multiple RRTs[END_REF]- [START_REF] Flavigné | Improving motion planning in weakly connected configuration spaces[END_REF]. In this context, trees are initialized and grown rather independently of one another. Other approaches aim at maintaining a road-map of RRTs over the space [START_REF] Li | An incremental learning approach to motion planning with roadmap management[END_REF]- [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. In this case, trees can be created or modified as a result of merging, splitting or pruning operations. Other approaches make use of sub-trees produced by previous queries [START_REF] Li | An incremental learning approach to motion planning with roadmap management[END_REF], [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. Others build RRTs in different subspaces, independently of each other [START_REF] Esposito | Adaptive RRTs for validating hybrid robotic control systems[END_REF], [START_REF] Belta | Computational techniques for analysis of genetic network dynamics[END_REF]. Finally, RRTs can be reduced to local connections between components of a large road-map [START_REF] Morales | Improving the connectivity of PRM roadmaps[END_REF]. In this work, we focus on growing several T-RRTs rooted at given waypoints.

When building several trees, controlling their number and the timing of the connection attempts are difficult issues [START_REF] Lavalle | Planning Algorithms[END_REF]. First, the number of trees can be unbounded [START_REF] Li | An incremental learning approach to motion planning with roadmap management[END_REF]- [START_REF] Gayle | Lazy Reconfiguration Forest (LRF) -an approach for motion planning with multiple tasks in dynamic environments[END_REF]. It can also be subjected to a pre-defined bound [START_REF] Ettlin | Randomised rough-terrain robot motion planning[END_REF]- [START_REF] Wang | An adaptive roadmap guided Multi-RRTs strategy for single query path planning[END_REF] or implicitly limited at runtime [START_REF] Flavigné | Improving motion planning in weakly connected configuration spaces[END_REF], [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. Second, the tree roots can be sampled a priori [START_REF] Ettlin | Randomised rough-terrain robot motion planning[END_REF] or at runtime [START_REF] Plaku | Samplingbased roadmap of trees for parallel motion planning[END_REF]- [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. They can also be strategic states discovered at runtime, such as configurations in narrow passages [START_REF] Strandberg | Augmenting RRT-planners with local trees[END_REF]- [START_REF] Flavigné | Improving motion planning in weakly connected configuration spaces[END_REF]. We focus here on the case where the number of trees is fixed and equal to the number of waypoints.

A. Multi-T-RRT

The pseudo-code of the Multi-T-RRT is presented in Algorithm 3. Instead of building a single tree, we build n trees rooted at n given waypoints q k init , k = 1..n. At each iteration, a tree T ′ is chosen for expansion in a round-robin fashion among the trees T k , k = 1..n. Then, an extension is attempted toward a randomly sampled configuration q rand , starting from its nearest neighbor, q ′ near , in T ′ . We use an Extend function and not a Connect one, as recommended in [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF]. If the extension succeeds and the transition test is satisfied, the new configuration q new is added to T ′ and connected to q ′ near . Then, we look for the configuration q ′′ near (and the tree T ′′ containing it), which is the closest to q new within all trees other than T ′ . A connection between q new and q ′′ near is attempted in both directions, by calling twice the attemptLink function. The exploration continues until all Algorithm 3: Multi-T-RRT input : the configuration space C, the cost function c : C → R+ and the waypoints q k init , k = 1..n output: the tree T

1 for k = 1..n do 2 T k ← initTree(q k init ) 3 while not stoppingCriteria({T k | k = 1..n}) do 4 T ′ ← chooseNextTreeToExpand() 5 q rand ← sampleRandomConfiguration(C) 6 q ′ near ← findNearestNeighbor(T ′ , q rand ) 7 qnew ← extend(q ′ near , q rand ) 8 if qnew = null and 9 transitionTest(T ′ , c(q ′ near ), c(qnew)) then 10 addNewNodeAndEdge(T ′ , q ′ near , qnew) 11 (T ′′ , q ′′ near ) ← findNearestTree(qnew) 12 T ← attemptLink(T ′ , qnew , T ′′ , q ′′ near , n) 13 if T = null then 14 T ← attemptLink(T ′′ , q ′′ near , T ′ , qnew , n) Algorithm 4: attemptLink(T 1 , q 1 , T 2 , q 2 , n) input : the extension step-size δ output: the tree T 1 if distance(q1, q2) < 10 • δ then 2 qcur ← q1 ; qnext ← extend(q1, q2) 3 while qnext = null and c(qnext) ≤ c(qcur) do 4 qcur ← qnext ; qnext ← extend(qcur, q2) 5 if qcur = q2 then 6 T ← linkAndMerge(T1, q1, T2, q2) ; n ← n -1
trees are merged or another stopping condition (number of nodes, number of expansions, running time) is met.

The attemptLink function (shown in Algorithm 4) was developed to attempt connections between both trees of the Bidirectional T-RRT [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF]. If the configurations q 1 and q 2 are closer than ten times the extension step-size δ, and if the cost along the local path between them decreases monotonically (which is checked after every step of size δ), the trees T 1 and T 2 are merged, and the number of trees is decreased by 1. A distance threshold of 10 • δ represents a good trade-off between 1) wasting time checking edges that are unlikely to be valid (if the threshold is too high) and 2) having difficulties connecting trees (if the threshold is too low), as shown by the experiments presented in [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF].

Using the transition test of T-RRT and testing tree connections based on cost constraints enables the Multi-T-RRT to favor low-cost regions of the space, and thus to yield lowcost paths. The cost of a path can be defined in several ways based on the costs of the configurations along the path, as we will show in the examples.

B. Other Multi-Tree Variants

To develop the Multi-T-RRT we addressed several points for which we had to choose among various alternatives. The first point was to decide which tree(s) to expand at a given step of the exploration process. The simplest strategy is to grow all trees at each iteration toward the same configuration q rand [START_REF] Li | An incremental learning approach to motion planning with roadmap management[END_REF], [START_REF] Gayle | Lazy Reconfiguration Forest (LRF) -an approach for motion planning with multiple tasks in dynamic environments[END_REF]. By having a single tree grown at each iteration, chosen in a round-robin fashion [START_REF] Ettlin | Randomised rough-terrain robot motion planning[END_REF], [START_REF] Wang | Triple RRTs: an effective method for path planning in narrow passages[END_REF], the trees are expanded toward different configurations q rand , which appears to work better. Another strategy is to expand the tree that is the closest to q rand [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. However, when testing this approach, we have found that it can be difficult to expand trees that are growing close to the boundaries of the configuration space. A more sophisticated approach consists of choosing the tree to be expanded based on some probabilities that can be fixed [START_REF] Strandberg | Augmenting RRT-planners with local trees[END_REF] or adaptive [START_REF] Wang | An adaptive roadmap guided Multi-RRTs strategy for single query path planning[END_REF]. But, we have found that such strategies show no clear benefit in terms of improving running time or path quality.

The second point was to decide when to attempt to link trees. The simplest strategy is to try after each successful expansion of a tree [START_REF] Clifton | Evaluating performance of multiple RRTs[END_REF], [START_REF] Ettlin | Randomised rough-terrain robot motion planning[END_REF], [START_REF] Wang | Triple RRTs: an effective method for path planning in narrow passages[END_REF], [START_REF] Li | An incremental learning approach to motion planning with roadmap management[END_REF], [START_REF] Gayle | Lazy Reconfiguration Forest (LRF) -an approach for motion planning with multiple tasks in dynamic environments[END_REF], [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF]. Other, more sophisticated approaches consist of attempting a connection only when the bounding box of the expanded tree has increased in size [START_REF] Strandberg | Augmenting RRT-planners with local trees[END_REF], or when some stochastic test is satisfied, based on fixed or adaptive probabilities [START_REF] Wang | An adaptive roadmap guided Multi-RRTs strategy for single query path planning[END_REF]. However, we have found that these approaches lead to many missed good opportunities for connection.

The third point was about how to perform the link attempt after a tree has been successfully expanded. This can involve a single tree, usually the nearest one [START_REF] Ettlin | Randomised rough-terrain robot motion planning[END_REF], [START_REF] Wang | An adaptive roadmap guided Multi-RRTs strategy for single query path planning[END_REF], [START_REF] Zucker | Multipartite RRTs for rapid replanning in dynamic environments[END_REF], or a randomly chosen one. It can also involve all the other trees [START_REF] Clifton | Evaluating performance of multiple RRTs[END_REF], [START_REF] Strandberg | Augmenting RRT-planners with local trees[END_REF], [START_REF] Flavigné | Improving motion planning in weakly connected configuration spaces[END_REF], [START_REF] Li | An incremental learning approach to motion planning with roadmap management[END_REF], [START_REF] Gayle | Lazy Reconfiguration Forest (LRF) -an approach for motion planning with multiple tasks in dynamic environments[END_REF] or a subset of these trees, containing, e.g., some of the closest ones and some randomly chosen ones [START_REF] Plaku | Samplingbased roadmap of trees for parallel motion planning[END_REF]. When a tree is chosen for the link attempt, we have to decide which node in this tree we will try to connect the new node of the expanded tree to. Again, it can be the nearest one or a randomly chose one. After evaluation, we have found that random choices are not beneficial. It works better to attempt a connection between the new node and its nearest neighbor within the nearest tree.

C. Useful-Cycle Addition

Even though the paths the Multi-T-RRT returns have low cost, from a higher-level ordering perspective, they might not represent the most efficient way to visit a set of waypoints. To address this issue, we propose a simple approach based on the anytime paradigm and the addition of useful cycles: after all trees are connected, we allow the exploration to continue and we activate a cycle-addition procedure. This leads to the appearance of new paths that can be of better quality (with respect to the sequence of visited waypoints) than the one found so far. Adding cycles works as follows: When a new configuration q new is added to the graph, we consider all other configurations within a pre-specified distance in C as potential candidate targets for new edges. Among these candidates, we are interested in those that are "close" to q new in C, but "far" from q new in the graph: for each candidate q c , if the cost of the local path between q new and q c in C is strictly less than the cost of the lowest-cost path between q new and q c in the graph, we add an edge between q c and q new , thus creating a useful cycle. For more details on the anytime variant of T-RRT, and for a theoretical analysis, the interested reader is referred to [START_REF] Devaurs | Efficient sampling-based approaches to optimal path planning in complex cost spaces[END_REF].

IV. EVALUATION RESULTS

We have evaluated the Multi-T-RRT on several academic path planning problems that differ in terms of C-space dimensionality, geometrical complexity and cost-function type. We report results for three of them here. To fairly compare all algorithms, we first set aside the enhancement involving useful-cycle addition. For each example, we define ten waypoints that have to be visited in a pre-defined order (only to facilitate the evaluation of the algorithms). The Landscape problem is the 2D cost-map illustrated by Fig. 2, in which the cost is the elevation. The Stones problem (presented in Fig. 3) is a 2-degrees-of-freedom (DoF) problem in which a disk goes through a space cluttered with rectangular-shaped stones. The objective is to maximize clearance, so the cost function is the inverse of the distance between the disk and the closest obstacle. The Inspection problem (shown in Fig. 4) involves a 6-DoF manipulator arm holding a sensor with a spherical extremity, used to inspect a car engine. The objective is to keep the sensor as close as possible to the engine, so the cost function is the distance between the sphere and the engine surface.

The Multi-T-RRT has been implemented in the motion planning platform Move3D [START_REF] Siméon | Move3D: a generic platform for path planning[END_REF]. To fairly assess it, no smoothing is performed on the solution paths. On all problems, we record the running time t (in seconds), the number of expansion attempts X, the number of nodes N in the produced tree, and several quality criteria evaluated on the extracted path (with steps of size δ): the average cost avgC, the maximal cost maxC, the mechanical work M W , and the integral of the cost IC. The mechanical work of a path is the sum of the positive cost variations along this path [START_REF] Jaillet | Sampling-based path planning on configuration-space costmaps[END_REF]. For all variables, we give values averaged over 100 runs. Results were obtained on an Intel Core i5 processor at 2.6 GHz with 8 GB of memory.

We have compared two variants of the Multi-T-RRT to the basic version presented in Algorithm 3. We report the results of this comparison in Table I. The first variant involves having a local temperature associated to each tree, as opposed to having a global temperature. After evaluation, it seems that this modification has barely any influence on the results. The second variant is based on ensuring that all trees remain balanced (in terms of number of nodes) during the exploration. After evaluation, it is unclear whether this modification is advantageous or not. It appears to have sometimes a positive impact (e.g., on the Stones problem) and sometimes a negative impact (e.g., on the Landscape problem) on performance.

We have compared the Multi-T-RRT to the Bidirectional T-RRT [START_REF] Devaurs | Enhancing the transition-based RRT to deal with complex cost spaces[END_REF] in two ways. First, in a simple scheme involving the Bidirectional T-RRT, we compute paths between pairs of consecutive waypoints, starting from scratch each time, and concatenate them to obtain the full path visiting all waypoints. Second, in an incremental scheme involving the Bidirectional T-RRT, we compute paths between pairs of consecutive waypoints, while keeping the tree built so far instead of deleting it as in the simple scheme. Results obtained with these two schemes are reported in Table I. As expected, Therefore, the performance improvement is not achieved at the expense of path quality.

V. APPLICATION OF THE MULTI-T-RRT

In Section IV, we used the Multi-T-RRT to compute a high-quality path visiting an ordered set of waypoints, but only for evaluation purposes. In practice, the waypoints are not ordered a priori. Such ordering-and-pathfinding problems encompass two levels: a low-level path planning problem aiming at connecting the waypoints (which is solved by a geometric path planner) and a high-level ordering problem aiming at finding an efficient way to visit all the waypoints, based on the costs of the low-level paths (which is a simple kind of task planning problem). No symbolic task planner is required if we consider that the latter problem is an instance of the Traveling Salesman Problem (TSP) involving the complete graph (when path-cost is the integral of the cost) or digraph (when path-cost is the mechanical work) whose nodes are the waypoints. The distance associated with an edge of this graph (or digraph) can be estimated as the cost of the lowest-cost path between two waypoints in the tree built by the Multi-T-RRT. When only few waypoints are defined, the TSP is solved by an exhaustive search among all sequences. When more waypoints are involved, the Nearest-Neighbor or Multi-Fragment heuristics are used [START_REF] Bentley | Experiments on traveling salesman heuristics[END_REF].

Based on this approach, we present an industrial inspection problem involving an aerial robot in a dense environment, as illustrated by Fig. 1. This example is typical of those addressed by the ARCAS project (http://www.arcas-project.eu). One of the goals of this project is to develop robot systems for the inspection and maintenance of industrial installations difficult to access for humans. In this example, a quadrotor is used to inspect an oil rig, going through eight waypoints defined a priori without explicit order (cf. Fig. 1). The quadrotor is modeled as a 3-DoF sphere (i.e. a free-flying sphere) representing the security zone around it. For safety reasons, it has to move in this environment trying to maximize clearance. The cost function is thus the inverse of the distance between the quadrotor and the obstacles. Assuming that the motions of the quadrotor are performed quasistatically, we restrict the problem to planning in position (controllability issues lie outside the scope of this paper). Even though this example features a large-scale workspace, the Multi-T-RRT can quickly provide a first solution path: in about 5 s on average over 100 runs.

As already mentioned, a path produced by the Multi-T-RRT does not necessarily provide the best-quality solution visiting all the waypoints (see Fig. 3). The variant of the Multi-T-RRT involving useful-cycle addition was developed to enable continual improvement of the solution quality. As an example, Fig. 5 shows the achieved benefit on the Stones problem, especially when compared to Fig. 3. The path in Fig. 5 is representative of what we obtain for a running time of 5 s (as observed over 100 runs). Its cost is about half the cost of the path in Fig. 3 (obtained in 0.3 s). As another example, Fig. 1 shows a high-quality solution path obtained in 50 s on the Oil rig problem. Its cost is about half the cost of the first path obtained in 5 s.

To quantify the benefits of adding useful cycles, we evaluate the evolution of the quality of the solution path over time on the four examples (see Fig. 6). On the Landscape problem, path quality is measured by the mechanical work, M W , and on the other problems it is measured by the integral of the cost, IC. We compare the rate of convergence of the quality of paths produced by the Multi-T-RRT to the rates of convergence observed when planning with versions of PRM [START_REF] Nieuwenhuisen | Useful cycles in probabilistic roadmap graphs[END_REF] and RRT obst way [START_REF] Ettlin | Randomised rough-terrain robot motion planning[END_REF] creating cycles. Fig. 6 shows that the Multi-T-RRT yields a slightly better rate of convergence than RRT obst way . The poor results of PRM, especially in high-dimensional spaces, are due to the fact that it does not involve any cost constraint. As it features a similar convergence rate, IRS [START_REF] Marble | Towards small asymptotically near-optimal roadmaps[END_REF] would probably not perform better. RRT* [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF] would provide a better point of comparison, but it would require a multi-tree extension, which is out of the scope of this paper.

VI. CONCLUSION

We have presented a multi-tree variant of T-RRT named Multi-T-RRT. To achieve the highest efficiency, we have selected the best techniques involved in other multi-tree path planners. We have studied several path planning problems with different cost functions, geometrical complexity, and configuration-space dimensionality. When looking for a path going through an ordered list of waypoints, we have shown that the Multi-T-RRT is faster than planners based on the Bidirectional T-RRT, and that it yields paths of similar quality. When the planning problems involve both ordering and path planning aspects (i.e. when the order of the waypoints is not defined a priori), an anytime variant of the Multi-T-RRT enhanced with useful-cycle addition enables the solution to be continually improved. This allows us to visit the waypoints following a high-quality path computed (without the use of a symbolic task planner) based on the costs of the low-level paths connecting pairs of waypoints. Finally, we have applied the Multi-T-RRT to a realistic industrial inspection problem. We have also shown that the convergence rate of the solution quality is better with the Multi-T-RRT than with PRM and RRT obst way .

The approach we have proposed to solve ordering-andpathfinding problems in continuous cost spaces is a general one. This paper has focused on enhancing the T-RRT algorithm to solve such problems, by applying the multiple-tree, anytime and useful-cycle paradigms. It would be interesting to enhance other path planners in a similar way and compare their performance to that of the Multi-T-RRT. The poor results obtained with PRM highlight that such path planners should involve cost constraints. Therefore, a good candidate seems to be RRT* [START_REF] Karaman | Sampling-based algorithms for optimal motion planning[END_REF], in a multiple-tree version.

Fig. 1 .

 1 Fig. 1. Top: eight waypoints (shown in red and circled in blue) defined for a quadrotor (whose close-up is shown in yellow) inspecting an oil rig. Bottom: example of a trajectory (going through these waypoints) produced in 50 s by the Multi-T-RRT enhanced with useful-cycle addition.

Algorithm 1 : 3 q

 13 Transition-based RRT input : the configuration space C, the cost function c : C → R+ and the initial configuration qinit output: the tree T 1 T ← initTree(qinit) 2 while not stoppingCriteria(T ) do rand ← sampleRandomConfiguration(C) 4 qnear ← findNearestNeighbor(T , q rand ) 5 qnew ← extend(qnear , q rand ) 6 if qnew = null and 7 transitionTest(T , c(qnear), c(qnew)) then 8 addNewNodeAndEdge(T , qnear , qnew) Algorithm 2: transitionTest (T , c i , c j ) input : the current temperature T and the increase rate Trate output: true if the transition is accepted, false otherwise

Fig. 2 .

 2 Fig. 2. Search tree built by the Multi-T-RRT on the Landscape problem. On this 2D cost-map, the cost is color-coded (from blue to red) and represented by the elevation. The ten waypoints are materialized by red disks.

Fig. 3 .

 3 Fig. 3. Path computed by the Multi-T-RRT on the Stones problem, going through all ten waypoints. The cost is the inverse of the distance between the 2-DoF yellow disk and the blue rectangular-shaped obstacles.

Fig. 4 .

 4 Fig. 4. Ten waypoints defined for the Inspection problem. The 6-DoF manipulator arm holds a sensor (the red sphere) that has to follow the surface of the car engine.

Fig. 5 .

 5 Fig. 5. Path produced by the Multi-T-RRT enhanced with useful-cycle addition on the Stones problem, for a running time of 5 s. The labels show the order in which the waypoints are visited.

Fig. 6 .

 6 Fig. 6. Evolution of the quality of paths produced by the Multi-T-RRT, RRT obst way and PRM over time on problems of increasing dimensionality.

TABLE I EVALUATION

 I OF THE MULTI-T-RRT AND BIDIRECTIONAL T-RRT ON THE Landscape, Stones AND Inspection PROBLEMS.

			avgC	maxC	M W	IC	t (s)	N	X
		-basic	11	22	240 10,000	0.06 1,100	6,000
	Multi-T-RRT	-balanced	11	22	230	9,900	0.18 1,300 13,000
	Landscape	-local temperature	11	22	240 10,000	0.08 1,200	9,000 Average values over 100 runs
	simple Bidirectional T-RRT	11	22	230	9,800	0.12 2,700 20,000	are given for:
	incremental Bidirectional T-RRT	11	22	240	9,900	0.12 1,400 10,000
		-basic	2.4	5.9	110 27,000	0.38 2,000 18,000	avgC	average cost
	Multi-T-RRT	-balanced	2.4	5.9	110 29,000	0.28 1,100	8,000 maxC	maximal cost
	Stones	-local temperature	2.5	5.9	110 29,000	0.38 2,100 18,000	M W	mechanical work
	simple Bidirectional T-RRT	2.3	5.9	97	25,000	0.57 4,500 38,000	IC	integral of the cost
	incremental Bidirectional T-RRT	2.2	5.9	104 28,000	0.43 2,100 18,000
		-basic	4.4	19	470 20,000	0.8	500 14,000	t	running time
	Multi-T-RRT	-balanced	4.1	19	480 20,000	0.9	600 16,000	N	number of nodes
	Inspection	-local temperature	4.3	19	470 19,000	1	600 18,000	in the tree
	simple Bidirectional T-RRT	3.8	19	400 16,000	2 1,100 39,000	X	number of
	incremental Bidirectional T-RRT	3.7	19	540 20,000	1.5	700 28,000	expansion attempts
	the Multi-T-RRT is faster than the planners involving the
	Bidirectional T-RRT. Moreover, in spite of performing a
	quicker exploration of the space, the Multi-T-RRT produces
	paths whose costs are only slightly worse than those of paths
	produced by the planners involving the Bidirectional T-RRT.