Fawzi Nashashibi -Ieee Member

Paulo Resende
email: paulo.resende@inria.fr

François Charlot
email: francois.charlot@inria.fr

Carlos Holguin

Michel Parent

Laurent Bouraoui
email: laurent.bouraoui@mines-paristech.fr

A Cooperative Personal Automated Transport System A CityMobil Demonstration in Rocquencourt

Keywords: CityMobil, Cybus shuttles, cooperative navigation, vehicle coordination, intelligent transportation, urban transport. I

This article tackles the problem of the autonomous navigation and coordination of multiple driverless vehicles for the transport of persons or goods in outdoor environments. The system composed of fully automated road vehicles, capable of providing an effective transportation service, was recently tested at the city of La Rochelle. This same system was further improved, and a new demonstration was performed at Inria Rocquencourt, in order to demonstrate the validity of the concepts for a coordinated navigation in the presence of ambiguous and conflictual situations in a mixed environment. The originality of the approach relies on the use of new cooperative concepts and their combination with advanced perception tasks operating simultaneously on several robots. This system was developed in the context of the European project CityMobil.

INTRODUCTION

The main objective of the CityMobil project was to design and develop new advanced concepts for efficient urban transportation capable of reducing energy consumption, optimizing used space, reducing pollution, improving security, comfort and quality of living. In order to educate the urban population to these concepts, scrutinize public acceptance and evaluate new transport modes in actual operating conditions, the project introduced systems composed of fully automated vehicles in the heart of some pilot cities bringing them closer to public observation.

It was in this context that between May and July 2011, Inria putted into service an on-demand door-to-door transport system in the city of La Rochelle [START_REF] Bouraoui | An On-demand Personal Automated Transport System: The CityMobil Demonstration in La Rochelle[END_REF], using fully automated electric and communicating road vehicles, better known as cybercars [START_REF] Parent | Cybercars: A solution for urban transport[END_REF]. These cybercars were demonstrated to the city population and their performance evaluated during the period of 3 months service. This demonstration showed how autonomous road vehicles can be effectively used in public transport, especially in low to medium demand areas and periods, to make the public transport service more reliable and frequent where conventional collective public transport cannot meet the mobility needs of the population in an environmentally and economically sustainable manner.

A new demonstration, organized after the end of the service in La Rochelle, was performed to validate the concept of collaboration between autonomous vehicles. This paper presents this new demonstration that took place at Inria's Rocquencourt campus in November 2011. As in La Rochelle, this site has presence of pedestrians and limited road traffic. However, unlike La Rochelle, Inria's campus is a private site with restricted access conditions. Therefore, this extended demonstration was not accessible to the general public.

The paper is composed of two main sections: the first will describe the system architecture from both the vehicles and the infrastructure elements. The second section will describe the operational aspects on the site. Results and conclusions will be provided at the end of the paper in the section III.

II. SYSTEM ARCHITECTURE

A. Overview

An on-demand door to door transportation service accessible to all the users of the INRIA Rocquencourt campus was implemented using fully automated electric and communicating vehicles. The vehicles behavior is similar to a horizontal elevator that can be called from a station and serve multiple destinations.

The monitoring and coordination of the vehicles was performed by a vehicles management system (VMS) through a dedicated wireless IPv6 mesh network [START_REF] Choi | IPv6 support for VANET with geographical routing[END_REF].

B. Site description

The demonstration site is a 365 m long two-lane road, where of 5 stations equipped with specially conceived wireless routers were installed.

The distance between the stations is as follows: In order to guarantee the wireless coverage over the whole site, an additional wireless router installed between Stop 3 and Stop 4, and served as a relay. For practical reasons, the "End" stop was not used during the demonstration.

All the stations are located in only one side of the road, to maintain the same site configuration that was implemented in La Rochelle, due to its physical infrastructure constraints (not enough space for the vehicles to turn around in both ends).

The Figure 1 shows an aerial view of the demonstration site.

C. Stations

Five stations where installed in the demonstration site where the users could call a vehicle via a touch screen interface. This screen also provided information such as the vehicles availability or the pick-up arrival time. All the stations where equipped with a wireless IPv6 communication device with corresponding antenna that relayed the data exchanged between the stations, vehicles and the vehicles management center.

Additional wireless relays needed to be added between some of the stations to have full coverage of communications in the demonstration area. These stations are installed only on one side of the road, forcing the vehicles to change lanes in order to approach the station platform.

D. Vehicle Platform

The vehicles, named Cybus, are represented on Fig. 2. These vehicles have the capacity to transport 4 or 5 seated persons at an average speed of 8 km/h and with a maximum speed of 12 km/h. They are basically the same than the ones used during the demonstration of La Rochelle in 2011.

Thanks to the gained experience during the 3 months of this demonstration, the vehicles were modified in order to improve control, localization and safety.

1.

Control During the La Rochelle demonstration, tedious oscillations were noticed when the Cybus was moving on long straight trajectories. After numerous trials and experiments, it has been shown that the low-level control algorithm was not responsive enough and that the current implementation was error-prone. Once this incorrect behavior was identified, the control system has been reviewed and re-configured in less than a day thanks to the highly modular design of the whole hardware/software system.

2.

Localization Initially the localization and obstacles detection tasks were performed by using the same 2 Ibeo Alasca XT laser sensors connected to a single fanless computer via an Ibeo ECU dedicated to read, synchronize, calibrate and merge the lasers data. This configuration used during the demonstration in La Rochelle highlighted some problems that were not experimented before.

First of all, the whole system needed to be restarted too often due to the lack of reliability of the ECU. Despite different changes in the ECU configuration, it was not possible to solve this problem with such kind of sensor configuration.

In addition, the laser based SLAM method [START_REF] Xie | A real-time robust global localization for autonomous mobile robots in large environments[END_REF] was not able to find the vehicle position in a very crowded environment. This was due to the obstruction of the lasers by people surrounding the vehicle as the lasers position was at less than 1 meter from the ground. In such situation, the current local map cannot match the recorded global map.

To solve these problems, a new sensor and processing configuration was implemented: two Ibeo Lux laser scanners, connected to a dedicated computer, have been integrated on the roof of the Cybus (around 2m height) at the front and at the rear ends. They provided from 2 to 4 layers range data with a maximum view angle of 110° and a maximum range of 200 A RTK (Real Time Kinematics) GPS has also been integrated into the vehicle to provide ground truth and to evaluate the accuracy of the localization.

3.

Safety Autonomous vehicles may need to stop in a progressive way in the case of obstacles in the way, or to stop immediately in the case of an imminent collision or sudden vehicle malfunction. It is for those reasons that the Cybus vehicles are equipped with several safety devices, managed by the highlevel software as well as by the low-level embedded electronics.

In the low-level, five emergency stop buttons, connected in series and arranged on each side of the vehicle and near the control screen, may be pushed at any time to immediately bring the vehicle to a full stop. Furthermore, each vehicle has a retractable front bumper with two emergency switches. These two switches are also connected in series with the emergency buttons. Finally, every vehicle is equipped with a radiocontrolled emergency cut-off remote of category 3 according to the French norm NF EN954-1for safety devices.

An emergency stop is activated when the emergency remote button is engaged or there is a loss of communication between the remote and the vehicle. The activation of the emergency stop implies the cut off the power supply of the drive components and the release of the vehicle parking brake, resulting in an immediate stop.

In the high-level software, obstacle detection is performed using 2 Alasca XT lasers scanners connected to an Ibeo ECU. The algorithm of detection defines a customizable safety area, depending among others on the driving direction, the speed and the dimensions of the vehicle. In the present of objects in this safety area, the vehicle will decrease its speed and eventually stop if the object gets closer to the vehicle. The lack of reliability of the Ibeo ECU still caused some problems to the system, but the system safety was not greatly affected.

E. Communications

Communications are an essential element in automated systems. They allow the cooperation between systems by the exchange of relevant information for an augmented perception, planning and control of these systems.

By means of communications we can have for example, more consistent and complete local and global maps, data collection, user requests management and anticipation, avoids unnecessary stops and deadlocks, or improves the traffic flow.

A specific multi-hop Mobile Ad hoc Network (MANET) was used to perform V2V and V2I wireless communications using standardized modules compliant with the IPv6 protocol. These modules where packaged into a dedicated hardware called Cube that was specially designed and conceived for ITS applications. The Cubes, shown in Figure 4, form the nodes of the mesh network. They handle all the data exchanged between the automated vehicles and the infrastructure elements such as the stations or the control center (VMS).

The Cube is a Linux embedded PC with 2 mini PCI communication cards. The radio frequency used was in the licensed ITS band of 5.9 GHz of the standard 802.11p for wireless access in vehicular environments (WAVE). The implemented services discovery based architecture provided an IPv6 network based on a B.A.T.M.A.N. mesh network and a Cables framework that manages the services. The diagram in Figure 5 is a representation of the mesh network and the data quality in the links between nodes. All traffic is encapsulated and forwarded until it reaches the destination, hence emulating a virtual network switch of all nodes participating. Therefore, all nodes appear to be a local link and are unaware of the network's topology and are unaffected network changes.

F. Vehicles Management System

The Vehicle Management System (VMS) it is an application handles the stations requests, the conflict areas and monitors the automated vehicles. The information about the state of these vehicles is stored in a pool and kept up to date by the VMS. The knowledge of the vehicles environment is limited to the conflict areas (mainly the stations) and the vehicles path. From the VMS an operator can stop all the vehicles at once. In this implementation the management of the stations and vehicles is centralized. Certain decisions such as the coordination between 2 vehicles could be decentralized with one of them taking the leading role and making the decision.

The VMS, represented in Figure 6, is composed of two distinct parts:

A web server that provides a user interface to the operator and stations, and an API that allows a vehicle to be called from the station or by a mobile device.

A server that communicates with the vehicles using a library called Cables.

The interface between the two applications is done through a common database.

1.

Station calls management Vehicle calls at the stations are done through the stations graphical interface and are stored in a central database. These dynamic web pages are handled by the VMS. Once there is a new call in the database the attribution of the station to a vehicle it is done by choosing the vehicle closest to the calling station and that, preferably, already travels in the direction of that station. Then, the VMS sends the station position as a new destination the selected vehicle that will add it to its destination queue in the correct order.

2.

Vehicles management A server communicates directly with the vehicles. This server monitors and manages the vehicles priorities when approaching conflict areas (e.g. stations) and when overtaking other communicating vehicles.

By setting priorities to access a conflict area and by automatically redistributing the vehicles along the network according to the demand it avoid deadlocks and optimizes the vehicles availability time. If only obstacle detection was used and no coordination was performed, deadlocks would happen. Besides, by doing an automatic redistribution of the vehicles, the waiting time is reduced in the stations with the highest demand.

Supervision

The VMS allows an operator to monitor and partially control the entire fleet of vehicles, stations and communications network. For example, an operator can send a command to stop or slow down all vehicles. The stations webbased graphical interfaces can be modified directly by the operator to display some specific information to the users at specific station. Through the VMS the operator can also see the status of the secured wireless mesh network and perform the adequate action in the case, for example, of a faulty node.

The operator can always override the automatic supervisor of the VMS. For example, the operator may disable all overtaking manoeuvres, or send predefined destinations to a specific vehicle (e.g. depot or charging station).

Through a web interface an operator could easily monitor the status and position of all the vehicles in the network.

III. SYSTEM OPERATION

