Ferréol Soulez
email: ferreol.soulez@obs.univ-lyon1.fr

Loïc Denis

Éric Thiébaut

Corinne Fournier

Charles Goepfert

This paper proposes a micro-particles detection scheme in digital holography. In our inverse problems approach, we estimate the optimal particles set which best models the observed hologram image. Such a method can deal with data which have missing pixels. By considering the camera as a truncated version of a wider sensor, it becomes possible to detect particles even out of the camera field of view. We tested the performances of our algorithm against simulated and experimental data for diluted particle conditions. With real data, our algorithm can detect particles far from the detector edges in a working area as large as 16 times the camera field of view. A study based on simulated data shows that, compared to classical methods, our algorithm greatly improves the precision of the estimated particle positions and radii. This precision does not depend on the particle's size nor location (i.e. inside or ouside the detector field of view).

Inverse problems approach in Particle Digital Holography: out-of-field particle detection made possible

Introduction

Experimental flow has long been studied by imaging particles that are carried by the flow, either naturally present or artificially seeded. The associated techniques are known as Particle Tracking Velocimetry (PTV) when each particle is tracked from one image to another, and Particle Image Velocimetry (PIV) where local flow velocity is computed on a correlation criterion basis, suitable only for high particle concentrations. Velocity and size measurements are still often limited to the plane of the light sheet by the lack of accurate truly 3-D imaging techniques. The 3-D recording capability of holography makes it suitable for 3-D particle imaging and a lot of work has been done in the field of particle holography [START_REF] Kreis | Methods of Digital Holography: A Comparison[END_REF][START_REF] Murata | Potential of digital holography in particle measurement[END_REF][START_REF] Pan | Digital holography of particle fields: reconstruction by use of complex amplitude[END_REF][START_REF] Hinsch | Special issue: Holographic particle image velocimetry[END_REF][START_REF] Sheng | Digital holographic microscope for measuring three-dimensional particle distributions and motions[END_REF][START_REF] Ooms | Optimizing Fourier filtering for digital holographic particle image velocimetry[END_REF]. Digital holography has greatly improved the processing of holograms by suppressing the wet-chemical processing step and opening new frontiers for automated volume analysis. Hologram analysis still remains largely influenced by traditional optical processing and most authors describe analysis of 3-D reconstructed volume obtained by diffraction simulation (Fresnel transform [START_REF] Kreis | Methods of Digital Holography: A Comparison[END_REF], fractional transforms [START_REF] Zhang | Applications of fractional transforms to object reconstruction from in-line holograms[END_REF] or wavelets transform [START_REF] Buraga-Lefebvre | Application of wavelet transform to hologram analysis: three-dimensional location of parti-cles[END_REF][START_REF] Liebling | Fresnelets: New Multiresolution Wavelet Bases for Digital Holography[END_REF]). We refer to those approaches as classical approaches as they are widely used and accepted.

Classical digital holography is performed in two steps: first, a hologram is recorded with a holographic setup directly on a digital camera; then, this hologram is processed by reconstructing a 3-D volume and analyzing the volume. In the case of particle holograms, the most efficient setup is the Gabor setup (Fig. 7) which requires few optical components: a coherent source, either a continuous laser or a pulsed laser synchronized with the camera, and a beam expander. The particles create a set of diffraction patterns named a Gabor hologram. The use of such a setup limits the particle size and the particle concentration to small particles whose projections on the sensor cover less than 1/100 to 1/10 of its surface [START_REF] Royer | Holographic velocimetry of submicron particles[END_REF]. The hologram is then processed in the second step by first reconstructing a volumetric image using diffraction simulation transforms (e.g. volume restitution by Fresnel transform of the zero padded hologram image); and then the reconstructed volume is segmented and the particle 3-D location and size can be deduced from the centroid and mean diameter of each segmented object. Such ap-proaches, however, suffer from various limitations: (i) the depth resolution is about 10 times worse than the transverse resolution; (ii) the field of view is limited and, in practice, must be restricted to the center of the hologram to reduce the border effects; (iii) spurious twin-images of the particles get reconstructed; (iv) multiple focusing [START_REF] Fournier | Digital in-line holography: influence of the reconstruction function on the axial profile of a reconstructed particle image[END_REF] can occur around the actual depth location of each particle.

In a companion article [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF], we recently suggested an inverse problems approach to overcome the issue of accurate location of particles. It however remains that digital holograms are today technologically limited to a few cm 2 , that are several times smaller than the size of holographic plates. Digital holograms therefore suffer from truncation which leads to strong artifacts in the reconstructed volume for particles found close to the borders of the hologram. Out-of-field particles (i.e. particles whose projection is outside the hologram) can not be recovered by classical techniques unless they are a few tenths of a pixels out of the hologram support.

Our preceding method detects particles in three steps: pre-localization using a classical algorithm, location refinement, and particle cleaning. The limitations of classical algorithms can be partially avoided by expressing the particle detection problem as a global optimization problem that is solved locally. However, the first step, based on classical algorithms, restricts the method to particles in the field of the hologram. This field is smaller in the case of digital holography than that of optical holography using a holographic plate and remains a drawback for flow measurements by digital holography. Since our previous algorithm is based on an inverse problems approach, it could be applied to out-of-field particles if only they could be detected in the first step of the method. In the present paper we propose to extend our inverse problems approach to the pre-localization stage of our method so as to achieve detection of out-of-field particles.

The paper is organized as follows. In Sect. 2, we briefly recall the general principle of our method [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF] which detects particles in three steps: pre-localization using a classical algorithm, location refinement, and particle cleaning. In Sect. 3, we also recall the equations of the model of the hologram used in our simulations and in the inverse problems approach . In Sect. 4, we derive an alternative technique to perform the volume exploration. In Sect. 5, we present results obtained with our algorithm on both simulated and experimental data showing that out-of-field particle detection was effectively achieved with almost no loss in precision. The benefits of the new algorithm are summarized in Sect. 6.

General Algorithm Scheme

The particle detection method presented here, takes place in a general particle localization scheme described in the companion article [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF]. The different stages of this algorithm are depicted in Fig. 1. First, a detection algorithm gets an approximate location of a particle, then this location is refined using an optimization procedure, and finally the diffraction pattern of the detected particle is subtracted from the hologram. The hologram is iteratively processed by the three steps and is gradually cleaned from all the previously detected particles. We stop the iterative method and consider that all particles have been detected when the next most likely particle found by the algorithm has aberrant parameters (negligible or negative contribution, too large or too small radius...).

We have shown on real data that the coarse localization step followed by the local optimization step effectively solve the global optimization problem of finding the best size and location of a particle in the (residual) hologram image, even for particles located near the hologram borders. Repeating the detection and localization steps after removal of previously detected particles, allows particles with fainter signatures to be found, in particular particles distant from the camera center.

Since our algorithm is derived following an inverse problems approach, it does not require any explicit direct inversion of the observed hologram: all estimations are made in the data space. As a result, our algorithm has a number of advantages over other existing methods. For instance, our method is insensitive to the twinimage problem inherent in in-line holography. Moreover, our algorithm is even able to deal with truncated field of view and with holograms with saturated or dead pixels. However, the first detection step (based on a classical approach) must be improved in order to detect out of field particles. Fig. 2. Notations used in the hologram model. The parameters p j = {x j , y j , z j , r j } are the position and radius of the j-th particle, (x ′ k , y ′ k ) is the location of the k-th pixel and ρ j,k is the distance between the projection of the j-th particle on the detector and the k-th pixel.

Model of the hologram formation

In this section, we recall the model of the observed holograms derived in our previous paper [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF]. We consider an in-line holography setup (see Fig. 7) where studied particles are illuminated with a collimated laser beam, and the digital camera records both the object wave (scattered light) and the reference wave (laser). We work with small spherical particles and under Fresnel's diffraction approximation [START_REF] Goodman | Introduction to Fourier Optics[END_REF], that is for z 3 j ≫ 4π r 4 j /λ, which is verified for the conditions considered in our simulations and experimental results: laser wavelength λ = 532 nm, particle distances to the detector z j ≈ 250 mm and particle radii r j 55 µm.

In our previous paper [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF], we have shown that, under such conditions and for n particles of radii r j and 3-D positions (x j , y j , z j ), the intensity measured by the detector at position (x, y) is given by:

I(x, y) = I src + I bg -2 I src n j=1 η j Re f j (x, y) + I src n i=1 n j=1 η i f i (x, y) η j f j ⋆ (x, y) (1)
where I src is the image level due to the laser source, I bg is the background level, the real factors η j ∈ [0, 1] account for possible variation of incident energy seen by a particle due to non-uniform laser illumination and f j is the amplitude diffraction term by j-th particle [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF][START_REF] Tayler | Fraunhofer holography applied to particle size analysis: a reassessment[END_REF]:

f j (x, y) = r j 2 i ρ j (x, y) J1 2 π r j ρ j (x, y) λ z j × exp i π ρ 2 j (x, y) λ z j (2)
where i = √ -1, J1 is the first order Bessel function and where ρ j (x, y) = (x -x j ) 2 + (y -y j ) 2 is the distance between the point (x, y) in the detector plane (at z = 0) and the projection (x j , y j ) of the location of the j-th particle on the detector.

Under the considered conditions (low particle density, z j ≈ 250 mm and r j 55 µm), the second order interference terms in Eq. ( 1) can be neglected [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF] and the model simplifies to:

I(x, y) ≃ I 0 - n j=1 α j g j (x -x j , y -y j ) , (3)
with I 0 = I src + I bg , α j = 2 I src η j and:

g j (∆x, ∆y) = r j 2 ρ(∆x, ∆y) J1 2 π r j ρ(∆x, ∆y) λ z j × sin π ρ 2 (∆x, ∆y) λ z j , (4)
where ρ(∆x, ∆y) = ∆x 2 + ∆y 2 is the projected distance.

The model in Eq. ( 3) will be the basis of our particle detection approach in Sec. 4. The resulting hologram expression is a sum of terms depending on the location and size of each diffracting particle. The notations and coordinate system used in our model are summarized in Fig. 2.

Particle detection

Assuming this model of hologram formation, we suggest in this section a new particle detection scheme achieving detection of out of field particles. Determining the optimal set of particles that best models the intensity distribution of the hologram is equivalent to finding the n particle parameters {x j , y j , z j , r j ; j = 1, . . . , n} that minimize the weighted least-squares penalty:

P n = x ′ k y ′ k w(x ′ k , y ′ k ) [m n (x ′ k , y ′ k ) -d(x ′ k , y ′ k )] 2 (5)
where

(x ′ k , y ′ k ) are the coordinates of the k-th pixel, m n is the model, d(x ′ k , y ′ k ) is the k-th pixel value in the observed hologram and w(x ′ k , y ′ k ) is its statistical weight: w(x ′ k , y ′ k ) = 1 Var(d(x ′ k , y ′ k )) if k-th pixel is measured, 0 otherwise. ( 6
) where Var(d(x ′ k , y ′ k ))
is the noise variance at pixel (x k , y k ). The noise considered here is meant to account for the detector noise and for the photon noise but Var(d(x ′ k , y ′ k )) can also be set to account for artifacts (such as fringes due to the detector glass window) and for approximation errors (such as digitization and such as speckle noise due to undetected objects at this stage or due to the neglected interferences between particles).

In the case where the noise is assumed constant all over the known pixels,

1 Var(d(x ′ k ,y ′ k )
) can be factorized out of Eq. ( 5). In this case w(x k , y k ) can be viewed as a binary mask on the data:

w(x ′ k , y ′ k ) = 1 if k-th pixel is measured, 0 otherwise. (7)
The model m n (x ′ k , y ′ k ) for n particles is given by Eq. ( 3):

m n (x ′ k , y ′ k ) = I 0 - n j=1 α j g j (x j -x ′ k , x j -y ′ k ) . (8)
As our algorithm proceeds iteratively for each particle (as described in Fig. 1), we consider the previously detected particles parameters as fixed. In our iterative approach, locating the n-th particle consists of determining the parameters of the n-th particle and the constant level that minimizes P n . This is achieved by fitting the parameters {α n , I 0 , x n , y n , z n , r n } on the centered residual image r n-1 :

r n-1 (x ′ k , y ′ k ) = d(x ′ k , y ′ k )+ n-1 j=1 α j g j (x j -x ′ k , y j -y ′ k )-c n-1 , (9)
which is the hologram image after having subtracted the contribution of the previous n-1 particles and a constant c defined as the mean gray level of the (n -1)-th residual image:

c n-1 = x ′ k y ′ k w(x ′ k , y ′ k )   d(x ′ k , y ′ k ) + n-1 j=1 α j g j (x j -x ′ k , y j -y ′ k )   .
(10) In other words, this constant is such that:

σ w r def = x ′ k y ′ k w(x ′ k , y ′ k ) r n-1 (x ′ k , y ′ k ) = 0 , (11)
which is used to simplify the equations in what follows.

For the detection/location of the n-th particle, our misfit criterion writes:

P n = x ′ k y ′ k w(x ′ k , y ′ k ) [I n -α n g n (x n -x ′ k , x n -y ′ k ) -r n-1 (x ′ k , y ′ k )] 2 (12)
with

I n = I 0 -c n-1 .
For the first particle, we take r 0