The demonstration performed at the INRIA Rocquencourt campus consisted of an on-demand door-to-door transport service using two Cybus, .in a circuit with two lanes with opposite traffic directions like represented in the following figure :   The station's where positioned in only one side of the road that leads to conflictive situations between the vehicles. The Several scenarios of conflicts between the driverless vehicles and other non-communicating obstacles were considered but only partially solved. An example of these conflicts is the presence of a truck or pedestrian obstructing the vehicle path. To solve this scenario in an optimal way a more complex set of sensors, then the one available in the vehicles would be required. In particular, it would be necessary to determine the type of obstacle and its properties (e.g. length, speed) so that the driverless vehicle could perform the adequate overtaking manoeuvre respecting at the same time the road traffic rules. With the current implementation and given the previous situation, the Cybus would just stop before the obstacle and wait for further instructions from the VMS.

The conflict management between the communicating driverless vehicles was centralized in the VMS although the vehicles add enough intelligence on board to solve these types of conflicts by themselves.

In order to automatically redistribute the automated vehicles in the circuit it would be necessary to know the vehicles occupancy and the doors status (closed or opened).

In the absence of these features the vehicles where redistributed by an operator inside the vehicle. Some representative conflict scenarios between communicating driverless vehicles were identified and solved:

Vehicle stopped at station

In this scenario, a vehicle 2 is stopped at a station.

In the situation where another vehicle 1 goes to the same station where vehicle 2 is stopped, the VMS will order the vehicle 1 to stop and wait for vehicle 2 to leave the station before allowing it to resume its movement. See Figure 8-a.

In the situation where vehicle 1 passes, without stopping, to the station where vehicle is stopped, the VMS will order vehicle 2 to remain stopped at the station and wait for vehicle 1 to pass. See Figure 8-b.

Vehicle leaving station

In this scenario, a vehicle 1 is leaving a station. Another vehicle 2 needs to stop until vehicle 1 has left the station by making a lane change maneuver. Only after this vehicle 2 can advance to the free station.

Vehicle arrive to the same station

In this scenario, a vehicle 1 and a vehicle 2 drive at the same time towards the same destination station. In this case, the VMS will determine, based on the proximity, which vehicle will have the priority. The closest vehicle will be allowed to continue and dock at the station, while the other vehicle will be ordered to stop and wait.

Vehicle arrive to the same station

In this scenario, a vehicle 2 has to overtake a vehicle 1 stopped in the same lane (not necessarily at a station). In this case, the VMS will provide vehicle 2 the authorization to overtake and order vehicle 1 to remain stopped until vehicle 2 has finished the overtaking manoeuvre. The vehicle 2 will then overtake vehicle 1 based on the information that vehicle 1 provides (e.g. pose and dimensions).

IV. CONCLUSION

We presented in this paper the extensions done at the Inria Rocquencourt demonstration for the European project CityMobil that allowed the collaboration between driverless vehicles. These modifications were made to the vehicles management system and certain elements of vehicles architecture to allow the vehicles cooperation. These vehicles were the same used in the on-demand personal automated transport system previously demonstrated in La Rochelle also for the CityMobil project.

The consecutive successful technological demonstrations performed and the very positive public acceptance demonstrates that these systems are a valid mobility solution for the cities of tomorrow and can already be used in today's private sites.

The vehicles cooperation is performed using a centralized architecture where the VMS plays the decisional role. The VMS communicates with all the operating vehicles and resolves all conflict situations and bottlenecks. However in some situations it may be necessary that the vehicles rely only on their own decisional system especially in case of a VMS shortage or while performing a reactive navigation (e.g. obstacle avoidance or overtaking). We expect to further develop this system and explore the use of a decentralized decision system to solve these later issues.

Stop 1 (

 1 Entrance) -Stop 2 (Building 1): 50m Stop 2 (Building 1) -Stop 3 (Building 7): 70m Stop 3 (Building 7) -Stop 4 (Building 25): 155m Stop 4 (Building 25) -Stop 5 (End): 90m

Fig 2 .

 2 Fig 2. Cybus vehicles at Inria Rocquencourt

Fig 5 .Fig 4 .Fig 3 .

 543 Fig 5. Mesh network diagram

Fig 7 .

 7 Fig 7. Two lane traffic schema at INRIA Rocquencourt

Fig 10 .

 10 Fig 10. Vehicle 1and vehicle 2 arriving in the same station

ACKNOWLEDGMENT

This work was supported in part by the European projects CityMobil and CityNetMobil.


Bezza Hafidi
email: b.hafidi@ucam.ac.ma

Nourddine Azzaoui

Criteria for longitudinal data model selection based on Kullback's symmetric divergence

Keywords: Model selection, AIC, KIC, Longitudinal data analysis, Kullback's symmetric divergence Sélection de modèles, AIC, KIC, Analyse des données longitudinales, divergence symétrique de Kullback

showed that (AIC) criterion and its corrected versions cannot be directly applied to model selection for longitudinal data with correlated errors. They proposed two model selection criteria, AICc and RICc, by applying likelihood and residual likelihood approaches. These two criteria are estimators of the Kullback-Leibler's divergence distance which is asymmetric. In this work, we apply the likelihood and residual likelihood approaches to propose two new criteria, suitable for small samples longitudinal data, based on the Kullback's symmetric divergence. Their performance relative to others criteria is examined in a large simulation study. RÉSUMÉ. Récemment, Azari et al. (2006) ont montré que le critère (AIC) ainsi que ses versions corrigées ne peuvent pas être directement appliqués aux données longitudinales avec des erreurs corrélées. Ils ont proposé deux critères, AICc et RICc, en utilisant la notion de vraisemblance et de vraisemblance résiduelle. Leurs critères sont des estimations de la distance asymétrique de divergence de Kullback-Leibler. Dans ce travail, nous proposons deux nouveaux critères adaptés aux échantillons de petites tailles de données longitudinales en se basant sur la divergence symétrique de Kulback et les approches de maximum de vraisemblance et de vraisemblance résiduelle. Les performances de ces critères sont examinées dans une étude de simulations.

Introduction

Model selection criteria play an important role in applied statistical data analysis especially in parametric and non parametric regression, mixture models, time series analysis... The most known is the Akaike Information Criterion (AIC) which is applicable in different arrays of modeling, see Akaike [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]. The AIC criterion was designed as an unbiased estimator of a variant of the Kullback's directed divergence between the model which presumably gave rise to the data and a fitted approximating one. The directed divergence, also known as the Kullback-Leibler information, is a measure of separation between two statistical models. It is asymmetric, meaning that an alternate directed divergence can be obtained by reversing the roles of the two models in the definition of the measure. The sum of the two directed divergences is the Kullback's symmetric divergence, which combines the information in both measures. Therefore, it functions as gauge of model disparity which is arguably more balanced than either of individual directed divergences. Model selection criteria based on the symmetric measure was first investigated by Cavanaugh [START_REF] Cavanaugh | A large-sample model selection criterion based on Kullback's symmetric divergence[END_REF]. He proposed an Akaike type criterion, named KIC, as an asymptotically unbiased estimator of the symmetric divergence in the case of large sample data. When the sample size is small or the number of fitted parameters is large to moderate fraction of the sample size, both AIC and KIC criteria suffer from a negative bias and result in a serious overfitting. In order to overcome this problem, many corrected versions of AIC was proposed in different special situations. For instance, Hurvich and Tsai [START_REF] Hurvich | Regression and time series model selection in small samples[END_REF] proposed a corrected AIC (AICc) for linear and non-linear regression and autoregressive modeling. Their work was extended to the case of autoregressive moving average by Hurvich [START_REF] Hurvich | Improved estimators of Kulback-Leibler information for autoregressive model selection in small samples[END_REF], to vector autoregressive modeling by Hurvich and Tsai [START_REF] Hurvich | A corrected Akaike information criterion for vector autoregressive model selection[END_REF] and multivariate regression by Bedrick et al. [START_REF] Bedrick | Model selection for multivariate regression in small samples[END_REF], Fujikoshi and Satoh [START_REF] Fujikoshi | Modified AIC and Cp in multivariate linear regression[END_REF].... On the other hand, using the Kullback's symmetric divergence Cavanaugh [START_REF] Cavanaugh | Criteria for Linear Model Selection Based on Kullback's Symmetric Divergence[END_REF] proposed a corrected Kullback's information criterion (KICc) for linear models. Hafidi and Mkhadri [START_REF] Hafidi | A corrected Akaike criterion based on Kullback's symmetric divergence: applications in time series, multiple and multivariate regression[END_REF] generalized the corrected KICc for multiple to multivariate regression and for univariate or multivariate autoregressive modeling. Longitudinal data analysis has been a great deal of interest in the fields of clinical trials, epidemiology, agriculture and medicine over the last decade. These data arise when repeated measurements are obtained for an individual or more outcome variables at successive points in time. The successive measurements for each individual tend to be correlated. Therefore, it is necessary to take into account this correlation in order to produce proper analysis. The comprehensive synthesis of both theoretical, applied aspects, model structure details and parameter estimation about longitudinal analysis is given in Diggle [START_REF] Diggle | Analysis of Longitudinal Data[END_REF], Hedeker and Gibbons [START_REF] Hedeker | Applied Longitudinal Data Analysis[END_REF], Fitzmaurice et al. [START_REF] Fitzmaurice | Applied Longitudinal Analysis[END_REF], Hand and Crowder [START_REF] Hand | Practical Longitudinal Data Analysis[END_REF] and Jones [START_REF] Jones | Longitudinal Data With Serial Correlation: a state-space approach[END_REF]. Recently Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF] showed that the classical corrected AICc criteria can not be directly applied to model selection for longitudinal data with correlated errors. They derived two model selection criteria: the first one is obtained by applying the maximum likelihood approach and the second is RICc derived by using the residual (restricted) likelihood approach. Both these two criteria are estimator of the Kullback-Leibler's divergence distance which is asymmetric. Cavanaugh [START_REF] Cavanaugh | A large-sample model selection criterion based on Kullback's symmetric divergence[END_REF][START_REF] Cavanaugh | Criteria for Linear Model Selection Based on Kullback's Symmetric Divergence[END_REF] suggested that the Kullback's symmetric distance is a preferable tool for model selection than the asymmetric one. In this paper, we use the kullback's symmetric divergence and we apply the likelihood and residual likelihood approaches to derive corrected version suitable for small sample longitudinal data. A large simulation investigations for these criteria show their appropriateness for model selection especially when the sample size and the signal to noise ration are small. In this case, they can be used as an alternative to classical criteria. The paper is organized as follows. Section 2 is devoted to the description of model structure of longitudinal data. Derivation of the corrected KIC c , is presented in Section 3. In Section 4, we presented the criteria obtained by using the residual likelihood. We end the paper with simulation results and with a small conclusion.

Preliminaries and notations

Let y ij represents the measurement for the i th subject at the time points t ij . we denote by Y i = (y i1 , . . . , y ini ) t a (n i × 1) vector of all repeated observations on the i th subject. The sign ' t ' stands for the transpose. For simplicity, we assume that n i = n for all i. Suppose that the true model is given by

Y i = X i0 β 0 + ε i i = 1, . . . , m, [1]
where ε i = (ε i1 , . . . , ε in ) is a (n × 1) vector at time n. It is assumed to follow a multivariate normal distribution with mean zero and covariance matrix σ 2 0 Σ 0 . β 0 is a (p 0 × 1) vector of unknown parameters. Finally X i0 is an (n × p 0 ) matrix with the j th row (x 0 ij1 , . . . , x 0 ijp 0

). Define the observation (N × 1) vector Y as Y = (Y t 1 , . . . , Y t m ) t where N = nm. Let X 0 = (X t 10 , . . . , X t m0 ) t be a (N × p 0 ) matrix of explanatory variables and ε 0 = (ε t 1 , . . . , ε t m ) t the error vector. With this notation, we can form a general regression model Y = X 0 β 0 + ε 0 , where ε 0 ∼ N (0, σ 2 0 V 0 ) and V 0 is an (N × N ) block diagonal matrix with (n × n) blocks Σ 0 . In practice, we do not know the true model described below, thus we fit the data to a candidate one

Y i = X i β + ε i i = 1, . . . , m, [2]
where β is a (p × 1) vector and X i is an explanatory (n × p) matrix with j th row (x ij1 , . . . , x ijp ) and ε = (ε t 1 , . . . , ε t m ) t is the vector error with ε i ∼ N (0, σ 2 Σ). As previous, the model [START_REF] Azari | Longitudinal data model selection[END_REF] can be written as

Y = Xβ + ε,
where Y is a dependent (N × 1) vector, X is an (N × p) matrix of explanatory and ε ∼ N (0, σ 2 V ) and V is an (N × N ) block diagonal matrix with (n × n) blocks Σ. With these previous considerations, we denote θ = (β, σ 2 , Σ, p) and θ 0 = (β 0 , σ 2 0 , Σ 0 , p 0 ) the parameter vectors of the candidate and the true models respectively. Let f θ0 (Y ) and f θ (Y ) represent the generating and candidate densities for data respectively. In this case, ignoring the constant term, the log-likelihood for the candidate model is given by

log f θ (Y ) = - 1 2 N log σ 2 + m log |Σ| + (Y -Xβ) t V -1 (Y -Xβ) σ 2 . [3]
Similarly the log likelihood function of the true model is obtained by replacing θ, σ, β, Σ and V in equation ( 3) by θ 0 , σ 0 , β 0 , Σ 0 and V 0 respectively.

REMARK. -It should be noted that, in longitudinal data analysis, the covariance matrix Σ is often assumed to have a special morphology. For example:

-The uniform structure: where Σ = (σ jl ) j,l=1...n is defined by

σ jl = ρ for t ij = t il 1 for t ij = t il
and ρ is a positive correlation coefficient between two measurements on the same subject.

the exponential correlation structure: in this case σ jl = exp(-γ|t ij -t il |) where γ is the correlation decay rate between two measurements on the same unit.

-The autoregressive correlation structure: it is a particular case of the latter where the observation times are equally spaced for all j. In this setting σ jl can be expressed

as σ jl = ρ |j-l| or σ jl = ρ |j-l| 1-ρ 2 .
The parameter ρ is the correlation between successive observations on the same subject. In this paper we assume that the correlation structure of Σ and Σ 0 is expressed as a function of an unknown parameter vector φ and φ 0 respectively: meaning that Σ = Σ(φ) and Σ 0 = Σ 0 (φ O ). As Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF], we assume the consistency of the maximum likelihood and the restricted maximum likelihood estimators of Σ. In all the rest of the paper we replace the precedent notations of θ and θ 0 by θ = (β, σ 2 , φ) and θ 0 = (β 0 , σ 2 0 , φ 0 ).

Derivation of the corrected (KICc) based on the Kulback's symmetric divergence

A measure of separation between the generating and a candidate model is given by the symmetric divergence (Kullback [18]). It is defined by

J(θ 0 , θ) = {d(θ 0 , θ) -d(θ 0 , θ 0 )} + {d(θ, θ 0 ) -d(θ, θ)}, [4]
where

d(θ 0 , θ) = IE θ0 {-2 log f θ (Y )
} and IE θ0 denotes the expectation with respect to f θ0 (Y ). Ignoring d(θ 0 , θ 0 ) which does not depend on θ, for the purpose of discriminating among various models depending on θ we propose a substitution of J(θ 0 , θ) as,

K(θ 0 , θ) = d(θ 0 , θ) + {d(θ, θ 0 ) -d(θ, θ)}. [5]
In practice, the exact computation of this quantity is not easily accessible. To overcome this problem, Cavanaugh [START_REF] Cavanaugh | A large-sample model selection criterion based on Kullback's symmetric divergence[END_REF] proposed an asymptotically unbiased estimator of

Ω(θ 0 ) = IE θ0 {K(θ 0 , θ)}. [6]
in the case of large sample data, where θ is the maximum likelihood estimator of θ, K(θ 0 , θ) has the same expression as in ( 5) by replacing θ by θ. This estimator is given by

KIC = N log σ2 + m log | Σ| + 3(p + 1). [7]
The purpose of this section is to derive a corrected version of this criterion. For this, we assume that the candidate class of models includes the true model. This assumption is also used in the derivation of KIC [START_REF] Cavanaugh | A large-sample model selection criterion based on Kullback's symmetric divergence[END_REF] and its other versions see for instance [START_REF] Azari | Longitudinal data model selection[END_REF][START_REF] Cavanaugh | Criteria for Linear Model Selection Based on Kullback's Symmetric Divergence[END_REF][START_REF] Hafidi | A small-sample criterion based on Kullback's symmetric divergence for vector autoregressive modeling[END_REF][START_REF] Hafidi | A corrected Akaike criterion based on Kullback's symmetric divergence: applications in time series, multiple and multivariate regression[END_REF].

In this setting the columns of X can be rearranged so that X 0 β 0 = Xβ * , where β * = (β t 0 , β t 1 ) t and β 1 is a (p -p 0 ) × 1 vector of zeros. Under this assumption, we show the forthcoming proposition.

Proposition 1 The criterion defined by

KICc = N log σ2 + m log | Σ| + (p + 1)(3N -p -2) N -p -2 [8]
is an approximate unbiased estimator of Ω(θ 0 ).

Proof: Ignoring the constant term, the log likelihood of the candidate model is given by

log f θ (Y ) = - 1 2 N log σ 2 + m log |Σ| + (Y -Xβ) t V -1 (Y -Xβ) σ 2 . [9]
The maximum likelihood estimators of β, σ 2 and V (φ) are given by

β = (X t V -1 X) -1 X t V -1 Y, σ2 = (Y -X β) t V -1 (Y -X β) N , [10]
where V -1 is V -1 (φ) evaluated at φ. This latter is obtained by maximizing the profile log-likelihood:

(φ) = - 1 2 (N σ2 (φ) + m log |Σ(φ)|)
, see for instance Azari [START_REF] Azari | Longitudinal data model selection[END_REF] and for more details see Diggle ([7] p. 64-65). From (5), the expectation of K(θ 0 , θ) has the form

Ω(θ 0 ) = IE θ0 {d(θ 0 , θ) + d( θ, θ 0 ) -d( θ, θ)}.
[11]

For clarity, we emphasize that

d(θ 0 , θ) = IE θ0 {-2 log f θ (Y )} θ= θ and d( θ, θ 0 ) = IE θ {-2 log f θ0 (Y )} θ= θ .
The notation θ= θ means that we calculate the expectation and then we replace θ by θ in the resulting formula. Now, we will compute each term of [START_REF] Hafidi | A corrected Akaike criterion based on Kullback's symmetric divergence: applications in time series, multiple and multivariate regression[END_REF]. By adding and subtracting Xβ * , we have

d(θ 0 , θ) = IE θ0 {-2 log f θ (Y )} θ= θ = IE θ0 N log σ 2 + m log |Σ| + (Y -Xβ) t V -1 (Y -Xβ) σ 2 θ= θ = IE θ0 N log σ 2 + m log |Σ| + (Y -Xβ * ) t V -1 (Y -Xβ * ) σ 2 + (β -β * ) t X t V -1 X(β -β * ) σ 2 θ= θ .
It is easy to see that the expectation IE θ0 {N log σ 2 } θ= θ = N log σ2 and similarly we have

IE θ0 {m log |Σ|} θ= θ = m log | Σ|.
By using the same reasoning we also have,

IE θ0 (β -β * ) t X t V -1 X(β -β * ) σ 2 θ= θ = ( β -β * ) t X t V -1 X( β -β * ) σ2 . On the other hand, since (Y -Xβ * ) t V -1 (Y -Xβ * ) = ε t 0 V -1 ε 0 , consequently we deduce that IE θ0 (Y -Xβ * ) t V -1 (Y -Xβ * ) σ 2 θ= θ = σ 2 0 σ2 tr( V -1 V 0 )
, where "tr" stands for trace. Finally,

d(θ 0 , θ) = N log σ2 + m log | Σ| + σ 2 0 σ2 tr( V -1 V 0 ) + ( β -β * ) t X t V -1 X( β -β * ) σ2 . [12]
Similarly, we have

d( θ, θ 0 ) = IE θ {-2 log f θ0 (Y )} θ= θ = IE θ N log σ 2 0 + m log |Σ 0 | + (Y -Xβ * ) t V -1 0 (Y -Xβ * ) σ 2 0 θ= θ = IE θ N log σ 2 0 + m log |Σ 0 | + (Y -Xβ) t V -1 0 (Y -Xβ) σ 2 0 + (β -β * ) t X t V -1 0 X(β -β * ) σ 2 0 θ= θ . [13]
With the same technique as in [START_REF] Hand | Practical Longitudinal Data Analysis[END_REF], the expectation calculation leads to the equality

IE θ (Y -Xβ) t V -1 0 (Y -Xβ) σ 2 0 θ= θ = IE θ ε t V -1 0 ε σ 2 0 θ= θ = σ2 σ 2 0 tr( V V -1 0 ). This implies that d( θ, θ 0 ) = N log σ 2 0 + m log |Σ 0 | + σ2 σ 2 0 tr( V V -1 0 ) + ( β -β * ) t X t V -1 0 X( β -β * ) σ 2 0 . [14]
The same reasoning leads to

d( θ, θ) = IE θ {-2 log f θ (Y )} θ= θ = IE θ N log σ 2 + m log |Σ| + (Y -Xβ) t V -1 (Y -Xβ) σ 2 θ= θ = IE θ N log σ 2 + m log |Σ| + ε t V -1 ε σ 2 θ= θ = N log σ2 + m log | Σ| + N. [15]
Substituting ( 12), ( 14) and ( 15) into [START_REF] Hafidi | A corrected Akaike criterion based on Kullback's symmetric divergence: applications in time series, multiple and multivariate regression[END_REF] we obtain

Ω(θ 0 ) = IE θ0 N log σ2 + m log | Σ| + σ 2 0 σ2 tr( V -1 V 0 ) + ( β -β * ) t X t V -1 X( β -β * ) σ2 [16] +IE θ N log σ 2 0 σ2 + m log |Σ 0 | | Σ| + σ2 σ 2 0 tr( V V -1 0 ) + ( β -β * ) t X t V -1 0 X( β -β * ) σ 2 0 -N . [17]
REMARK. -Under the assumption that φ is a consistent estimator of φ 0 , Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF] approximate V by V 0 ; i.e. V = V 0 + o p (1). They have given a corrected version of the Akaike criterion and they showed that

AICc = N log σ2 + m log | Σ| + 2 N (p + 1) N -p -2 , [ 18
]
is an approximate unbiased estimator of the first term in Ω(θ 0 ) which is in our case the right hand side of equation ( 16).

Let us remark that