(x ′ k , y ′ k ) = d(x ′ k , y ′ k ), and c 0 = x ′ k y ′ k w(x ′ k , y ′ k ) d(x ′ k , y ′ k ).
The expression in Eq. ( 12) can be expanded as:

P n = σ w r 2 + 2 α n σ w r g -2 α n I n σ w g + I 2 n σ w + α 2 n σ w g 2 (13
) where it is accounted for the fact that σ w r = 0 (see [START_REF] Fournier | Digital in-line holography: influence of the reconstruction function on the axial profile of a reconstructed particle image[END_REF] and where:

σ w = x ′ k y ′ k w(x ′ k , y ′ k ) ( 14
)
σ w r 2 = x ′ k y ′ k w(x ′ k , y ′ k ) r 2 n-1 (x ′ k , y ′ k ) (15)
σ w r g = x ′ k y ′ k w(x ′ k , y ′ k ) r n-1 (x ′ k , y ′ k ) g n (x n -x ′ k , y n -y ′ k ) (16) σ w g = x ′ k y ′ k w(x ′ k , y ′ k ) g n (x n -x ′ k , y n -y ′ k ) ( 17
)
σ w g 2 = x ′ k y ′ k w(x ′ k , y ′ k ) g 2 n (x n -x ′ k , y n -y ′ k ) . (18)

A. Optimal Parameters Determination

Given x n , y n , z n and r n , the optimal offset I + n and factor α + n for such a model are obtained by solving the linear system:

       ∂P n ∂I n In=I + n ,αn=α + n = 0 ∂P n ∂α n In=I + n ,αn=α + n = 0
which yields:

I + n = -σ w g σ w r g σ w σ w g 2 -σ 2 w g = α + n σ w g σ w , (19)
α + n = -σ w σ w r g σ w σ w g 2 -σ 2 w g , (20)
Replacing I n and α n by their optimal I + n and α + n values in Eq. ( 13) yields the partially optimized criterion:

P + n def = P n | In=I + n ,αn=α + n = σ w r 2 - σ w (σ w r g ) 2 σ w σ w g 2 -σ 2 w g . (21)
It remains to optimize P + n with respect to the other parmaters, namely x n , y n , z n and r n . This is equivalent to maximizing:

Q n def = σ w (σ w r g ) 2 σ w σ w g 2 -σ 2 w g (22)
under the condition (by definition α j = 2 γ A 0 2 η j > 0):

α + n = -σ w σ w r g σ w σ w g 2 -σ 2 w g > 0 .
(23)

B. Fast Computation

It is possible to achieve fast determination at pixel resolution of the projected coordinates (x ′ n , y ′ n ) of the n-th particle given its other parameters z n and r n . The terms σ w r g , σ w g 2 and σ w g (given respectively by equations Eq. ( 16), Eq. ( 17) and Eq. ( 18)) are involved in the expression of the criterion Q n to maximize. When z n and r n are chosen, these terms are simply 2-D discrete convolutions (the term σ w is a constant which can be computed once for all). Under circulant approximation (as CCD is a rectangular grid), it is possible to compute these terms by a few fast Fourier transforms (FFT's) for any projected coordinates (x n , y n ) matching the locations of the CCD pixels. In other words, given the depth z n and radius r n of the n-th particle, it is possible to quickly compute a 2-D map of the criterion Q n with respect to the projected coordinates and with a resolution equal to the size of a pixel. The maximum of this map indicates the most likely projected positions of the n-th particle if it were located at distance z n with a radius r n . In the following, we detail the computation of the 2-D map of Q n .

We denote g the vector of values of g n for given depth z n and radius r n and sampled at projected coordinates equal to the CCD coordinates:

gk def = g n (x ′ k , y ′ k ) . ( 24
)
Similarly we denote by r the 2-D map of residuals:

rk def = r n-1 (x ′ k , y ′ k ) d k is measured, 0 otherwise, ( 25
)
and by ẘ the 2-D map of non-stationnary weights:

ẘk def = Var(d k ) -1 d k is measured, 0
otherwise (mainly outside the field of view).

(26) We purposely defined r and ẘ for a field of view larger than the CCD support. The next section explains how support and aliasing issues are solved by using such an extended field of view. Note that the value of r outside the field of view is irrelevant that's why the weights ẘ must be set to zero outside the field of view (and for any bad pixels).

Under the circulant approximation, the 2-D map of σ w g can be approximated by:

å def = F -1 • diag(F • ẘ) • F • g ( 27
)
where F is the 2-D discrete Fourier transform (DFT) matrix operator and where diag(u) is the diagonal matrix with diagonal elements equal to the elements of vector u; hence diag(u) • v is simply the element-wise product of vectors u and v:

[diag(u) • v] k = u k v k . ( 28
)
Under circulant approximation, å represents the correlation between ẘ and g. Likewise, the 2-D maps of σ w r g and σ w g 2 can be respectively approximated by:

b def = F -1 • diag(F • diag(ẘ) • r) • F • g , (29)
c def = F -1 • diag(F • ẘ) • F • diag(g) • g . (30)
Finally, the 2-D map of the criterion Q n is approximated by:

Q k ≈ Qk = σ w bk σ w ck -å2 k (31)
which is the value of Q n for a projected position (x n , y n ) = (x ′ k , y ′ k ). For every particle and for every sampled depth and particle radius, a circulant approximation of the 2-D map of Q n must be computed. This requires 7 FFT's (4 direct FFT's plus 3 inverse FFT's). In fact, the DFT of ẘ can be computed once for all, whereas the DFT of diag(ẘ) •r must be computed for every new particle, hence only 5 more FFT's are required for every particle and for every sampled depth and particle radius. With such a formulation, it is possible to seek for particles outside of the field of view in a working area T which is larger than the CCD. However, for fast computation of P n , we use fast Fourier transforms which yield spatial aliasing because of the circulant approximation. To limit the aliasing effects, we perform all computations on a region R which is larger than the working area T where particles are to be detected. The width and height of the region R are chosen so that W R ≥ W T + W and H R ≥ H T + H where H T × W T and W × H are the dimensions of the working area T and of the detector (cf. Fig. 3). The residuals and statistical weights over the region R are set according to Eq. (25) and Eq. (26).

C. Reconstruction Outside of the Field of View

D. Segmentation

Laterally, P n is sampled with the same size as pixel width inside and outside of the field of view. P n is computed for several slices (z n , r n ) with a coarse sampling (typically every mm in depth and every 5µm for radius). This four dimensional space of solutions is exhaustively explored in order to extract parameters {x ′+ n , y ′+ n , z + n , r + n } that minimize P n . In practice, to reduce computation, only a three dimensional space of solutions (x ′ , y ′ , z) is evaluated using an average radius. For improved efficiency, the mean radius can be estimated using the direct hologram analysis described in reference [START_REF] Denis | Direct extraction of mean particle size from a digital hologram[END_REF]. Eventually, at the optimal depth z + n , P n can be computed for a few r n values depending on particle radius dispersion (see Sec. 5).

Results

A. Simulated Data

To assess the performance of our algorithm, we first processed simulated data. We compare the results obtained with the new algorithm and the results obtained with a classical approach. A study of the accuracy versus transversal location of the particle was conducted. This data set is composed of 32 holograms simulated with the same experimental conditions. Each hologram records the diffraction pattern of 20 disks with random radii between 40 and 50 µm and randomly distributed in a box of 13.72×13.72×50.00 mm 3 . The wavelength of the laser is 0.532 µm. The camera has 1024 × 1024 pixels spaced by 6.7 µm and is placed at about 250 mm from the studied volume in order to stay in the approximation range defined in section 3. The field of view of this camera (6.86 × 6.86 mm 2 ) represents a quarter of the section of the studied volume and is centered on it.

Our inverse problems approach is based on a simplified model of the hologram image formation. To account for the effects of such simplifications on the detection and localization of the particles, we use a more accurate model to simulate the hologram images used in this test. More precisely, we use the image model given by Eq. ( 1) to simulate the holograms whereas the simplified model in Eq. ( 3) is assumed during the detection and localization steps of the algorithm. Hence the interference pattern between waves diffracted by different particles is accounted for in the simulated images whereas it is neglected by the algorithm. In order to properly simulate quantization by the detector, the dynamics of the simulated hologram images is converted into 8-bit integer values prior to their processing.

In the particle detection step, a reconstruction every millimeter in z was sufficient for a successful coarse detection. These 50 detections were performed for an average radius of 45 µm. At the depth z + that minimizes P n , 5 detections were performed for equally sampled radii between 40 and 50 µm. This leads to (55 × 5 + 2) FFT computations per particle.

The results of this inverse problems approach were compared with a classical particle detection algorithm (based on hologram reconstruction followed by 3-D segmentation). As such conventional methods can only deal with particles inside (or very close to) the camera field of view, we conducted the tests on simulated holograms under exactly the same conditions to have approximately the same number of detected particles. To achieve meaningful statistical results, we performed all the tests on 100 simulated holograms.

Our method detected almost every particle: on average, only 20 particles out of 640 were not detected (detection ratio of 97 %) and 2 were falsely detected (false positive ratio of 0.3 %). The conventional method detects 631 particles out of 2000 (detection ratio of 32 %). As expected, we noticed that such a method cannot properly deal with out of field particules (detection ratio of 100 % in the field of view and 6 % outside). Comparison of lateral error (on x and y) for classical and "inverse problems" approach (pixel size: 6.7 µm).

Errors on particle parameters (x, y, z, r) are plotted as a function of the distance of the particle to the closest edge of the CCD. This representation seems to underline well algorithm artifacts. Results are shown in Fig. 4, Fig. 5 and Fig. 6.

• Transversal errors of both methods are plotted in Fig. 4. Sub-pixel precision is reached for both methods but our method is unbiased and more precise. It presents an equal standard deviation on both x and y axis (∆x ≃ ∆y ≃ 0.3 µm or 1/20 pixel). This standard deviation does not significantly depend on the size nor the location of the particles.

• Fig. 5 shows the depth error for both methods. Our method improves greatly the classical one: the bias is negligible and the sub-pixel standard deviation (∆z ≃ 0.9 µm or 1/7 pixel) is almost constant whatever the particle parameter.

• Only errors on radius for our approach are plotted in Fig. 6 as radius estimation by the classical method presents errors too large to be plotted on the same figure (bias is -8.8µm and ∆r ≃ 1.8 µm for classical method). Radius estimation by our method is very precise with a standard deviation of ∆r ≃ 0.04 µm for an average radius of 45µm. Although radius estimation is unbiased for particles outside of the field of view, a small negative bias can be noticed for particles inside. This bias is due to the simplifications made in our model: the neglected interference term (second term in Eq. ( 1)) is larger when particles are inside of the field of view.

B. Experimental Data

We carried out an experimental test of our algorithm using real data from an in-line holography setup. The experimental test data are the same as in the companion article [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF]. The experimental layout is shown by Fig. 7.

The camera is a 12-bit CCD with 1280 × 1024 pixels of size 6.7 µm × 6.7 µm. It is located at about 25 cm from the injector in order not to disturb the flow experiment. This leads to a small but realistic numerical aperture of Ω = 0.014. The laser wavelength is 0.532 µm. Studied particles are water droplets injected by a piezoelectric injector (conceived by Bremen University) which generates monodisperse droplets. Their sizes are adjustable between 50 to 100 µm.

The experimental dataset consists of a set of 100 pairs of holograms with four or five droplets on each image. An example of such holograms can be seen in the central box of Fig. 9. The 3-D locations and radii of droplets were extracted from this dataset by our algorithm. Each hologram was processed as described Sec. 2. At each iteration, P is computed to detect the next particle. Paste-up Fig. 8 shows an example of a slice (x, y) with calculated P at each iteration. Fig. 9 illustrates the result of the processing of one hologram of the dataset. Let us point out that although four particles were in the camera field of view (in the box on the figure), up to 18 were detected.

The measured particles locations are shown in Fig. 10. The average trajectory clearly is close to the straight line expected from the experimental conditions. Detected out-of-field particle locations and their spacings describes well their ballistic trajectory. At each extremity of the reconstructed jet, the low signal to noise ratio increases the observed deviations from the ideal straight trajectory. For particles far from sensor edges, assumptions made in Sec. 3 become less relevant and model noise is added to measurement noise. The jet divergence which can be seen in Fig. 10 is mostly due to vibrations of the injector [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF]. Hence, the effective precision of the measured locations is smaller than the variations due to these physical effects.

The droplet sizes estimated by our algorithm have a bell-shaped distribution (see Fig. 11) with a mean diameter of 94.1 µm and a standard deviation of 0.3 µm in agreement with the settings of the droplet injector.

Conclusion

In a companion article [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF], we have introduced a new algorithm that accurately estimates the parameters of particles in digital holography. In this paper, we complete this scheme by describing a new algorithm for the detection of particles outside of the field of view of the sensor. The most important difference from other existing techniques is that our processing is based on an inverse problems approach and thus does not require any direct inversion. In this framework, we introduce a simplified model of the hologram images which directly incorporates the sizes and locations of the diffracting particles. We solve the problem by determining the set of particle parameters for which the difference between the model and the data is minimized. The objective function turns out to have multiple local minima and thus a global optimization procedure is required to properly solve the problem. Our algorithm effectively achieves global minimization by performing an approximate detection of the particles in the whole parameter space followed by a local refinement. By repeating these steps on the residual images, our algorithm is able to detect particles even if their signature is faint compared to the speckle noise due to the other particles in the original hologram.

We have tested our algorithm on both simulated and real data. Our results show that the precision along the depth direction is greatly improved and is much better than the optical resolution in such conditions (δz ≥ λ/Ω 2 = 2.6 mm). In simulations our algorithm achieves unbiased localization with constant accuracy even far from the sensor edges. Comparison with classical methods shows very strong improvements (laterally ∆x = ∆y = 0.3µm, ∆z < 1µm in depth and ∆r < 0.05µm for radius). When dealing with experimental data, it appears that our algorithm is also robust with respect to non-homogeneous illumination and to spurious patterns due to experimental noise, as seen in Fig. 9. As the estimation is done in the data space, our algorithm can cope with truncated holograms or bad pixels: it is sufficient to set to zero the weight of unmeasured or bad pixels. Such a formulation allows us to detect particles outside the classical field of view, effectively extending the area of measurement. For instance, Fig. 9 shows that, even though only six particles signatures could be seen on the CCD, our algorithm was able to accurately localize 12 additional particles, using non-local information in the detected image.

In its current implementation, the computation burden of our algorithm is quite heavy: for our simulation, the full processing (detection, refinement and cleaning) takes about 7 minutes per particle on a Pentium IV CPU at 3.60 GHz. For routinely use, we plan to accelerate the method. This could be done if we manage to perform multiple-particle detection per pass. This improvement but also a more extensive analysis of perfomances and possible limitations of our algorithm are currently under study and will be part of forthcoming paper. For instance, as the density of particles augments, one expects that interferences between the wave diffracted by close particles to become more important which could break the approximation made in the detection step. We have checked [START_REF] Soulez | Inverse problem approach for particle digital holography: accurate location based on local optimisation[END_REF] that our method works for at least 100 particles over 1000 × 1000 pixels but we did not tried to push this limit. Nevertheless, in-line digital holography methods are clearly devoted to diluted particle conditions (because they involve the recording of 3-D information on a 2-D detector with no phase information), the occurrence of pairs of very close particles is very unlikely in these conditions.

Fig. 1 .

 1 Fig. 1. Synopsis of the method.

Fig. 3 .

 3 Fig. 3. Computing region R and working area T compared to the CCD size (W×H)

  Fig.[START_REF] Hinsch | Special issue: Holographic particle image velocimetry[END_REF]. Comparison of lateral error (on x and y) for classical and "inverse problems" approach (pixel size: 6.7 µm).

Fig. 5 .

 5 Fig.5. Comparison of depth error for classical and "inverse problems" approach (pixel size: 6.7 µm).

Fig. 6 .

 6 Fig.6. Error on radius for the "inverse problems" approach (pixel size: 6.7 µm).

Fig. 7 .

 7 Fig. 7. In-line holography setup

Fig. 8 .

 8 Fig. 8. Paste up representing a slice (x, y) of P computed on an experimental hologram. It is composed of 13 images (numbered from 1 to 13). Each image is computed on the residuals after the detection of the particles (pointed by arrows) detected in the previous images.

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Superimposition of an experimental hologram (in the box) on an synthesized hologram of 18 reconstructed particles.

Acknowledgments

The authors are very grateful to the anonymous reviewers whose fruitful remarks have helped to improve the paper.

The algorithms and the simulations presented in this article have been implemented in Yorick (http://yorick.sourceforge.net), a free data processing language written by D. Munro.


Antoine Clement
email: aclement@mines-albi.fr

Daouda Kamissoko
email: dkamissoko@mines-albi.fr

Franc ¸ois Marmier

Didier Gourc
email: dgourc@mines-albi.fr

The robustness of the tour planning face to risks of road freight transport

Keywords: risk, assessment, decision, tour, robustness

Currently, over 90% of goods are transported by road terrestrial mode. To progress, companies of Road Freight Transport sector (RFT) are subject to a number of issues, whether economic, regulatory, environmental, societal, occupational health and safety ... The planner has a crucial role in achieving these challenges, it is he who undertakes the construction method of the tours, for the future of the company. Thus he bears a heavy responsibility on his shoulders. Its ability to integrate the company's constraints and to master the balance between these issues in the construction of tour is vital. However a list of risks is potentialy present and can affecte his planning. This lack of vision about risks impact prevents him to achieve his mission. Thus we propose in this article a mean to qualify the tours generated in relation to the issues mentioned but also in terms of potential hazards that could come impair the quality of tours. In this article, we propose an evaluation process of road transport of goods tours and an indicator to measure the robustness of touring.

I. INTRODUCTION

Nowadays, more than 90% of goods are conveyed by road. Road Freight Transport (RFT) companies therefore embody a sector of capital importance to the French economy; they represent an annual turnover of 44 billion euros, which is generated by more than 37,200 companies (more than 90% of companies in this sector are SMEs/VSEs with less than 50 employees) and roughly 420,000 employees. This sector is characterized by strong competition and growing pressure from customers and suppliers, tight delivery times, exacerbated flexibility, etc. (economic stakes). RFT companies are also subjected to specific national and community standards and regulations (regulatory and environmental stakes). In parallel, they are required to fulfil performance duties in terms of preventing risks of occupational accidents and diseases (article L4121-1 of the French labor code) (occupational health-safety issues). Simultaneous consideration of these three obligations has led to adoption of the term "sustainable transport systems".

There are Transport Management System (TMS) aids that ease the pressure on the planner, but these only partially integrate the area of sustainable transport operations. TMSs are generic solutions mainly aimed at transport companies that deliver to regular customers. These solutions are unsuitable for integrating specific characteristics and for the diversity of customer demands, to which the transport company must respond, if it is not to lose the contract. Such solutions do not allow fleets of mixed vehicles to be managed (optimization of vehicle selection based on its energy autonomy and/or city/town center, suburban or regional customer locations). Driver allocation to delivery rounds is not systematically integrated into TMSs, although it is crucial to limit driver overtime and to equitably distribute workload among drivers. Ultimately, TMSs mostly ignore both health and safety aspects and consideration of different combined energies. Moreover, these aids are relatively expensive because they are usually divided into modules or systems; a purchase or subscription is required for each module, function or hardware item (the average cost of a rented module being 30 euros/month/truck and a module does not operate separately). Yet, cost is a major constraint, especially for VSEs/SMEs. All these curbs (unsuitability to customer diversity, to mixed fleets, no overtime management, cost of systems) oblige many VSEs/SMEs to perform their planning using a spreadsheet, a word processing program or even pencil and paper.

New proposals are required for planning sustainable transport services and these must take into account the financial, environmental and social issues.

Besides taking criteria into account for the quality evaluation of the generated tours, we would also propose solutions to integrate the measurement of the robustness of the solution to potential risks that could alter the planning done. In front of the complexity of its task, we should be able to suggest the planner some indicators to facilitating the selection and the choice of the best solution. Among these indicators, we propose to build a representative indicator of the confidence that can be attributed to a solution. By confidence we mean the ability of this solution to keep the commitments made, for example with the customer or compliance with the maximum working time (time amplitude of the tour).

This article is divided into 5 parts. In section 2, we present a state of art for the selected issue. In section 3, the part of the methode to describe the work developed to address the problem. Section 4 contains an example to illustrate the proposals made. To finish, the conclusion takes stock of the situation and provides some ideas for complete this work.

II. STATE OF THE ART

A. Related work about risk analysis

For several years, we have led, in our laboratory, research dealing with the theme of analysis and control of risks in the project and process management. These works led to the development of innovative proposals based on the R-IO model [START_REF] Nguyen | A decisionmaking tool to maximize chances of meeting project commitments[END_REF][START_REF] Marmier | A risk oriented model to assess strategic decisions in new product development projects[END_REF][START_REF] Marmier | Strategic decision-making in NPD projects according to risk: application to satellites design projects[END_REF]. These works aim to take into account all combinations of occurrence risks and their influence on schedules. Scenarios of schedules thus generated are evaluated according several criteria: probability of occurrence, duration, cost, and any other specific criteria for a type of evaluated schedule.

In the literature, risk management methods refer to a standard process with well-known steps: identification, assessment, quantification, treatment and risk monitoring [START_REF] Pmi | PMBOK Guide[END_REF]. [START_REF] Tixier | Review of 62 risk analysis methodologies of industrial plants[END_REF] propose a classification with 62 existing approaches. Methods are sorted depending on whether they are deterministic and/or probabilistic, but also qualitative or quantitative. Gourc et al. offer a grid of different reading, classifying approaches of risk management whether they are symptomatic or analytical [START_REF] Gourc | Vers un modèle général du risque pour le pilotage et la conduite des activités de biens et de services : Propositions pour une conduite des projets et une gestion des risques intégrées[END_REF]. The first group of approaches, also called risk-uncertainty is associated with approaches where the project risk management is converted into managing uncertainty [START_REF] Ward | Transforming project risk management into project uncertainty management[END_REF]. The second group of approaches views risk as an event that may affect the achievement of project objectives [START_REF] Carter | Introducing RISKMAN: the European project risk management methodology[END_REF].

Two themes are well known for their reference to innovation and by the omnipresence of risk: (1) the project management of new product development (NPD) at an operational level. It consists of choosing an exclusive orientation as a development strategy. It considers a level of global tolerance of risk. (2) NPD portfolio management at a tactical level. A first definition of portfolio management is given by [START_REF] Cooper | New Product Portfolio Management: Practices and Performance[END_REF]: a dynamic of decision making that enables the list of projects to be always updated. In this process, new projects are evaluated, selected and sorted. Past projects can be accelerated, stopped, or paused and resource assignments can be modified in order to build a balanced portfolio of innovative projects, on a balance risk / profitability.

Representations of project plan and tour are wery simimar, all of them can be modelised with activity networks. So, we hypothesize that concepts and software tools developped on project risk management by our research team since a few years are totally relevant for assessment of the robustness of tour in transport sector.

B. Metric to assess the planning robustness

In the literature, the robustness notion is sometime used to talk about algorithm. But this review is based on the robustness of solutions. [START_REF] Ferrario | Evaluation of the robustness of critical infrastructures by Hierarchical Graph representation, clustering and Monte Carlo simulation[END_REF] point out that there is no unique definition of the term robustness. To qualify solutions, several authors propose definitions of robustness ([11, 12, 13]) define the robustness as the aptitude to which a system can function correctly in the presence of inputs different from those planned. [START_REF] Jen | Stable or robust? What's the difference?[END_REF][START_REF] Ali | Definition of a robustness metric for resource allocation[END_REF] describe that the concept of robustness should be defined for a given set of features and under a given set of possible disturbances. There is important to note that the Society for Risk Analysis recently proposed a definition of robustness [START_REF] Sra | SRA develops glossary of risk-related terms[END_REF]. A system is robust to perturbations when goals can be achieved despite disturbances occur.

In the context of transport plannification, robust planning is defined as a planning that is insensitive to disturbances that might occur during the implementation of planned activities. The definition of a robustness metric aimed to measure the sensitivity/insensitivity of the planning to potential disturbances. The robustness can be assess using metrics. Those metrics must enable to reflect a-priori the planning's capacity of the tour to respect objectives or requirements that we define as the contractual space. Thus, in our context, the contractual space will be defined, for instance, by the tour's maximum working time, the maximum time for each driving period, etc.

The directed acyclic graphs can be used to model schedules and, by extension, tour transport's schedules. In the context of planning for directed graphs acyclic, [START_REF] Shi | Robust task scheduling in non-deterministic heterogeneous computing systems[END_REF] offer two definitions of the robustness and propose two expression of the robustness. The first is expressed by the multiplicative inverse of the expected value of relative schedule tardiness. The second is calculated by the multiplicative inverse of schedule miss rate.

This work addresses primarily deterministic graphs. For our problem, we are led to consider the probability of occurrence of hazards during the tour. These hazards can modify some characteristics of the graphs (duration of activities, modification of some links).

[18] and [START_REF] Canon | Evaluation and Optimization of the Robustness of DAG Schedules in Heterogeneous Environments[END_REF] propose a comparative analysis of several metric robustness. Metrics are:

• "makespan standard deviation" • "makespan differential entropy"

• "slack mean"

• "probabilistic metric"

• "lateness likelihood"

• "makespan 0,99-quantile" [START_REF] Smith | Measuring the Robustness of Resource Allocations in a Stochastic Dynamic Environment[END_REF] and [START_REF] Shestak | A Stochastic Approach to Measuring the Robustness of Resource Allocations in Distributed Systems[END_REF] also provide metrics taking the dynamic nature of environment into account for which schedules are generated. The robustness is expressed as the probability that the task performed, i.e. the tour, is finished in accordance with constraints, for instance the total duration or the delivery time.

III. PROPOSALS

Based on the results we obtained earlier, particularly in the ANR project PLAS'O'SOINS [START_REF] Bastide | Plas'O'Soins: A software platform for modeling, planning and monitoring homecare activities[END_REF] related to the tour's planning of health workforce for the maintenance of patient Figure 1. Proposed approach at home [START_REF] Gourc | An approach based on Tabu search technique for solving a multi time window home healthcare scheduling problem[END_REF] and our works on R-IO. The aim of the project PLAS'O'SOINS was to improve the actor's planification in health field with take in consideration risks . We propose to evaluate the quality of schedules generated by relying on the measure of possible scenarios induced by the possible occurrence of the identified hazards. For this purpose, we implement the ProRisk model on the schedule of activities concerning a tour.

A. Suggested approach

The proposed approach (Figure 1) consists of 5 steps. The first step consists in building schedules tours to serve the identified orders. Thus, on the basis of received orders, this step enables to produce a tour schedule set represented by process P described by their tasks T i , (i = 1...t). Based on the identified generic risks and planned activities in these processes, step 2 produces a list of specific risks for each tour and generates all the possible risk scenarios. In step 3 the planner choose plannings thus he wants to evaluate accord criteria selected. In step 4, the robustness are calculated to determine a other indicator, probability of commitments' respect, for each tour initially generated by step 1. Finally, step 5, conducted by the planner, aims to select the tour that best meets these requirements

B. Model proposal

A tour is presented as a process P composed by its tasks T i , (i = 1...t), where t is the number of tasks (Figure 2). A tour, mobilizing a vehicle and a driver, is a process involving a succession of activities alternating with:

• starting activity, corresponding at the begin of the tour with loading of the vehicle and preparing; • driving activities, corresponding at all mouvment to connect two points, warehouses or on a delivery points; • waiting activities, corresponding at the time to wait before delivery a site; • activities on site (AoS), corresponding at loading and unloading goods on site and administratif part; • closing activity, corresponding at the end of the tour with unloading goods and cleanning the vehicle. A tour is considered feasible or acceptable if it meets a set of constraints C. The constraints can include for example the respect of a maximum duration of the tour, the respect of a time slot to a client, etc.

A tour is represented as a process P described by its tasks T t , (t = 1...T ), where T is the number of tasks. The planning process provides a baseline P i which does not include the risks. A set of critical resources (ERC) is required to process certain tasks. A process is also described by its set ER of identified risks R i , (i = 0...n), where n is the number of risks. Each R i is characterized through the risk management process. A risk R i is also characterized by its occurrence period, namely the tasks where the risk can occur. Its probability proba(R i ) is the probability that the event related to Ri happens. Its impacts on costs are rated CI(R i ) and the project duration DI(R i ). The impacts may affect, by the domino effect, a different task from the one characterizing the period of occurrence. These probabilities and impacts are also named initial probability and initial impact. A risk scenario ScR corresponds to the combination of S risks occurring during a process P. A process with n risks leads to 2 n risk scenarios. ScR s , (s = 1, ..., 2 n ) is one possible implementation with k risks (0 ≤ k ≤ n). The total number of risk scenarios, presenting k risks is equal to n!/k(n -k)!. Its probability is noted proba(ScR s ) (the probability that the events related to this risk scenario occur and other risks do not occur). Depending on type of tasks composing the process, certain generic risks can be identified, see Table I for an example.

A risk can have an impact on the strong constraints and thus make a scenario become not feasible but it also simply affect the indicators without modifying its feasibility.

Based on these hypothesis, we are implementing a ProRisk analysis process allowing to generate all possible scenarios and calculate the indicators characterizing each scenario. The aim is, first, to determine the probability of occurrence of risk scenario proba(proba(ScR s ) . This probability of the scenario depends on the combination of risks in this scenario (Equation 1).

proba(ScR

s ) = n i=1 proba(R i ) if R i ∈ ScR s 1 -proba(R i ) if R i / ∈ ScR s (1)
Beside the occurrence probability of a scenario, we characterize each scenario by other indicators. Very typically, the time and cost indicators are mobilized to effect this characterization. Thus, for each identified scenario, its duration and its cost are calculated. The expression of the cost calculation is shown in Equation 2. The project scenario cost is named C(ScRs). It includes the cost of tasks C(T t) which constitute risks present in ScR s , on the process tasks and the global cost GC initial (R i ) of occurring risks. GC initial (R i ) is the aggregation of risk impacts on the cost criterium. This includes a fixed portion of the total cost (materials, tools, parts, etc.) and indirect costs depending on the duration of action of each possible impact of risk

C(ScR s ) = T t=1 C(T t ) + Ri∈ScRs t GC initial (R i ) (2)

C. Expression of the level of robustness

As a first approach, we propose to define a robustness indicator of each tour planning as the measure of the ability to respect contractual space expressed by the defined constraints. Thus, robustness Ro(i) of the tour i, will represent the probability to respect this set of constraints Co. It is obtained from the occurrence probabilities of risk scenarios ScR s , which respect constraints Co. This set of scenarios will be named Sc in (Figure 3).

Ro(i) = ScRs∈Scin proba(ScR s ) (3)
At this stage, this expression represents a first formulation of a robustness indicator. We expect to offer a finer indicator to weight the importance of the scenarios that don't respect the constraints and maybe take in consideration indicators limits.

IV. NUMERICAL APPLICATION

A. Description of the example

For half a day planning, our partner may have to deliver three customers with one vehicle up to 200 customers with 20 vehicles. To illustrate our proposals, we consider an example situatein a low range of a real planning with 4 delivery points (Sit j ) and two couple vehicle, driver. We consider that vehicle are different the first one is an electrical vehicul and the second one is a diesel vehicle. These four delivery points are spread over a geographical area compatible with the vehicle's energy. The vehicle is stored in a warehouse (W 1 ), which also stores a portion of workpieces. Each point, warehouse and delivery point, can potentially be a place where loading and unloading's actions of the vehicle will be to achieve.

The list of orders to be processed is shown in the table (Table III). There are three orders starting at the warehouse and one where the picking point is the Sit 3 .

For each order are described:

• The place of loading of objects to transport (Origin);

• The place of unloading of objects to transport (Dest.);

• The loading time and unloading time calculated based on quantity; • Time window where the driver can delivery goods. Additional information can be mobilized, such as:

• Characteristics of the objects to be transported: volume in liter (Vol.), mass in kg, number of objects to be transported (Qty) per unit carton/ pallet. • Presence of dock. Two generic risks are considered, as described in Table II. The risk R 1 expresses the effects of the traffic jam when accessing to Sit 3 . The risk R 2 reported at a possible absence of the handler from his workplace if the delivery takes place between 10h00 am and 10h30 am. It is frequent, but not systematic, the handler in charge of loading / unloading the truck, takes his break at that time.The risk R 3 expresses the possibility that a driver does not know the localisation where goods need to be delivery just for the Sit 1 .

In a complementary manner, the constraints of the transport company, inherited from legal and regulatory constraints require that:

• Time allowed for a tour does not exceed 5h; • Driving periods shall not exceed not 1h;

• The autonomy of electical vehicle is 4h.

B. Results and discussion

To comply with the orders to be processed and any eventual constraints, for example it's an obligation to go at the delivery point Sit 3 before the delivery point Pli4, there are 40 solutions for the tour. A part of this 40 solutions are shown in Table IV. This table shows the succession of points to visit, starting from the warehouse W 1 and go back to the warehouse W 1 .

The initial schedule proposed for the solution S 1 is shown in Figure 4, it shows the details of the solution S 1 , indicating for each activity required:

• Its type: start, driving loading, unloading, ending of the tour when come back at the warehouse; • Duration of each activity ;

• The activity's start time. This solution can be characterized by many criteria:

• Time allowed;

O 1 W 1 Sit 1 - 00h10 [8h-13h] O 2 W 1 Sit 2 - 00h15 [8h-11h] O 3 W 1 Sit 3 - 00h10 [8h-13h] O 4 Sit 3 Sit 4 00h05 00h05 [10h-13h] Table IV SOLUTIONS TO REALIZE DELIVERY ORDERS Solution S 1 V 1 W 1 Sit 1 Sit 2 Sit 3 Sit 4 W 1 V 2 - - - - - - ... ... ... ... ... ... ... ... S 20 V 1 W 1 Sit 1 Sit 2 W 1 - - V 2 W 1 Sit 3 Sit 4 W 1 - - S 21 V 1 - - - - - - V 2 W 1 Sit 1 Sit 2 Sit 3 Sit 4 W 1 ... ... ... ... ... ... ... ... S 40 V 1 W 1 Sit 3 Sit 4 W 1 - - V 2 W 1 Sit 1 Sit 2 W 1 - -
• The maximum driving time;

• Total driving time;

• Autonomy of the vehicle V 1 ;

• Total time for loading and unloading;

• Environmental impact of the solution.

In the example discussed in this article, only some of these criteria are implemented and associated values are take in consideration.

From these results, the solutions do not respect the constraints set forth above may be dismissed: respect of the maximum amplitude of 4 hours, etc... For the second step of the proposal three risks are created. Eight risk scenarios may occur:

• None risk: ScR • Only the risk R 1 occurs: ScR R1 ;

• Only the risk R 2 occurs: ScR R2 ;

• Only the risk R 3 occurs: ScR R2 ;

• Risks R 1 and R 2 occur on the same tour: ScR R1,R2 ;

• Risks R 1 and R 3 occur on the same tour: ScR R1,R3 ; • Risks R 2 and R 3 occur on the same tour: ScR R2,R3 ;

• Risks R 1 and R 2 and R 3 occur on the same tour: ScR R1,R2,R3 . Each senario have a probability of occurence describe in the (Table VI). This probability is calcul with different value of all potenciel risk.

To simplify indicators calcul we propose in first time to considere that the time allowed is the cost, the driving time by the vehicle V 2 is the environmental impact and the driving time is the indicator of health and sécurity. On these first results,

S 2 V 1 W 1 Sit 1 Sit 3 Sit 4 Sit 2 W 1 V 2 - - - - - - S 18 V 1 W 1 Sit 3 Sit 4 Sit 2 W 1 - V 2 W 1 Sit 1 W 1 - - - S 40 V 1 W 1 Sit 3 Sit 4 W 1 - - V 2 W 1 Sit 1 Sit 2 W 1 - -
the planner have a lot of indicator to select few scenario as presented in the (Figure 1) step 3. To illustrate the robustess we chose three solutions S 2 , S 18 and S 40 (Table VII). S 2 is at first glance a good solution on the environmental indicartor because the planning is realized also by a electrical vehicle. S 40 is better than the solution S 2 on the health and security indicator because the repartion of the delivery is balanced. S 18 is maybe chose for find a compromise between all solution. Table V shows for each scenario, the same verification of the compliance constraints can thus be performed, as shown by the columns "validity". Finally, the aggregated criteria of measure of the probability to respect constraints is calculated as shown in Table VIII. The solution S 2 , previously identified as the best solution on the environemental level to be highly sensitive to disturbances: 82% of chance to not respect the constraints often the constraint doesn't respected was the autonomy of the vehicle but the time allowed at the tour is also exceeded. For solutions S 18 and S 40 they are not valide on few senario only cause by the autonomy of vehicle. The same solution with vehicle V 1 and V 2 swap will be a robust solution.

The robustness is an other indicator to help the planner to chose an appropriated solution. It's not because a solution is considerate less robust than an other that it is not possible to chose it. The Planner need to realize a lot of compromise and the robustness has to aim to help him in this chose. The robustness shows the solution unstable but this shows also that the first vision of indicator can be unstable and this view should be change the chose of the planner.

V. CONCLUSION

In the planning process, planners are generally facing with lack of indicators that are used to detremine the influnece of the risk on choosing the solution. In addition in these problems, there is a complexity which originates from the impact of risks on the genration of the solutions while trying to satisfy constraints and modify indicators of the first solution. In this study, we proposed a model and a generic approach to help the planners with better evaluation and comparasion of the obtained oriented solutions. For this purpose, we propose a set of indicators that are used to characterize each of the solutions. In a complementary manner, and to reflect the dynamic and stochastic character of execution environment of tours, we present an indicator that determines the robustness of the solutions. This indicator expresses the ability of the solution to respect the framework of constraints despite potential disruptions. We plan to continue the definition of indicators by adding precision, particularly using different ways of weighting the various solutions outside of framework. For example, a solution that does not respect one of the identified constraints will have a different weight of another solution by not observing none constraints.

Figure 2 .

 2 Figure 2. example of process

Figure 3 .

 3 Figure 3. Illustration of Sc in and Scout sets

Hour 8 :

 8 00 8:10 8:40 8:50 9:00 9:25 9:40 9:45 9:55 10:00 10:20 10:25 11:10 11:20 S1 00:10 0:30 0:10 0:10 0:25 0:15 0:05 0:10 0:05 0:20 0:05 0:45 0

Figure 4 .

 4 Figure 4. Example of planning for the solution S1

Table I EXAMPLES

 I OF GENERIC RISKS DEPENDING ON THE TYPE OF TASKS

	Task type	Example of generic risks	Possible cause(s)
	Activity on site Accident / injury of the driver	Missing operator in charge of handling
			Handling equipment not available
		Mistake	Client is absent
			Handling equipment not available
		Accident / injury of the driver	Missing operator in charge of handling
			Handling equipment not available
	Driving	Driver vigilance drop	Loading/Unloading realized by the driver more important than esti-
			mated
			Long driving period
		Insufficient autonomy of the vehicle, especially for
		vehicles with green power electric only, gas only).