β = (X t V -1 X) -1 X t V -1 (X 0 β 0 + ε 0 ) = (X t V -1 X) -1 X t V -1 Xβ * + (X t V -1 X) -1 X t V -1 ε 0 = β * + (X t V -1 X) -1 X t V -1 ε 0 . [19]
Using the assumption that φ is a consistent estimator of φ 0 and the fact that β -β * = (X t V -1 X) -1 X t V -1 0 , we compute the error term given in equation [START_REF] Jones | Longitudinal Data With Serial Correlation: a state-space approach[END_REF]. In the first term, we have the approximation

IE θ0 m log |Σ 0 | | Σ| ≈ 0. [20]
On the other hand, since V = V 0 + o p (1) then the third term in the right-hand side of ( 17) becomes

IE θ0 σ2 σ 2 0 tr( V V -1 0 ) ≈ IE θ0 t 0 ( V -1 -V -1 X(X t V -1 X) -1 X t V -1 ) 0 σ 2 0 = IE θ0 t 0 (V -1 0 -V -1 0 X(X t V -1 0 X) -1 X t V -1 0 ) 0 σ 2 0 = N -p.
[21]

Furthermore, using the expression [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF] and by a simple matrix calculus, the fourth term on the right-hand side of (17) becomes,

IE θ0 ( β -β * ) t X t V -1 0 X( β -β * ) σ 2 0 = IE θ0 t 0 V -1 0 X(X t V -1 0 X) -1 X t V -1 0 0 σ 2 0 = tr(X(X t V -1 0 X) -1 X t V -1 0 ), = p.
[22]

Now it remains to compute

IE θ0 N log σ 2 0 σ2 . It is known that the distribution of N ( σ2 σ 2 0 )
is a chi-square with (N -p) degrees of freedom. From lemma (3) in [START_REF] Cavanaugh | Unifying the derivations of the Akaike and corrected Akaike information criteria[END_REF], we have

IE θ0 N log σ 2 0 σ2 = (p + 1) + o(1)
.

[23]

According to [START_REF] Kullback | Information theory and statistics[END_REF] and substituting ( 20), ( 21), ( 22) and (23) in the right side of ( 17), we obtain the estimator of Ω(θ 0 ) given in the proposition.

REMARK. -It is important to note that if p held fixed and N tends to infinity, the second term in the expression [START_REF] Fitzmaurice | Applied Longitudinal Analysis[END_REF] of KICc tends to 3(p + 1). This yields the criterion KIC given in [START_REF] Diggle | Analysis of Longitudinal Data[END_REF].

Derivation of RIC sd criterion based on residual likelihood and symmetric divergence

The method of restricted (or residual) maximum likelihood is a way of estimating variance components in a general linear model. It was first introduced by Patterson and Thompson [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF]. By adopting the results of Diggle [START_REF] Diggle | Analysis of Longitudinal Data[END_REF] or Verbyla [START_REF] Verbyla | A conditional derivation of residual maximum likelihood[END_REF], and omitting irrelevant terms, the restricted likelihood for the candidate model ( 2) is defined by

log f (r) θ (Y ) = - 1 2 (N -p) log σ 2 + m log |Σ| + log |X t V -1 X| + Y t AY σ 2 [24]
where

A = V -1 -V -1 X(X t V -1 X) -1 X t V -1 .
The residual maximum likelihood estimator of (β, σ 2 , φ) are given by the equalities β = (X t Ṽ -1 X) -1 X t Ṽ -1 Y and σ2 = (Y -Xβ) t Ṽ -1 (Y -Xβ)/(N -p). Following Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF], Ṽ -1 is V evaluated at φ which is obtained in this setting by maximizing the residual likelihood given by r (φ) = (φ) -1 2 log |X t V -1 (φ)X|. Similarly, the restricted likelihood function for the true model [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF], can be obtained by replacing σ, β, Σ, X, p and V in equation ( 24) with σ 0 , β 0 , Σ 0 , X 0 , p 0 and V 0 respectively. A useful measure of the discrepancy between the residual log-likelihood function of candidate and true models is the Kullback's symmetric divergence defined by

Ω r (θ 0 ) = IE θ0 {d r (θ 0 , θ) + d r ( θ, θ 0 ) -d r ( θ, θ)}.
[25]

where

d r (θ 0 , θ) = IE θ0 {-2 log f (r) θ (Y )}.
Under the assumption that the true model is include in the family of candidate models, we show the following proposition: Proposition 2 The criterion defined by

RIC sd = (N -p) log σ2 + m log | Σ| + p log(N ) + (N -p) 2 N -p -2 +(N -p) log( N -p 2 ) -Ψ( N -p 2 ) [26]
is an approximate unbiased estimator of the discrepancy Ω r (θ 0 ), where Ψ denote the digamma function.

Proof: We have

Ω r (θ 0 ) = IE θ0 {d r (θ 0 , θ) + d r ( θ, θ 0 ) -d r ( θ, θ)}. [27]
From the model consideration, it is easy to see that X t A = 0. This implies the following:

d r (θ 0 , θ) = IE θ0 {-2 log f (r) θ (Y )} θ= θ = IE θ0 (N -p) log σ 2 + m log |Σ| + log |X t V -1 X| + (Y -Xβ * ) t A(Y -Xβ * ) σ 2 θ= θ = IE θ0 (N -p) log σ 2 + m log |Σ| + log |X t V -1 X| + (Y -X 0 β 0 ) t A(Y -X 0 β 0 ) σ 2 θ= θ = IE θ0 (N -p) log σ 2 + m log |Σ| + log |X t V -1 X| + ε t 0 Aε 0 σ 2 θ= θ = (N -p) log σ2 + m log | Σ| + log |X t Ṽ -1 X| + σ 2 0 σ2 tr( Ã-1 V 0 ), [28]
where Ã-1 is A evaluated at φ. Similarly, we have

d r ( θ, θ 0 ) = IE θ {-2 log f (r) θ0 (Y )} θ= θ = IE θ (N -p 0 ) log σ 2 0 + m log |Σ 0 | + log |X t V -1 0 X| + Y t A 0 Y σ 2 0 θ= θ .
Furthermore, we derive the following simple expectation equalities

IE θ Y t A 0 Y σ 2 0 θ= θ = IE θ (Y -Xβ * ) t A 0 (Y -Xβ * ) σ 2 0 θ= θ = σ2 σ 2 0 tr( Ṽ A 0 ).
With the same technique as in [START_REF] Hedeker | Applied Longitudinal Data Analysis[END_REF], we obtain

d r ( θ, θ 0 ) = (N -p) log σ 2 0 + m log |Σ 0 | + log |X t V -1 0 X| + σ2 σ 2 0 tr( Ṽ A 0 ) + ( β -β * ) t X t A 0 X( β -β * ) σ 2 0 .
[29]

As it was argued bellow, since X t A = 0 then we have the following

IE θ Y t AY σ 2 θ= θ = IE θ (Y -Xβ) t A(Y -Xβ) σ 2 θ= θ = IE θ ε t Aε σ 2 θ= θ = tr( Ṽ A).
The third term of ( 27) is then given by:

d r ( θ, θ) = IE θ {-2 log f (r) θ (Y )} θ= θ , = IE θ (N -p) log σ 2 + m log |Σ| + log |X t V -1 X| + Y t AY /σ 2 θ= θ = (N -p) log σ2 + m log | Σ| + log |X t Ṽ -1 X| + tr( Ṽ Ã).
[30] From (28), ( 29) and (30), we obtain,

Ω r (θ 0 ) = IE θ0 (N -p) log σ2 + m log | Σ| + log |X t Ṽ -1 X| + σ 2 0 σ 2 tr( Ã-1 V 0 ) [31] +IE θ0 (N -p) log σ 2 0 σ2 + m log |Σ 0 | | Σ| + log |X t V -1 0 X| |X t Ṽ -1 X| -tr( Ṽ Ã) + σ2 σ 2 0 tr( Ṽ A 0 ) + ( β -β * ) t X t A 0 X( β -β * ) σ 2 0 -p 0 log σ 2 0 . [32]
By adopting the same reasoning as Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF], we can approximate equation (31) by

N log σ2 + m log | Σ| + p log(N ) + (N -p) 2 N -p -2 .
[33]

By using the equality [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF] and simple matrix simplifications, we have

∆ IE θ0 ( β -β * ) t X t A 0 X( β -β * ) σ 2 0 = IE θ0 ε t 0 V -1 0 X t (X t V -1 0 X) -1 X t A 0 X(X t V -1 0 X) -1 X t V -1 0 ε 0 σ 2 0 = tr(X(X t V -1 0 X) -1 X t A 0 ) = tr(X(X t V -1 0 X) -1 X t V -1 0 ) -tr(X(X t V -1 0 X) -1 X t V -1 0 X 0 (X t 0 V -1 0 X 0 ) -1 X t 0 V -1 0 ). = p -p 0 .
[34]

Under the same assumption used in the previous section, φ is a consistent estimator of φ 0 , we can approximate (in probability) Ṽ by V 0 and Σ by Σ 0 . Hence, we have

IE θ0 m log |Σ 0 | | Σ| ≈ 0 and IE θ0 log |X t V -1 0 X| -log |X t Ṽ -1 X| ≈ 0.
On the other hand,

IE θ0 σ2 σ 2 0 tr( Ṽ A 0 ) = IE θ0 t 0 ( Ṽ -1 -Ṽ -1 X(X t Ṽ -1 X) -1 X t Ṽ -1 ) 0 (N -p)σ 2 0 × tr{(V -1 0 -V -1 0 X(X t V -1 0 X) -1 X t V -1 0 ) Ṽ } ≈ IE θ0 t 0 (V -1 0 -V -1 0 X(X t V -1 0 X) -1 X t V -1 0 ) 0 (N -p)σ 2 0 × tr{(V -1 0 -V -1 0 X(X t V -1 0 X) -1 X t V -1 0 )V 0 } = IE θ0 t 0 (V -1 0 -V -1 0 X(X t V -1 0 X) -1 X t V -1 0 ) 0 σ 2 0 . = N -p.
[35] Furthermore, we have

IE θ0 tr( Ã Ṽ ) ≈ IE θ0 tr(V -1 0 -V -1 0 X(X t V -1 0 X) -1 X t V -1 0 )V 0 ) = N -p.
[36]

Now it remains to compute IE θ0 (N -p) log

σ 2 0 σ2 . Since, the distribution of (N -p)( σ2 σ 2 0
) is a chi-square with (N -p) degrees of freedom, then a simple calculus shows that:

IE θ0 (N -p) log σ2 σ 2 0 = (N -p){log 2 N -p + Ψ( N -p 2 )} [37]
ignoring the constant term in (32) (p 0 log σ 2 0 ) and according to (33) and substituting (37), ( 35) and (36) in the right side of (32), we obtain the estimator of Ω r (θ 0 ) given in the proposition.

Simulation and conclusions

The purpose of this section is to study and compare the performance of both criteria KICc and RIC sd , introduced in this paper, to those studied in Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF], namely RIC and AICc, and the well known criteria AIC, BIC (Akaike [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF], Schwarz [START_REF] Schwarz | Estimating the dimension of a model[END_REF]). The performances are evaluated as a function of the number of subjects, the repeated measurements, the SNR and the correlation structure of longitudinal models. We generate 1000 realizations from model [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF], with p 0 = 3, β 0 = (1, 2, 3) t and variables X 0 is a three column explanatory matrix. The explanatory variables of the candidate model were stored in a p = 7 column matrix X; where the first three columns of X are those of X 0 . The explanatory variables were randomly generated from the standard normal distribution. The number of repeated measurement was fixed at n = 10 and the number of subjects m have been changed to take values m = 1, m = 5, m = 10 and m = 30. The SNR ratio was also controlled at 1, 5, and 10, where SNR = var(x 0 t ij β0) var( ij ) . We have considered the uniform correlation structure of longitudinal data with ρ tacking 0.5 and 0.9. The results are summarized in Table 1 in which we give the percentages of the correct model order selected by each criterion considered here. We recall that:

BIC = N log σ2 + m log | Σ| + p log(N ) AIC = N log σ2 + m log | Σ| + 2(p + 1) AICc = N log σ2 + m log | Σ| + 2 N (p + 1) N -p -2 RIC = N log σ2 + m log | Σ| + p log(N ) + (N -p) 2 N -p -2 .
From Table 1, we see that when SNR = 1 and m = 1 the KICc criterion performs better than the others criteria followed by AICc. However, RIC sd is superior than RIC, BIC, AIC and KIC. If the number of the subjects increases, RIC sd outperforms all criteria, except when SNR = 1 and ρ = 0.5 where BIC is the best. We note also, that in all setting KICc is better than AICc and when m is large (m = 30), BIC and RIC give approximatively the same results. This latter remark is mentioned also by Azari et al. [START_REF] Azari | Longitudinal data model selection[END_REF]. Moreover, when SNR increases both RIC and RIC sd are improved and the performance of these criteria is better as ρ increases. On the other hand the performance of KICc and AICc decays as m increases.

Other examples of simulations not reported here, with the autoregressive correlation structure of longitudinal data, give the same results as in the preceding example.

In conclusion, when both the number of subjects and the repeated observations are small, it is preferable to use the KICc criterion. However, when m or n is large one should use RIC sd or BIC. This latter is favored for a small SNR ratio and RIC sd for a moderate to large SNR ratio.

Conclusion

In this paper, we have derived two model selection criteria, KICc and RIC sd , for application in longitudinal data analysis. Our criteria are based on the Kullbak's symmetric divergence. A small simulation study is undertaken to compare the performance of our criteria to other well known criteria. Moreover, our simulation studies show that the KICc criterion outperforms all other criteria when the sample size and the SNR ration are small. RIC sd is favored when the sample size and the SNR ration are moderate to large. Furthermore, KICc and RIC sd are superior to AIC and RIC, respectively, which are based on the asymmetric divergence.

Table 1 :

 1 Percentages of correct model order selection for the uniform correlation structure

			m = 1	m = 5 m = 10 m = 30
			n = 10 n = 10 n = 10	n = 10
		AIC	20.9	66.7	69.7	70.5
		AICc	95	79.8	74.4	71.8
	ρ = 0.5	KIC	34.2	85.2	86.2	86.8
	SN R = 1	KICc	96.4	89.7	88.1	87.2
		BIC	32.4	93.4	96.1	98.4
		RIC	35.4	64.3	82.8	92.8
		RICsd	45.8	83.8	92.6	96.1
		AIC	21.1	67.3	69.7	72.7
		AICc	95.1	77.7	76.3	74.3
	ρ = 0.5	KIC	33.9	84.8	87.3	88.2
	SN R = 5	KICc	96.6	90.4	89.4	88.9
		BIC	32.1	89.2	90.0	94.9
		RIC	65.1	92.6	96.5	98
		RICsd	77.9	95.2	98.1	99
		AIC	20.0	64.8	69.9	74.2
		AICc	95.7	76.2	75.8	76.1
	ρ = 0.5	KIC	33.2	82.9	85.2	87.8
	SN R = 10	KICc	97.3	89.7	88.1	88.0
		BIC	32.0	92.4	96.2	99.0
		RIC	82	95.8	96.0	97.1
		RICsd	91	97.8	98.5	98.9
		AIC	24.9	66.4	69	72
		AICc	96.2	75.8	75.4	74.6
	ρ = 0.9	KIC	39.6	83.4	87.1	88.5
	SN R = 1	KICc	97.1	89.2	88.2	87.1
		BIC	34.2	92.6	95.6	98.7
		RIC	50.0	85.3	91.7	98.9
		RICsd	60.6	94.4	95.7	99.7
		AIC	21.2	67.7	70.7	75.5
		AICc	95.7	78	75.7	73.3
	ρ = 0.9	KIC	39.2	84.4	86.7	87.5
	SN R = 5	KICc	97.0	89.	88.7	88.1
		BIC	33.4	93.2	95.8	99.0
		RIC	89.4	95.7	98.5	99.0
		RICsd	94.6	98.1	99.1	99.9
		AIC	22.9	65.8	74.4	74.4
		AICc	97.2	77.5	79.3	75.2
	ρ = 0.9	KIC	39.7	83.9	87.3	88.2
	SN R = 10	KICc	97.8	89.6	90.8	88.7
		BIC	32.9	93.5	97.5	99.2
		RIC	95.1	98	99.1	99.8
		RICsd	97.6	98.6	99.9	100

Acknowledgements

This article was partially written while the first author was in post-doctoral at Mathematical Institute of Burgundy. He wishes to thank Veronique Maume-Deschamps for her helpful discussions. He also wants to thank the regional council of Burgundy for his research grant. The authors would like to thank the Associate Editor and the anonymous referee for their comments which helped improving the presentation of this paper.


Jean-Antoine Désidéri

MGDA Variants for Multi-Objective Optimization

Keywords: multi-objective optimization, descent direction, convex hull, Gram-Schmidt orthogonalization process, BFGS quasi-Newton method

, in comparison with an evolutionary strategy. Recently, a variant, MGDA-II, has been proposed in which the descent direction is calculated by a direct procedure [4] based on a Gram-Schmidt orthogonalization process (GSP) with special normalization. This algorithm was tested in the context of a simulation by domain partitioning, as a technique to match the different interface components concurrently [3]. The experimentation revealed the importance of scaling, and a slightly modified normalization procedure was proposed ("MGDA-IIb"). In this new report, two novel variants are proposed. The first, MGDA-III, realizes two enhancements. Firstly, the GSP is conducted incompletely whenever a test reveals that the current estimate of the direction of search is adequate also w.r.t. the gradients not yet taken into account; this improvement simplifies the identification of the search direction when the gradients point roughly in the same direction, and makes the Fréchet derivative common to several objective-functions larger. Secondly, the order in which the different gradients are considered in the GSP is defined in a unique way devised to favor an incomplete GSP. In the second variant, MGDA-IV, the question of scaling is addressed when the Hessians are known. A variant is also proposed in which the Hessians are estimated by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) formula.

Variantes de l'algorithme MGDA pour l'optimisation multiobjectif

Résumé : Ce rapport fait suite à plusieurs publications dans lesquelles on a proposé et testé un Algorithme de Descente à Gradients Multiples (MGDA) pour traiter les problèmes d'optimisation différentiable multiobjectifs. La méthode a été introduite originellement dans [2], et à nouveau formalisée dans [START_REF]Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF]. Sa capacité à identifier le front de Pareto a été mise en évidence dans [START_REF] Zerbinati | Comparison between MGDA and PAES for Multi-Objective Optimization[END_REF], en comparaison à une stratégie évolutionnaire. Récemment, une variante, MGDA II, a été proposée dans laquelle la direction de descente est calculée par une procédure directe [START_REF] Mgda Ii | A direct method for calculating a descent direction common to several criteria[END_REF] s'appuyant sur le processus d'orthogonalisation de Gram-Schmidt (GSP), effectué avec une certaine normalisation. On a testé l'efficacité de l'algorithme dans le contexte d'une simulation par partitionnement de domaine, comme technique pour raccorder concouramment les différentes composantes d'interface [3]. On a observé l'importance des facteurs d'échelle, ce qui a conduit a une légère modification de la procédure de normalisation ("MGDA-IIb"). Dans ce nouveau rapport, deux nouvelles variantes sont proposées. La première, MGDA-III, réalise deux améliorations. Premièrement, le GSP est exécuté incomplètement dès lors qu'un test révèle que l'estimation courante de la direction de recherche convient aussi aux vecteurs gradients qui n'ont pas encore été pris en considération; cette amélioration simplifie l'identification de la direction de recherche lorsque les gradients ont une direction commune approchée, et présente également l'avantage d'augmenter la dérivée de Fréchet commune aux fonctions objectifs. Deuxièmement, l'ordre dans lequel les differents gradients sont pris en compte dans le GSP est défini d'une manière unique conçue pour favoriser l'interruption rapide du GSP. Dans la seconde variante, MGDA-IV, la question du choix des échelles appropriées pour normaliser les gradients est abordée en supposant les hessiens connus. Une variante est également proposée dans laquelle les hessiens sont estimés par la formule de Broyden-Fletcher-Goldfarb-Shanno (BFGS).

Mots-clés : optimisation multiobjectif, direction de descente, enveloppe convexe, processus d'orthogonalisation de Gram-Schmidt, méthode BFGS de quasi-Newton

Introduction

We consider the simultaneous minimization or reduction of n objective-functions, {J i (Y )} (i = 1, . . . , n), assumed to be smooth (say C 2 ) functions of the design-vector Y ∈ R N , where n ≤ N .

Our analysis is developed to identify an appropriate direction of search ω to update the design vector from a given design-point Y 0 , center of an open ball B in which the objective-functions are well-defined and smooth:

Y 1 = Y 0 -ρω (ρ > 0, step-size) (1)
The above iteration is a descent step, if two conditions are met. Firstly, the vector ω ∈ R N should be such that ∀i = 1, . . . , n

: ∇J i (Y 0 ), ω > 0 . (2)
Then, -ω is a descent direction common to all objective-functions. Secondly, the step-size ρ should be adjusted appropriately. In this report, we only focus on the first condition.

In [2] and [START_REF]Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF], we have introduced the notion of "Pareto-stationarity":

Definition 1 (Pareto-stationarity)
The design-point Y 0 is said to be Pareto-stationary if there exists a convex combination of the gradients, ∇J i (Y 0 ), that is equal to 0:

∃α = {α i } (i = 1, . . . , n) such that : α i ≥ 0 (∀i); n i=1 α i = 1; n i=1 α i ∇J i (Y 0 ) = 0 . (3)
We have shown that Pareto-stationarity is a necessary condition to Pareto-optimality. Thereafter, we restrict ourselves to the case where the initial design-point Y 0 is not Pareto-optimal, thus not Pareto-stationary.

Clearly, the above condition (2), as it only involves scalar products, can be applied to projected gradients, in case of constrained minimization. More specifically, suppose that the active scalar constraints at Y = Y 0 are the following:

g 1 (Y 0 ) = g 2 (Y 0 ) = • • • = g K (Y 0 ) = 0 . (4)
Then define the vectors

v k = ∇g k (Y 0 ) (k = 1, . . . , K) (5)
normal to the constraint surfaces. Apply the GSP to these vectors to get the family {w k } (k = 1, . . . , K) of orthonormal vectors that collectively span the same subspace, and define the following projection matrix:

P = I N - K k=1 [w k ] [w k ] T (6)
where the bracketed vector [w k ] stands for the column-vector of its components viewed as a N × 1 matrix, and the superscript T indicates transposition. Then, the forthcoming MGDA construction should be applied after replacing the original gradients, ∇J i (Y 0 ), by their projections onto the subspace tangent to the constraint surfaces, that is, P ∇J i (Y 0 ). This remark being made, without loss of generality, we are considering thereafter the unconstrained formulation.

In the original formulation of MGDA [2]- [START_REF]Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF], the vector ω is defined as the minimum-norm element in the convex hull of the gradients:

U = u ∈ R N / u = n i=1 α i ∇J i (Y 0 ) ; α i ≥ 0 (∀i); n i=1 α i = 1 (7)
that is:

ω = Argmin u∈U u ( 8
)
This definition is the most general and in particular, it is applicable whether the gradient vectors are linearly independent or not. The element ω can be identified by numerical minimization in the convex hull, which can be parameterized isomorphically to the hypercube [0, 1] n-1 (see [START_REF]Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF]). This minimization can however be numerically delicate, and in fact, not necessary.

The convex hull can also be viewed as an affine structure, since:

∀u ∈ U : u -u n = n i=1 α i u i - n i=1 α i u n = n-1 i=1 α i u n,i (u n,i = u i -u n ) (9)
Hence, U ⊆ A n-1 (or using affine-space notations, U ⊆ Ȧn-1 ), where A n-1 is a set of vectors pointing to an affine sub-space Ȧn-1 of dimension at most n -1.

Let us examine these affine and vector structures, with the support of U (affine)

U (vector) Ȧ2 O O ⊥ u 1 u 2 u 3 u ∈ U
u i , ω = ω 2 (∀i) (10)
This is a highly favorable situation in which all the Fréchet derivatives are equal.

Inria

The element ω being defined, the MGDA iteration generalizes the classical steepest-descent method [START_REF]WikipediA: The Free Encyclopedia[END_REF] to multi-objective optimization by using the vector -ω as the direction of search. We have established that, certain weak provisions on the problem formulation being made, if the step-size is adjusted optimally, the iteration accumulates at a Pareto-stationary design-point. Whenever ω = 0 (or, in a numerical procedure, when ω < T OL), the current design-point is Pareto-stationary, and the optimization is interrupted. Hence, in [START_REF] Zerbinati | Comparison between MGDA and PAES for Multi-Objective Optimization[END_REF], the efficacy of MGDA to identify the Pareto front has been demonstrated, and comparisons with an evolutionary strategy (PAES) have been made.

Recently, a variant, MGDA-II, has been proposed in which the descent direction is calculated by a direct procedure [START_REF] Mgda Ii | A direct method for calculating a descent direction common to several criteria[END_REF] based on a Gram-Schmidt orthogonalization process (GSP) with special normalization. This version is valid only when the gradient vectors are linearly independent, and it proceeds as follows. Due to the importance of scaling, user-supplied scaling factors {S i } (i = 1, . . . , n), are assumed given, and the following scaled gradients are defined:

J i = ∇J i (Y 0 ) S i (11)
(S i > 0; e.g. S i = J i for logarithmic gradients). The GSP is performed as follows:

• Set u 1 = J 1
• For i = 2, . . . , n, set:

u i = J i -k<i c i,k u k A i (12)
where:

c i,k = J i , u k u k , u k (∀k < i) (13)
and:

A i =      1 - k<i c i,k if nonzero ε i otherwise (ε i arbitrary, but small). ( 14
)
As a result of this construction, a new element ω is defined, as the minimum-norm element in the convex hull of the orthogonal vectors

{u i } (i = 1, . . . , n) ω = n i=1 α i u i (15)
in which the coefficients {α i } are strictly positive and less than 1:

α i = 1 u i 2 n j=1 1 uj 2 = 1 1 + j =i ui 2 uj 2 < 1 . ( 16
)
Due to the orthogonality of the family {u i } (i = 1, . . . , n), ω = ---→ OO ⊥ , and (10) holds. This is illustrated by Fig. 2. Consequently

J i , ω = A i + k<i c i,k ω 2 = ω 2 (17)
by definition of the normalization constant A i . In conclusion, the Fréchet derivatives of all objective-functions are equal.

It should be emphasized that in general the newly-defined element ω is distinct from the former, except in the particular case of two objective-functions (n = 2), when the gradient vectors form an obtuse angle. It was later observed that situations in which a normalization constant A i for some i was negative was to be avoided. To do this, in MGDA-II b, the definition

u 1 u 2 O u 3 O ⊥ u ∈ U
A i = 1 - k<i c i,k (18)
was maintained only if this number is strictly-positive. Otherwise, we have used the following modified scale:

S i = k<i c i,k S i (19)
so that:

c i,k = k<i c i,k -1 c i,k (20)
and

k<i c i,k = 1 , (21)
and A i = ε i , for some small ε i . (This procedure was referred to as the "automatic rescale procedure".)

This has led us to the same formal conclusion: the Fréchet derivatives are equal; but the value is much larger, and (at least) one objective-function has been rescaled. This variant was tested on a somewhat peculiar model problem of domain-partitioning, in which all objective-functions tend to 0, which results in a Pareto set restricted to a single point. In this application rather remote from the context for which MGDA has been devised, using logarithmic scaling of the gradients (S i = J i ) and automatic rescale, MGDA-II b was found to converge satisfactorily; in fact, at a rate only twice less than the quasi-Newton method; additionally, the iteration indicated asymptotically an interesting trend to convergence acceleration [3].

Nevertheless, these developments have brought up some open questions to which this report brings certain answers. In particular, the following three:

1. Since the element ω provided by MGDA-II is in general different from the original one, how can we guarantee the convergence of MGDA to a Pareto-stationary design-point?

Inria 2.
In which order should we arrange the gradients to perform the GSP ?

3. Can we devise an adequate scaling of the gradient inspired from quasi-Newton methods?

MGDA-III provides an answer to the first two questions, and MGDA-IV to the fourth.

2 Ordered and economical GSP : MGDA-III

The driving idea is that in case of numerous gradients, trends might emerge among them, permitting to account for the general direction of subgroup by a unique vector in the orthogonal basis. Hence the GSP could be interrupted as soon as a direction is found to be a descent direction common to all objective functions while being constructed on the basis of only I < n gradients.

To achieve this purpose, in the following algorithm, at the stage of computing a new orthogonal vector, the gradient that is elected among those not yet accounted for, is the one for which the scalar product with the current estimate of the element ω is algebraically smallest. In this way, it is the vector for which the construction so far is the least satisfactory. Thus, computational economy is achieved through the specification of the ordering in which the gradients are used to perform the GSP , with the expectation of a rapid interruption of the process. Further comments on the expected gain in efficiency will be made subsequently.

Algorithm

Again, one starts from scaled gradients, {J i }, and duplicates, {g i },

J i = ∇J i (Y 0 ) S i (S i : user-supplied scale) (22)
and proceeds in 3 steps: A,B and C.

A: Initialization

• Set 1 u 1 := g 1 = J k / k = Argmax i min j J j , J i J i , J i (23)
• Set n × n lower-triangular array c = {c i,j } (i ≥ j) to 0.2

• Set, conservatively, I := n. 3• Assign some appropriate value to a cut-off constant a : (0 ≤ a < 1).

B: Main GSP loop; for i = 2, 3, . . . , (at most) n, do:

1. Calculate the i -1st column of coefficients:

c j,i-1 = g j , u i-1 u i-1 , u i-1 (∀j = i, . . . , n) (24)
and update the cumulative row-sums:

c j,j := c j,j + c j,j-1 = k<i c j,k (∀j = i, . . . , n) (25)
2. Test:

-If the following condition is satisfied

c j,j > a (∀j = i, . . . , n) (26)
set I := i -1, and interrupt the GSP (go to 3.).

-Otherwise, compute next orthogonal vector u i as follows:

# Identify index = Argmin j {c j,j / i ≤ j ≤ n}. 4# Permute information associated with i and :

g-vectors: g i g rows i and of array c and corresponding cumulative row-sums

c i,i c , # Set A i = 1 -c i,i ≥ 1 -a > 0 and calculate 5 u i = g i -k<i c i,k u k A i (27)
# If u i = 0, return to 1. with incremented i; otherwise:

g i = k<i c i,k u k = k<i c i,k g k (28)
where the {c i,k } are calculated by backward substitution. Then, if c i,k ≤ 0 (∀k < i):

Pareto-stationarity detected: STOP MGDA iteration; otherwise (exceptional ambiguous case):

STOP GSP ; compute ω according to original definition and go to C.

3.

Calculate ω as the minimum-norm element in the convex hull of {u 1 , u 2 , . . . , u I }:

6 ω = I i=1 α i u i = 0 (29)
where:

α i = 1 u i 2 I j=1 1 uj 2 1 1 + j =i ui 2 uj 2 . ( 30
)
C: If ω < T OL, STOP MGDA iteration; otherwise, perform descent step and return to B.

Properties

Case I = n. In this case, the GSP is performed completely, and the algorithm is equivalent to the former MGDA-II with the enhancement that the rescale of the b-version is never necessary, since the specified ordering implies that ∀i :

A i ≥ 1 -a > 0 . (31)
Case I < n (incomplete GSP ). Here, the Fréchet derivatives satisfy different bounds according to two subcases:

• First I Fréchet derivatives:

g i , ω = u i , ω = ω 2 > 0 (∀i = 1, . . . , I) (32)
• Subsequent ones (i > I):

By construction, the vectors {u 1 , u 2 , . . . , u I } are orthogonal, and ω is given by (29), so that:

g i = I k=1 c i,k u k + v i (33)
where v i ⊥ {u 1 , u 2 , . . . , u I }. Consequently,

g i , ω = I k=1 c i,k u k , ω = I k=1 c i,k ω 2 = c i,i ω 2 > a ω 2 > 0 (34)
2.3 A posteriori justification of the choice of g 1

At initialization, we have set

u 1 := g 1 = J k / k = Argmax i min j J j , J i J i , J i (35)
This was equivalent to maximizing c 2,1 = c 2,2 , that is, maximizing the least cumulative row-sum, at first estimation. This makes at start, the worst case less severe. One anticipates that the favorable situation for which all cumulative row-sums are positive (or > a), is more likely.

Expected benefits

According to the section above, the specified ordering has been devised to permit the GSP to be performed incompletely. When gradients exhibit a general trend, ω is found in fewer steps and this realizes a computational economy.

Secondly, the rescale procedure has been circumvented.

Thirdly, an incomplete GSP results in an element ω of larger norm since it realizes the minimization in a smaller subset, namely the convex hull of an incomplete orthogonal basis. This corresponds to larger directional derivatives, since

g i , ω ω = ω or a ω (36)
and to the greater efficiency of the subsequent MGDA descent step.

3 Using Hessians to better scale the gradients: MGDA-IV

Adressing the question of scaling when Hessians are known

In single-objective optimization, when both gradient and Hessian are known, Newton's method is the most effective unless additional information is provided.

For the optimization of the objective J i (Y ) alone, Newton's method writes:

Y 1 = Y 0 -p i ( 37
)
where the vector p i is given by the solution of the system:

H i p i = ∇J i (Y 0 ) (38)
where H i is the Hessian matrix of objective function J i at Y = Y 0 . Hence the preconditioning by the inverse Hessian realizes a form of optimal scaling. However, in general, the vector p i is not parallel to the gradient itself. Thus to ensure that the iteration remains a descent step, only its projection should be retained.

Thus, we propose to split the vector p i into orthogonal components

p i = q i + r i (39)
where:

q i = p i , ∇J i (Y 0 ) ∇J i (Y 0 ) 2 ∇J i (Y 0 ) (40)
and r i ⊥ ∇J i (Y 0 ), and to define the scaled gradient as follows:

J i = q i (41)
We define MGDA-IV as MGDA-III applied to the gradients scaled as above. With the previous notations, this is equivalent to defining the scaling constant S i as follows:

S i = ∇J i (Y 0 ) 2 p i , ∇J i (Y 0 ) (42)

BFGS-inspired variant: MGDA-IV b

When the Hessians are not known exactly, they can be approximated by the BFGS iterative estimate:

(∀i = 1, . . . , n) H (0) i = Id (43) H (k+1) i = H (k) i - 1
s (k) T H (k) i s (k) H (k) i s (k) s (k) T H (k) T i + 1 z (k) T i s (k) z (k) i z (k) T i (44)
in which k is the MGDA iteration index, and:

s (k) = Y (k+1) -Y (k) (45) z (k) i = ∇J i (Y (k+1) ) -∇J i (Y (k) ) (46)

Conclusion

The different proposed variants of the MGDA are summarized in Table 1 with the indication of the major enhancements realized by each version, and references to publications and test-cases considered.

At this stage, further research is necessary to experiment MGDA-III and MGDA-IV and assess in practical engineering applications the improvements actually achieved by the theoretical enhancements. In particular, test-cases involving more complex Pareto fronts corresponding to a larger number n of objective functions, and possibly discontinuities, should be examined. Scaling with exact or approximate Hessians is a very promising option, but it should be verified cautiously also. Additionally, the question of step-size adjustment should also be investigated thoroughly.

Fig. 1

 1 drawn in the case n = 3. Here vectors are represented in the R 3 affine space with a given origin O. The gradient vectors are here denoted {u i } (i = 1, 2, 3). The convex hull of the gradients is the set of vectors of origin O pointing onto the triangle made of the 3 endpoints of {u i }. This triangle lies in a plane (generally speaking a subspace of dimension at most n -1) denoted Ȧ2 . The orthogonal projection of O onto the plane Ȧ2 is denoted O ⊥ . The figure has been drawn in the case where O ⊥ / ∈ U.

Figure 1 :

 1 Figure 1: Affine and vector structures: here, three vectors {u i } (i = 1, 2, 3) are considered to define the convex hull U, and u is an arbitrary element in U; all four vectors are associated with representatives of origin O; the endpoints of {u i } define the green triangle, U, affine structure associated with the convex hull U; U is contained in the plane Ȧ2 ; O ⊥ is the orthogonal projection of O onto Ȧ2 ; the figure illustrates the case where O ⊥ / ∈ Ȧ2 . Now, in the particular case where O ⊥ ∈ U, or equivalently, ---→ OO ⊥ ∈ U, and since ---→ OO ⊥ ⊥ U, ω = ---→ OO ⊥ , and by orthogonality:

Figure 2 :

 2 Figure 2: Construction of orthogonal vectors in MGDA-II

  principle related to minimum-norm element in convex hull of gradients [2] [6] -CV proof to Pareto stationary designs -Insensitive to Pareto front convexity -Multiple quadratics -Fonseca testcase (non-convex Pareto front; comparison MGDA vs PAES) [7] -DDM for Poisson pb. [3] [5] b-version -Meta-model assisted gradient computation [8] -CV requires a few database enrichments -Eulerian flow about wing [8] -Navier-Stokes duct flow [8] MGDA-II -Direct computation of descent direction by GSP [4] [3] -Requires linearly independent gradients -Modified definition of descent direction, and Pareto-stationarity test necessary n! possible orderings -DDM for Poisson pb. (scaling essential; verified CV to unique Paretostationary solution) [3] b-version -Automatic gradient rescale when normalization coefficient is found <0 -More efficient (larger directional derivatives) -(id.) MGDA-III -Specific ordering in GSP -Incomplete GSP -Resort to standard MGDA when Paretostationarity test ambiguous -Not limited to linearly-independent gradients -Even larger directional derivatives -Pareto-stationary accumulation points MGDA-IV -Scaling inspired from Newton's method using Hessians -Step-size of order unity expected to be appropriate b-version -Uses BFGS approximations to Hessians

Table 1 :

 1 Variants of MGDA with details on progressive enhancements

RR n°8068

The choice made for g 1 will be justified afterwards.

The main diagonal of array c is to contain cumulative row-sums.

The integer I ≤ n is the expected number of computed orthogonal basis vectors.RR n°8068

Note that necessarily c , ≤ a < 1.

c i,i = former-c , ≤ a; g i = former-c ,k .

Note that ω is calculated on the basis of a smaller number of gradients if I < n; here all computed u i = 0, and 0 < α i < 1.Inria

Inria

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


MRS and electrical prospection in the context of weathered peridotite rocks in the South of New Caledonia.

Girard 1* J-F., J-M. Baltassat 1 , A. Legchenko 2* A., S. Morlighem 3 , I. Domergue-Schmidt 3 , P. Maurizot 4 , J-L. Folio 3 and B. François 1 1 BRGM, Orléans, France; 2 IRD / LTHE, Grenoble, France; 3 VALE, New Calédonia. 4 BRGM, Nouméa, New Caledonia; jf.girard@brgm.fr

In November 2011 and July 2012 were conducted geophysical prospections in the southern part of New Caledonia. Both electrical tomography (ERT) and Magnetic Resonance Sounding (MRS) were performed to obtain a better understanding of the water storage and circulation. The geological context is the weathering profile of peridotite rock [START_REF] Maurizot | Geological map of New Caledonia at 1/500 000[END_REF].

In such environment, the general distribution of geology is (from surface to depth), the iron cap, laterite (red and yellow), saprolite, fractured and then fresh bedrock. The lateritic profile is developed over a thickness ranging from 20m to 60 m.

Electrical imagery is a well-established method to detect variations in depth and laterally in this geological context [START_REF] Robineau | Geoelectrical imaging of a thick regolith developed on ultramafic rocks : groundwater influence[END_REF] and is routinely used in mining exploration and exploitation to estimate the thickness of laterite.

Hydrogeology in this environment is complex because water flows in heterogeneous media: infiltration trough the iron cap, unsaturated zone, fractured bedrock. As a consequence, various regimes of hydrological responses are observed, from low permeability (laterite) to high permeable zones (fractured zone and locally "pseudo-karst" behavior).

MRS as a tool to detect directly water and provide insight of higher permeable zone sounds attractive in such a context. The sites studied are well documented thanks to preexisting boreholes logs and hydrogeological studies.

Difficulty to perform MRS there is link to two major reasons. First, the rocks present nonnegligible magnetic susceptibility. Despite surface measurement of the geomagnetic field revealed to be relatively homogeneous at the loop scale (< 100 nT variation), standard Free Induction Decay (FID) MRS measurement appeared to be un-practicable (like observed by [START_REF] Roy | Widespread occurrence of aquifers currently undetectable with the MRS technique in the Grenville geological province[END_REF]. The average magnetic susceptibility is 5 10 -4 SI and it proved to be suitable to perform MRS measurement in Spin Echo (SE) mode [START_REF] Roy | MRS in spineco mode -2008 tests in the Greenville[END_REF][START_REF] Legchenko | Application of the magnetic resonance sounding method to the investigation of aquifers in the presence of magnetic materials[END_REF][START_REF] Vouillamoz | Characterizing aquifers when using magnetic resonance sounding in a heterogeneous geomagnetic field[END_REF].

The second difficulty is link with the very fine structure of laterite, where a large part of water is not detectable by MRS like bound water in clay. But the high yield zones revealed to produce a clear MRS signal in SE mode.

We present a review of the various ERT and MRS responses observed in this context.


Reza Farahbakhsh

Noel Crespi

Ángel Cuevas

Sraddha Adhikari

Mehdi Mani
email: mehdi.mani@itron.com

Teerapat Sanguankotchakorn
email: teerapat@ait.ac.th

socP2P: P2P Content Discovery Enhancement by considering Social Networks Characteristics

Keywords: Content Discovery, Peer-to-Peer networks, Social Network, Random Walk

Content management appears as an essential requirement in order to deploy enhanced P2P networks. In P2P networks, the content that is requested by the query node could be located at different locations/nodes; therefore an efficient search mechanism is required. The proposed search algorithm in this paper, called socP2P, relies on peers' social relationships (friendships, shared interests, shared background and experiences) to improve the content discovery compared to similar solutions. With socP2P nodes can improve searches by using knowledge gained by 'overheard information' during their stay in the network. In addition, our algorithm exploits peers' common interests, friendships, and capability of memorizing experiences (received and routed queries) by them. Simulation results show that socP2P is able to achieve a high success rate, low delay and low overhead. We have verified that our algorithm is not only useful in finding popular contents in the network but also good enough to locate rare files. The obtained results reflect that exploiting social information in P2P networks leads to a more efficient content search mechanism.

I. BACKGROUND AND RELATED WORK

A Peer-to-Peer (P2P) network is a virtual overlay network, built at application layer, on top of a physical network topology. Decentralization and self-organization are key features of P2P networks. AS the number and location of contents in a P2P system have a strong impact on the performance, it is not surprising that content management in P2P networks is one of the research challenges in this area. Peers do not have global knowledge about the network and they have to interact with each other to find the required resource. With the increasing number of peers, the rate of possible interactions may also increase exponentially and hence this strategy is not applicable. In addition, existing unstructured P2P search techniques create huge network traffic or they require high search overhead [START_REF] Liu | research towards efficient social in the Global Information Technology[END_REF].

Most of the search methods for unstructured P2P networks are focusing on solving their issues by means of blind flooding. Authors in [START_REF] Liu | S resource discovery[END_REF] present a new social-like P2P algorithm (Social-P2P) for resource discovery by mimicking human interactions in social networks where peer nodes are people and connections are relationships. Random Walk [START_REF] Bisnik | Modeling and Analysis of Random Walk Search Algorithms in P2P Networks[END_REF] is one popular searching method in P2P which a query is forwarded to a randomly chosen neighbor hop by hop until enough responses have been found. Also Gnutella [START_REF] Ripeanu | Peer-to-Peer Architecture Cas Network[END_REF] appears as one of the most popular P2P file-sharing applications which require nodes to broadcast messages to their neighbors in order to perform content lookup.

Given that recent studies of P2P systems have provided a series of useful solutions to improve the performance of content management in P2P networks, realizing an advanced method of content management by means of social relations of users in OSNs is still an open issue in ongoing research. The objectives, requirements and challenges for P2P social networks have been formulated in [START_REF] Buchegger | A case for p2p infrastructure for social networks -opportunities & challenges[END_REF].

In this paper, we use interest similarity gestures (exhibited by the fact that both nodes are interested in the same content) to establish P2P links rather than relying on random associations or DHT (Distributed Hash Table ) rules. The proposed algorithm, named socP2P, uses social relation of peers to improve content management. Finally, we evaluate the performance of the proposed algorithm via simulation in terms of average success rate and path length to desired contents. In addition, we compare socP2P with Random walk algorithm.

The rest of this paper is organized as follows. The proposed resource search algorithm, socP2P, is presented in section II. Section III contains the performance evaluation of the proposed algorithm. Finally, we conclude the paper and present future direction in Section IV.

II. SOCP2P, PROPOSED RESOURCE SEARCH ALGORITHM

A. Proposal Details Overview

As stated earlier, the driving concept in this paper is to extract the social characteristics of the human beings behind a network of peers in a P2P system, and use them to improve the content management. We want to identify nodes with similar interests and enable them to help each other to find desired content. Nodes do not need to declare their interests, because inferring the interest of other nodes in the network is a gradual process for all nodes in the network. The more time they spend in the network the more they learn. Hence, they can route queries towards the right nodes more efficiently. In addition, peers may develop friendship with nodes having similar interest, so peers can generate their friends list and use it to route queries in the network.

In our design, we create new links between peers if we find/detect an interest similarity between them. For example, when node 'A' makes a request for object 'X', and if it gets a successful reply from node B, both nodes learn that they have a matching interest. Our approach differs in detecting interests between nodes compared to [START_REF] Upadrashta | to-Peer Systems[END_REF] that use the concept of providing a strength value to each node in a friend list. Furthermore, a node that generates a query and the node that successfully responds to that query are considered to share similar interest. Thus each interaction between nodes helps to determine if two nodes have interest similarity or not. In addition, each node has the capability of memorizing those received queries even if they are not successful. Therefore, it is also able to recommend a "right node" (i.e. a node that likely stores the searched content).

We have tried to use interest similarity gestures to make P2P links rather than relying on random associations or DHT rules. The links are created in two ways based on interest similarity. In the first method, Directed Interest Links (DIL) is created based on interest similarity between nodes. For simplicity, we set that after nodes A and B learn that they have a matching interest, node A creates an Interest Link with node B. We avoid creating a bidirectional link here, because as it happens in our normal life, relationship links may not always be mutual. For example as OSNs have shown us, some people are very well known and popular in the society, so many people like to create a link with them. However, these kinds of links are not often bidirectional.

The second type of link creation is the Network Suggested Forwards (NSF) and then Created Links model. The above association is the simplest and it is quite limited because it does not exploit the fact that there may be other copies of object X in the network, which node A and node B may not be aware of. The network is usually flooded with queries directed to other nodes for searching purposes. A node may receive queries for the same content more than once. In this case, the node that receives the query has already learned or memorized that some other nodes were looking for the same content earlier and instead of flooding the query in the network or forwarding the query in a random way, the query node can be recommended to get the content directly from the node that was looking for the same content earlier. Moreover, it is now possible to learn that the two nodes share interest profiles and so there is a possibility for them to create an interest-link which otherwise would have not been discovered. A network node keeps records and uses search queries to automatically adapt and make note of objects that indicate related interest relationships.

Fig. 1 shows an NSF model example in which node A is looking for content X. Node A sends query to node E, which memorizes the received query and since node E does not have the content X, it forwards the query to node C. Node C has the content X, thus it sends back the requested content to querying node A. Later in the time, node B looks for the same content X. Then, it sends its query to node E, which has previously encountered the same query. Thus, it recommends node A to node B assuming that node A might have the content X as it was looking for it earlier. Node E also forwards node B's query to node D. Since node D has the content X, it gives back the desired content to node B. The recommended node A also provides the content to node B. In this scenario, node A and node B are interested in same content X. Node C and node D already had the content X. Thus this gives us the clue that all four of these nodes (A, B, C and D) may be interested in the same contents. Thus, Directed Interest Links (DILs) can be created among these four nodes, as shown in Fig. 2.

B. Three Ways Approach to Search

In our proposed design, we establish a restriction that a query node can send its query directly or via intermediate nodes to a maximum of 'M' nodes. Our objective is to restrict the number of content queries injected in the network for a particular content to M. However, there are two main options for forwarding the queries: in parallel or sequentially. We decide to use a parallel approach since it leads to a better response time [START_REF] Crespo | systems[END_REF]. If the query node is unable to find the content in the network (with this restriction of M number of nodes) it will get it directly from the central server. A query node can start its network search in three different ways as shown in figure 3.

The initial method and preferred one is named Recommended Based Search (RBS). The first action of the query node is to check its query record and check if it has already encountered a query for the same content that it is looking for. If it finds that some other node was looking for the same content earlier, it will try to get the content from that node, which we call a Recommended Node. The second way is by making use of useful friends, which is called Friends-Based Search (FBS). Useful friends is a list of the friends of each node that have previously helped that node to find content, either by directly providing the desired content or by forwarding the query to the right node (meaning a node that had the requested content). The third method called Neighbor-Based Search (NBS). If a node doesn't have any recommended node for the looking content and also it does not have any useful friend(s), it will send the query to its neighbors. If the number of Neighbors is greater than M, it randomly chooses M neighbors.

C. Handling Unnecessary Traffic

In the example of Fig. 4, node S sends a query to nodes L, M and P. Nodes L and P are unaware that M has already received the query from node S. Therefore, they forward the query to M again. Similarly, neither L nor M, are aware that they have already received the query from node S, thus they both forward the query to each other. Therefore, M could receive the same message up to four times which we consider as unnecessary traffic in the network. It is clear that eliminating the logical links LM, MP, LQ and MQ does not reduce the search scope of node S. In our algorithm, when S forwards a query to M, it also appends information about the total list of nodes that the query has been forwarded to. Thus, when node L receives a query from node S, it also gets the information. Thus L does not forward the same query again. This technique has two advantages: (i) it avoids to forward queries again and again to the same node, thus eliminating possible loops and unnecessary traffic from the network; (ii) it increases searching area by forwarding queries to new nodes.

D. Discovering New Friends

P2P systems are dynamic, which means that peers join and leave the system and also their interest changes over the time. Given this dynamism, nodes always have a need to discover new friends. In our model, nodes add a new node in its friend list if this node successfully replies to a node's query. This update is realized because there is an interest similarity between the two nodes. However, there is no guarantee that a friend will be always useful to the node since (i) the interest of a node or its friend might change, and (ii) some content might be so popular that it is stored by a major portion of the nodes, thus a query for this content could be replied successfully by many nodes. In our model, a rank (or weight) is provided to each friend based on their usefulness. The rank of a friend is increased every time they prove to be useful. To maintain a fix maximum number of friends, the friends that become the least useful based on their ranking are removed from the list when a new useful friend is found.

III. PERFORMANCE EVALUATION

In this section we present the simulation methodology as well as the performance evaluation of socP2P algorithm. In addition, we also present a summary that shows how socP2P compares with the Random Walk algorithm and Interest Based Shortcuts [START_REF] Sripanidkulchai | using interest-based locality in peer 40[END_REF].

A. Simulation Methodology

We create a random topology by means of random bidirectional links between peers. Initially peers don't have information regarding to interest of others, thus they keep an empty Friend List. In addition, initially peers present an empty query record, since nodes have not encountered any queries.

The number of resources in the network is denoted by F. Initially, content/resources are distributed randomly in the network based on a popularity parameter of each resource. These popularity parameters are unique and randomly assigned at the beginning. The queries are generated sequentially by random nodes. Then, a query node randomly chooses one of the contents as its query resource (the content that query node is looking for). Each query includes: information about the query node and the query resource, a TTL value of M. Then a query is forwarded to M/2 nodes at first. These M/2 nodes could be useful friends or neighbors. These nodes could forward the query to one neighbor or Useful Friend of theirs. Thus in total, query could be forwarded to M (1 x M/2 + M/2 x 1) nodes. The Success Rate is the rate at which our algorithm resolves queries. In other words, it is how often a query node has found the content it is looking for in the network. The success rate is calculated using the following formula in our algorithm:

Success Rate = Success NBS + Success FBS + Success RBS TQ × 100%

B. Results and discussion

The plot in fig. 5 shows that ASR is increasing with increasing value of M. From this analysis, we can determine a suitable value of M according to our success rate requirements. As indicated in Fig. 6, if the number of resources (F) increases in the network, the ASR tends to decrease slightly. However, for higher values of M, the effect is not significant.

978-1-4673-2713-8/12/$31.00 ©2012 IEEE Fig. 6. ASR for N=100, Q=1000 and different values of F Fig. 7 shows the average success rate wit total number of valid queries in the network (Q).

C. Comparison of socP2P and Random Walk

Results in figure 8 show that socP2P has the highest success rate for all values of M. The improvement socP2P is especially significant for lower values of M ( 2%), verifying that socP2P can achieve a without flooding the network with query mes Finally, we also found that the average path length terms of overlay hops) across several simulated 3.5 hops in case of using Gnutella, whereas it is reduced only 1.27 hops away in case of Gnutella with shortcuts. However, socP2P further improves both results the path length to only 1.08 hops. Average Success Rates for N=100, F=1000 and Q=100, 1000, 10,000 of socP2P and Random Walk show that socP2P has the highest he improvement shown by ignificant for lower values of M (i.e. a high success rate without flooding the network with query messages. for N=500, F=100, Q=10,000 depicts the Average Success Rate versus the average path length (in ) across several simulated scenarios is 3.5 hops in case of using Gnutella, whereas it is reduced to only 1.27 hops away in case of Gnutella with shortcuts. However, socP2P further improves both results by reducing The approach also exploits from received and routed queries queries can be routed intelligently towards them. simulation results evaluate the performance of the proposed algorithm and demonstrate that socP2P well-known search mechanisms. rate than Random Walk and at the same time message overhead. Finally, socP2P not only the popular content, but also Many possible directions for future on this advance. Creating mutual interest links between nodes that have interest similarity, increasing the search success rates of less-popular files, include overhearing part in forwarded node query and also reducing the search range from two hops to one are some of our possible future endeavors

Fig. 1 .

 1 Fig. 1. Network Suggested Forwards and thus Created Links (NSF)

Fig. 3 .

 3 Fig. 3. Three Ways Approach to Search

Fig. 4 .

 4 Fig.4. An Example of P2P Overlay[START_REF] Xiao | Systems by Adaptive Connection Establishment on Computers[END_REF]

Fig. 5 .

 5 Fig. 5. ASR for N=100, F=100 and different values of Q

Fig. 7 .

 7 Fig. 7. Average Success Rates for N=100, F=1000 and Q=100, 1000, 10,000

Fig. 8 .

 8 Fig. 8. Comparison of Average Success Rate for N=500, F=100, Q=10,000 Fig. 9 depicts the Average Success Rate versus Popularity of a file.Finally, we also found that the average path length terms of overlay hops) across several simulated 3.5 hops in case of using Gnutella, whereas it is reduced only 1.27 hops away in case of Gnutella with shortcuts. However, socP2P further improves both results the path length to only 1.08 hops.

  of nodes a query is forwarded to Average Success Rate N=100. F=1000 for N=100, Q=1000 and different values of F average success rate with respect to the queries in the network (Q).

Fig. 9 .

 9 Fig. 9. ASR vs. Popularity of Resources for Random wal IV. CONCLUSION AND In this paper a new search algorithm discover required content in P2P network proposed. It exploits the properties of social networks interest links are created between nodes that exhibit similar interest.The approach also exploits from received and routed queries queries can be routed intelligently towards them. simulation results evaluate the performance of the proposed algorithm and demonstrate that socP2P well-known search mechanisms. rate than Random Walk and at the same time message overhead. Finally, socP2P not only the popular content, but also Many possible directions for future on this advance. Creating mutual interest links between nodes that have interest similarity, increasing the search success rates of less-popular files, include overhearing part in forwarded node query and also reducing the search range from two hops to one are some of our possible future endeavors

978-1-4673-2713-8/12/$31.00 ©2012 IEEE

ACKNOWLEDGMENT

This work has been supported by the EU ITEA 10028 "Web-of-Object" (WoO) founded by DGCIS.