Wei Su

Rui Liu

Ting Liu

Jianxi Huang

Xiaodong Zhang

Junming Liu

The estimation of tree height based on LiDAR data and QuickBird imagery

Keywords: estimation of tree height, LiDAR, watershed segmentation, Weibull function

The estimation of tree height is advanced following the development of LiDAR technique. The estimation model of tree height considering suppressed trees is developed in order to extract tree height accurately using LiDAR data. Filtered LiDAR data and Quickbird imagery are segmented using watershed segmentation method based on mathematical morphology to get the boundary of trees. And the highest point in each canopy object is used to estimate tree height. Weibull distribution is used to estimate height distribution of the suppressed trees. The experiment results indicate that the watershed segmentation method based on mathematical morphology is an effective method to extract the boundary of trees. And the R 2 between the tree height estimated using estimation model of tree height considering suppressed trees and the tree height measured by field work is 0.93.

Introduction

LiDAR (Light Detection And Ranging) is a remote sensing technology that can be used to extract detailed 3D structure information. LiDAR is an active system which emits high frequency bursts of laser pulses. So they are able to penetrate relatively dense vegetation canopies, and can be used to get the 3D structure information of trees even for suppressed trees.

There are many researchers explore how to estimate tree height using LiDAR data in recent years, most of them do this using CHM (Canopy Height Model). This CHM is the difference between DSM (Digital Surface Model) and DTM (Digital Terrain Model), which is the height of non-ground objects such as trees, buildings, vehicles etc. Liu et al. develop smoothed CHM to estimate tree height. Firstly, they smooth CHM using Gaussian smoothing method. Then, they develop a double tangent crown edge recognition (DTCER) algorithm to detect the canopy boundary, and the local maximum search algorithm with fixed or variable window was used to recognize crown tops. Persson et al. create DCM (Digital Canopy Model) from subtracting DTM from DSM, then which is smoothed with different scales. A parabolic surface was fitted to the elevation data to determine which scale to choose for different parts of the image, and the height and crown diameter are estimated for the identified trees [START_REF] Persson | Detecting and measuring individual trees using an airborne laser scanner[END_REF] . Popescu et al. use regression models and cross-validation to estimate tree height based on LiDAR data. They found that estimating tree height for deciduous plots gave superior results without calibrating the search window size based on forest type [START_REF] Sorin | Seeing the trees in the forest: using lidar and multispectral data fusion with local filtering and variable window size for estimating tree height[END_REF] . Zimble et al. search potential canopy top points using a moving radius circular search window with 0.9m radius, and adjoining pixels that are higher than 85% of their neighbors are combined into clumps that are assumed to include the peaks of tree crowns [START_REF] Daniell | Characterizing vertical forest structure using smallfootprint airborne LiDAR[END_REF] . A common ground of all of these researches is that there should be some prior knowledge about their study area which is the basis of determining initial parameters used in model. Therefore, there are some difficulties in extending those methods to other areas.

We develop an estimation model of tree height considering suppressed trees in this study. The watershed segmentation algorithm is used to segment tree canopy objects, the highest point in each tree canopy object is searched using a moving window size, and the height of suppressed trees is estimated using Weibull function. The rest of the paper is organised as follows: Section 2 describes the study area and data source. Section 3 presents the methodology of the estimation model of tree height considering suppressed trees. And the rresults, discussion and conclusions are presented in Section 4 and Section 5.

Study area and data sources

The study area is an urban area (Kuala Lumpur City Centre, from 101º 40'46" to 101º 43'25"E, 3º 8'00" to 3º 10'39"N) located midway along the west coast of Peninsular Malaysia, at the confluence of the Klang and Gombak rivers. This region is characterized by rolling topography and flourishing vegetation, and the height of high trees range from 20m to 30m. A subset area covered with dense high trees in Kuala Lumpur City Centre is selected as study area in this study. And the QuickBird imagery and CHM image produced using LiDAR data of this subset are as Figure 1

Image smoothing

Subtracting DEM from DSM, we can get CHM of canopy. However, for the irregularity of LiDAR data and the twigs of the canopy surface, there will be false tree vertex if we use unsmoothed CHM image directly to identify the tree vertex and crown boundary. So we smooth the CHM image in order to remove these noise before identifying tree vertex and crown boundary [START_REF] Tomas | Detection and analysis of individual leaf-off tree crowns in small footprint, high sampling density lidar data from the eastern deciduous forest in North America[END_REF] . However, this smoothing will weaken the difference between tree vertex and around pixels, which will decrease the ability to identify the tree vertices. So he Gaussian filter method is used to smooth CHM image. The filtering kernel of Gaussian filter is as follows:

2 2 x + y -( ) 1 2 g(x,y, )= e 2 2 2    (1)
Where, x and y are the distant to the center of smoothing nuclear, σ is used to adjust the parameters of the Gaussian function. Figure 2 shows a schematic diagram of a 3×3 filtering kernel of Gaussian smoothing algorithm.

g(-1,1,σ) g(0,1,σ) g(1,1,σ) g(-1,0,σ) g(-1,-1,σ) g(0,0,σ) g(1,0,σ) g(0,-1,σ) g(1,-1,σ)

Fig. 3. 3×3 filtering kernel of Gaussian smoothing algorithm

Moving this filtering kernel, we smooth CHM image. For the pixel p(x,y) located in the centre of filtering kernel, its DN value is computed as:

d / 2-0.5 d / 2-0.5 i -d / 2 0.5 j -d / 2 0.5 p( X ,Y ) p( X i,Y j ) g( x i, y j, )             
(2) X and Y are coordinates of the pixel in the image, d is the size of the filtering template which is an odd number such as 3, 5, 7 etc., i and j are the distant to the studied pixel.

Watershed segmentation based on mathematical morphology

Principle of watershed segmentation

The watershed segmentation algorithm is a robust image segmentation tool, which can be used to extract interested objects in the image effectively. Vincent's simulate immersion algorithm is used in this study [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF] . The pixels in gradient image are classified as four types during the watershed segmentation: marked pixel, waterline pixel, unlabelled pixel and classified pixel. The watershed is finished when all pixels are classified pixels or waterline pixels. There are two steps in Vincent's method: sorting and immersion.

(1) Sorting: ascending the pixel value in line with the DN value of pixel, we can get a sorted pixel matrix. In the process of gradually submerged, we can call the pixels which need to be processed. But we needn't process all pixels in every loop, which can speed up the calculation.

(2) Immersion: this immersion is performed in ascending order to call each pixel based on sorted gradient image. Starting from the minimum DN value in the whole gradient image, immersion is done one pixel by one pixel. This loop queue is used to extend labeled catchment basins. Based on the distance to the neighborhood, we label finished watershed.

Fig.4. Flowchart of Vincent's segmentation algorithm

The study area is a tropical rain forest area, where there are dense and high tropical trees. This results in a difficulty in canopy segmentation because there is overlapping between one canopy and another. Watershed segmentation method is an effective way to extract canopy objects because which can take advantage of the height difference between trees fully. In addition, the watershed segmentation method is a method without any parameters, which is suitable in many situations. So we can get an ideal image segmentation result without prior knowledge. However, there will be a minor gradient change if there are two trees which is located near and the canopy top is relatively flat. In this case, it is difficult to identify these tree objects using watershed segmentation based on gradient. Fortunately, watershed based on the distribution of height is used to do further segmentation. The watershed segmentation algorithm based on the distribution of height is as follows: searching the maximum points in each objects produced from previous segmentation using fixed window, the further segmentation will be done if there are two or more maximum points in one objects.

Waterline pixel in the watershed segmentation is the intersection of two or more tree objects. The classification of these points is done by the inclination angle between this point and crown vertex which the studied point is possible belonged to. And the studied point should be classified as the tree objects which have a larger inclination angle. For an example in Figure 4, there is an intersection between tree A and tree B. r1 and r2 are the inclination angles between this intersection and tree A and tree B. And this point should be classified as tree A, otherwise, classified as tree B.

Merging of over-segmentation

There is over-segmentation in watershed segmentation sometimes, especially when there are larger branches highlighted in studied objects. So we should merge these over-segmented image objects. The Full Lambda Schedule algorithm of Robinson et al. is used in this study [START_REF] Robinson | Implementation of a fast algorithm for segmenting SAR imagery[END_REF] , which is as follows:  is the length of the public side of the image object i and the image object j. , ij t is computed using Formula (8), we will merge image object i and image object j when , ij t is less than threshold.

length i j 2 i , j i j i j i j oo t u u / ( ( o ,o )) oo       (8)

Extracting tree vertex

The segmented results using watershed image segmentation method is only the candidate tree objects, which are not necessarily the real tree objects. So we must identify if or not every classified tree objects is tree in the next step. And the NDVI produced from the spectral bands of Quickbird image is used to do this judgment process. The mean NDVI is used to classify the tree objects. The objects with mean NDVI value less than the threshold will be classified as non-tree, which will be deleted; otherwise, which should be classified as tree.

The highest point in an object is classified as tree vertex after this classification, these are used to estimate tree height. We use a 5×5 moving window to calculate the tree height if there more than on maximum value in one image object, and the mean DN value in this 5×5 window acts as the tree height in this case. If the maximum DN value is still equivalent in 5×5 window, we will extend window size until we can get the unique maximum DN value.

Estimating the height of suppressed trees using Weibull function

The height of suppressed trees can't be estimated using those existing methods [START_REF] Maltamoa | Estimation of timber volume and stem density based on scanning laser altimetry and expected tree size distribution functions[END_REF] , including the method described in this Section 3.3. So we should develop some method to do this. Weibull function is proposed by Weibull, who is the Swedish physicist, which is used widely in forestry measuring and research for the flexibility of its form. Weibull function is LiDAR usually only can find the trees in the dominant, and the distribution of the underlying wood is difficult to obtain [START_REF] Maltamoa | Estimation of timber volume and stem density based on scanning laser altimetry and expected tree size distribution functions[END_REF] , to make use of the Weibull function for modeling and prediction and extract the height distribution of the lower wood.The Weibull function is proposed by the Swedish physicist Weibull, because is widely used in forestry research, it is the main method for the prediction of the height distribution and the diameter distribution [START_REF] Meng | A study of relation between D and H-Distribution by using the Weibull funtion[END_REF] . The three-parameter Weibull function of the density function f (x) and distribution function F(x):

c 1 c c x a x a ( ) e( ) x a f ( x ) b b b 0 x a             (9) c xa ( ) b F( x ) 1 e  



(10) a is a location parameter, b is a scaling function, c is a shape parameter. In the application of using the Weibull function to simulate the trees high, a is the height of minimum spanning tree height, so that a becomes a set value .At this point, the three-parameter Weibull function is converted to two-parameter Weibull function used here are a set worthy of the two-parameter Weibull function. Because the determination of the laser radar, lidar accuracy combined with mutual shadowing between trees, and stacked cover, trees and undergrowth of wood is difficult to be identified, measured by the tree is not up to the minimum spanning tree, so only is the measured number of large trees, trees number is unknown, the corresponding can not be used to complete the Weibull function to the left truncated Weibull function to simulate the distribution of trees, get the value of b and c, then a complete weibull function simulate the distribution of the trees. Left truncated Weibull function of the density function f (x) as follows:

c 1 c c c x t x ( ) e(( ) ( ) ) x t f ( x ) b b b b 0 x t            (11)
t is the measured height of the laser radar. Weibull function parameter estimation There are several ways, the most common and accurate maximum likelihood estimates the parameters of the weibull function, 11 b and c, values, and then substituted into equation (10) to get a high degree of less than t probability F (t).

Results and analysis

Using the algorithms discussed in section 2 to segment LIDAR data of the two pilot areas, and the results as shown below. We can be seen from Figure 6 that trees are basically of two pilot areas have been extracted, and the border in line with the actual situation, some spread between the artificial construction of the trees are extracted. However, due to the spectrum appear synonym spectrum with the spectrum of foreign bodies, and some buildings were classified as trees, some trees are not extracted.

NDVI can remove a large number of non-vegetation information. Can be seen from Figure 5, trees and buildings are well separated. But can be found in Figure 5 (a): the top right of the red marked man-made structures has also been judged as trees, the average NDVI of these areas is 0.4, but the bottom left of the green marked regional NDVI trees region NDVI only 0.36. Staggered addition of artificial buildings and trees, mixed distribution, but when the building height to a certain extent, the shadow of the building will be obscured to some trees, then trees NDVI values also fall, leading to an average of some trees image object NDVI dropped to less than the threshold, thus false positives for non-trees objects. Artificially built-up areas in the second image, the trees did not get to some of the buildings shadow, but in the forest part of the overall NDVI values are relatively high, so the shadow of the forest part of the less affected. In order to retain more information of the crown, set low NDVI, the first image of NDVI threshold is relatively high with 0.35, while the second image of 0.3.

Image segmentation and image objects of the merger after the completion of the judgment, the process to complete the Full Lambda the Schedule method, this method taking into account the spectral information of the Quickbird image object, area information and public side length of the adjacent area. Spectrum and size of the difference is proportional to the object with two images, side length and area with the adjacent edges of the two image objects is inversely proportional to the differences. This method is more reasonable to merge some of the region has a long common edge, but because of the high spectral similarity of the two image trees region, when the trees were over-segmentation, some of this should be a whole tree is divided into multiple regions ,but adjacent region area does not necessarily have similar spectral information, if it belongs to the edge of the area from the whole partition, it would be more similar to the edge of the area and other adjacent tree spectrum rather than it should be subordinate the tree itself. And because of the generally over-segmentation regions of various shapes, so the shape parameter and the height distribution of the information is difficult to join .So we use spectral information, but the merger of the two adjacent regions, the spectral information is too close to the merger, only two adjacent regional spectral information over a certain range when the merger. In this paper, spectral information is also used NDVI values.

Extracted from the tree height and the measured points were compared, the results are as follows: We can see in Figure 7 that using LiDAR measured tree height and measured height has a high correlation, which is close to the 1:1 ratio. Measured tree height = 1.12 × the LiDAR tree height +0.087, R2 = 0.93.

Simulated trees using Weibull function, and predict the two images to detect the elevation distribution of the trees greater than 5m tree. 4057 and 1284, respectively, on the simulation results of the trees greater than 5m R20.96 and 0.90, predicted less than 5m of trees 26 and 12. As the first image is mainly forest areas ,so using the tree height distribution measured in the first image as a test ,set to detect minimum tree height is 15m, left-truncated Weibull function simulation to obtain between 5-15mthe tree height distribution, the results show that the 3702 tree greater than 15m 5-15m of the tree 443, the distribution of more than 5m tree obtained by the Weibull function and to detect the distribution of greater than 5m tree compared to the R2 is 0.93. The above three Weibull function simulation through the Kolmogorov -Smirnov test can account for that using left truncation Weibull function can be effective simulation and forecasting detected and undetected trees.

Conclusions

This paper combined the use of LiDAR data with Quickbird images and constructed the model of finding tree height which take into account the height of the trees of lower wood .This method relative to the maximum search approach, be able to very good use of the spatial variation of the height distribution, rational use of trees and between trees existing gap in the watershed segmentation process, bifurcation identification of trees overlap or tree good effect; After the first watershed segmentation using watershed segmentation on a high level and some smooth crown can be split effectively, and using the Full Lambda the Schedule algorithm can effectively merge the over-segmentation area, so not only guarantee the adequacy of the segmentation to reduce over-segmentation; Combination of LiDAR data and Quickbird images can effectively distinguish between forest and other surface features ;Using the Weibull function to simulate the distribution of tree height and predict the distribution of the trees. This method of height information extraction can obtain high extraction efficiency in the region of dense forest and artificial buildings and trees interlaced distribution. Proposed in this study to consider the underlying Measured tree height Estimated tree height wood to a high degree of tree height to find the height model extraction process does not require too much prior knowledge, application scalability. Of course, the proposed algorithm there is need for further improvements: 1. In the process of tree height extraction, the threshold of NDVI and canopy size need to manually set and the threshold size would affect the results, the NDVI impact on the distinction between forest and other surface features, and crown size affects the crown dense region segmentation results. 2. In the process of canopy object mergering, each region is relatively small and are independent, the use of the consistency of every region of the spectrum to the merger will bring some errors. If you can use the multi-echo or echo image of the laser radar, canopy and lower operating directly on point clouds, merged correctly.

  (a) and (b) respectively, which are resample to 1m spatial resolution.

Fig. 1 . 3 MethodologyFig. 2 .

 132 Fig.1. Study area (subset of Kuala Lumpur City Centre, Malaysia). (a) Quickbird image, (b) CHM image

Fig. 5 .

 5 Fig.5. Schematic diagram of ownership of the intersection point

o

  are the area of image object i and j respectively, i u and j u are the mean DN value of image objects i and j, distance between the image object i and the image object j, length i j ( ( o ,o ))

Fig. 6 .

 6 Fig.6. Segmentation results ((a) segmentation result of test area one (b) segmentation result of test area two)

Fig. 7 .

 7 Fig.7. tree height fitting results

Acknowledgements

This research is funded by Agricultural Research Elitist Fund of Ministry of Agriculture of the People's Republic of China, Open-end Fund of Key Laboratory of Agri-informatics, Ministry of Agriculture, P. R. China [2012003], and National Key Technology R & D Program of the Ministry of Science and Technology [2012BAD20B0103]. We would like to thank the reviewers for valuable comments and constructive advice.


R Mezhoud

F.-Z Ighezou

R J Lombard

Connecting the ground state mean square radius to the dipole excitation

Keywords: 03.65Ge, 31.10 +z, 32.70Cs 1

We study relationships between the dipole excitation and the ground state ms radius of a two-body system in the case of local potentials. We recall the inequality obtained long ago by Bertlmann and Martin, and discuss correction factors transforming the inequality in an approximate expression. Connecting the correction factor to the contribution of the lowest dipole state to the sum rule, we get a lower bound to the ms radius. Inverting the relationships yields a bound for the square of the dipole transition matrix element, and thus a bound to the lowest dipole state transition rate.

1 Introduction.

The Bertlmann-Martin inequality [1,2] yields an upper bound to the ground state mean square radius (ms radius) from the lowest dipole excitation energy. It concerns a particle moving in a local potential or two particles interacting via a scalar interaction. This relationship is exact in the case of the harmonic oscillator. A simple correction factor proposed by Bertlmann and Martin [1] extends its exact character to the Coulomb potential. It also constitutes a good approximation for a wide class of potentials [3,4]. First established in D = 3 dimensional space, this situation is valid in any D ≥ 2 space [5].

The correction factor proposed by Bertlmann and Martin [1] uses only spectroscopic data, i.e. the energy of the 2s state. In this respect it is very efficient. At the same time, it is legitimate to ask the question of a possible alternative, based on another observable. Because the saturation of the inequality is directly connected to the rate of the sum rule exhausted by the lowest excited state, linking the correction factor to this rate is a very natural choice. Thus the purpose of the present work is investigate what can be learned from this rate.

As we shall see, the present study establishes a lower bound to the ground state ms radius in terms of the dipole transition matrix element. Its inversion produces a bound to the transition probability or the lifetime of the lowest dipole state. The correction factor built from the lowest contribution to the dipole sum rule is found to be depending on the potential, but in such a way that it can still be used to improve the lower bound mentioned above.

The dipole excitation case.

We consider the Schrödinger equation in the D = 3 dimensional space, for a particle in a local potential, or equivalently two spin-less particles interacting via a scalar interaction. The spherical symmetry is also assumed. In this case, the Bertlmann-Martin bound [1,2] reads

< r 2 > ≤ 3 2(E 1,1 -E 1,0 ) . (1)
Here, < r 2 > is the ground state ms radius. Use is made of units h = m = 1. The states are labeled by n and , where n -1 is the number of nodes of the wave function, and the angular momentum. The equality is reached for the harmonic oscillator.

The bound (1) is obtained from an under-estimation of the dipole sum rule :

[E 1,1 -E 1,0 ] ∞ n=1 | < n, 1|z|1, 0 > | 2 (2) ≤ ∞ n=1 [E n,1 -E 1,0 ]| < n, 1|z|1, 0 > | 2 = 1 2 .
Here z = 2 π/3 r Y 1,0 (θ, ϕ) [START_REF] Edmonds | Angular momentum in quantum mechanics[END_REF]. The next step uses the closure ∞ n=1 |n, >< n, | = 1 in the left hand term of (2), together with < 1, 0|z|1, 0 >= 0.

Formally, the bound (1) can be transformed into an equality by introducing a correction factor

< r 2 >= 3 2(E 1,1 -E 1,0 ) ϕ (3)
Clearly, ϕ depends on the potential, and nothing is gained unless this factor can be reached in a model independent way. Bertlmann and Martin [1,2] have proposed the empirical expression :

ϕ ≈ 1 - 1 4 E 2,0 + E 1,0 -2E 1,1 E 2,0 -E 1,0 2 . ( 4
)
This form is exact for the Coulomb potential and preserves the equality in the harmonic oscillator case. It is an excellent approximation for a wide class of potentials [3], and in particular for the power-law potentials [4]. These ones are defined by

V (r) = sign(α) V 0 r α ; α > -2 and = 0 , (5)
Note that the correction factor (4) loses its efficiency as α < -1 and tends to -2.

The following general conclusions have been drawn from a large ensemble of tested potentials :

-the factor (4) constitutes a good approximation if the number of bound states is large or if the ground and first dipole states are well bound. The accuracy which is reached for < r 2 > is of the order of few percent or better.

-the evaluation of ϕ through (4)becomes questionable if the number of bound states is limited or if the ground and first dipole state are weakly bound. This situation occurs for short range potentials. We briefly recall that these potentials decay fast enough to be considered zero beyond a given radius. They are characterized by the existence of critical coupling constants, depending on , below which no bound state exists. It is clear, for instance, that the expression (4) loses its meaning if the coupling constant is so low that the lowest dipole level is not bound.

As stated in the introduction, the inequality (1) is saturated for the harmonic oscillator, i.e. it becomes an equality when the first excited state exhausts the sum rule. Consequently it is intuitively legitimate to search for an expression of ϕ based on the fraction of the sum rule carried away by the lowest excited state. It is the purpose of the present work to investigate this question, and to check how far such an expression is dependent on the potential.

Let us denote by

T = 2[E 1,1 -E 1,0 ]| < 1, 1|z|1, 0 > | 2 , (6)
the contribution of the lowest dipole excited state to the sum rule.

A first general property implies T ≤ ϕ, which is proved in the following way. From Eqs (3) and ( 6), we obtain the following ratio

T ϕ = 4π | < 1, 1|rY 1,0 |1, 0 > | 2 < r 2 > . (7)
The ground state ms radius can be written

∞ 0 R 2 1,0 (r)r 4 dr = ∞ o R 2 1,0 (r)r 4 dr Ω (Y 1,0 ) 2 dΩ (8) = 4π R 2 1,0 (r)(Y 1,0 ) 2 (Y 0,0 ) 2 r 4 drdΩ = 4π < 1, 0|r 2 Y 2 1,0 |1, 0 >
Here, R n, (r) is the radial part of the wave function Ψ n, ( r) = R n, (r)Y ,m (θ, φ). Then by using the closure, we get

< r 2 > = 4π n | < n, |rY 1,0 |1, 0 > | 2 (9) ≥ 4π| < 1, 1|rY 1,0 |1, 0 > | 2 .
It immediately proves that

T ≤ ϕ . ( 10
)
As before, the equality is reached for the harmonic oscillator. Because of this property, we get a strict lower bound for the ms radius. It reads

< r 2 > ≥ 3 | < 1, 1|z|1, 0 > | 2 . ( 11
)
As illustrative example, this can be checked with the Coulomb potential for a charge unity. In length units a 0 = 1, together with the Bertlmann-Martin bound, it gives

1.665 < < r 2 > = 3 < 4 .
To go beyond this result, it is necessary to calculate explicitly ϕ and T for few potentials.

We have started with the power-law potentials (5). Due to scaling properties, the results are independent on the coupling constant V 0 . The calculations have been performed for -1.5 ≤ α ≤ 2. This ensemble has been chosen as a typical set of long range potentials. Note that the case α = 0 corresponds to the ln r potential [START_REF] Hall | [END_REF]. Apart from the Coulomb potential and the harmonic oscillator, the Schrödinger equation is solved numerically. The results are displayed in fig 1, where the variations of T and ϕ are shown versus α. The two quantities are not proportional, they obey T ≤ ϕ, the equality being reached asymptotically as α → 2.

We have completed our study by considering the following typical short range potentials :

-Hulthén :

-V 0 e -r 1-e -r -Yukawa : -V 0 e -r r -polynomial : -V 0 1 (1+r) 6 4 -Pöschl-Teller : -V 0 1 cosh 2 (r) -Gaussian : -V 0 e -r 2
This sample yields an overview of decays at large r and behaviors near the origin. In the case of the Pöschl-Teller potential, we have added the D=1 case, which is solvable analytically.

In order to exhibit the connection between T and ϕ, we found convenient to express the results under the form

ϕ = T ν . ( 12
)
The exponent ν is displayed in fig 2 as a function of T for the above sample of potentials. In short range cases, the existence of a critical coupling below which the lowest = 1 level is not bound makes the calculations near the critical value very delicate. For this reason the curves start at values for which we are confident that the numerical accuracy is of the order of 0.1 %.

The spreading of the results is large, and clearly the connection between ϕ and T ν is depending on the potential. However, this dependence has a well defined structure. The different curves outline a hierarchy among the short range potentials according to the behavior at the origin and the decay at large r. A 1/r divergence at r = 0 yields results close to the ones of the powerlaw potentials, as exemplified by the the Hulthén and the Yukawa potentials. The situation is sensibly different for potential finite at r = 0. The discrepancy with respect to the power-law case is larger the faster the potential decays.

It may seem paradoxical that the shift of the results goes toward the large T as the decay of the potential becomes faster. Actually, it is due to the fact that for a bound state to exist requires the wave function to be more concentrated inside the potential when the decay is sharper. This produces a larger overlap between the 1s and the 1p wave function for rapidly than slowly decaying potentials.

The curve corresponding to the power-law potentials draws a limit. Denoting by ν pl the exponent corresponding to the power-law potentials, we obtain the following lower bound

ϕ ≥ T ν pl , (13)
which yields

< r 2 > ≥ 3 2(E 1,1 -E 1,0 ) T ν pl . (14)
Whereas the actual ν may differ from ν pl by a factor 2 or 3, the impact on ϕ does not exceed 10 -15 %. Consequently, this lower bound is useful, and it can be combined with the Bertlmann-Martin limit (1). We note also that for power-law potential, T ν overcomes the difficulty encountered with the correction factor (4) for -2 < α < -1.

The results for Pöschl-Teller in D = 1 and D = 3 are very close. Thus the number of spatial dimensions does not play an important role as long as D remains low. To be more specific, this conclusion is expected to remain valid for the three-body case in the Jacobi coordinates, which is equivalent to D = 6. Because of its analytical solutions, the D = 1 case offers a possibility to establish the limits T → 0, 1 on firm grounds. To proceed, we replace r by x in the above definition of the Pöschl-Teller potential, and set V 0 = s(s + 1), which is more convenient to express the analytical solutions. The potential strength is thus a function of s.

The states are labeled by a single quantum number n, the ground state wave function is given by Ψ 0 (x) = B(1/2, s) -1/2 cosh -s (x), where B(1/2, s) is the Beta function ; the spectrum is given by E n = -(s -n) 2 , which means that s > 1 in order to bind the lowest dipole state.. We get [8, 9, 10]

< x 2 >= 1 2 ζ(2, s) ≤ 1 2(E 1 -E 0 ) = 1 2s -1 (15) ϕ = 2s -1 2 ζ(2, s) (16)
T = s -1 2 (2s -1) Γ(s -1/2) Γ(s) 4 . ( 17
)
Here, ζ and Γ are the Riemann Zeta function and the Gamma function [11], respectively.

From ( 4), ( 14) and ( 15), it can be checked that ϕ and T → 1 as s → ∞. In the limit s → 1, which corresponds to the critical coupling for the n = 1 state, we found from ( 14) and (15)

T → s -1 2 π 2 ; ϕ → π 2 12 . (18)
To quote numbers, for s = 1.001 we get ϕ = 0.8225, T = 0.0049 and ν = 0.0367. It illustrates the fact that ν undergoes a vertical tangent as T → 0, as it can be checked from lim s→1 dν/dT ≈ -0.195/[(s -1)log 2 (s -1)].

3 A bound to the dipole state lifetime.

The inequality (11) can be inverted, and the measurement of the ground state ms radius provides a upper bound to the dipole transition matrix element. Corollary, it yields a lower bound to the lifetime of the dipole state.

If the correction factor is known, we have

| < 1, 1|z|1, 0 > | 2 ≤ 1 2(E 1,1 -E 1,0 ) ϕ . (19)
This is still a strict upper bound. On the other hand, to the extent that ϕ is well approximated by ( 4), it is tempting to write

| < 1, 1|z|1, 0 > | 2 ≤ 1 2(E 1,1 -E 1,0 ) 1 - 1 4 E 2,0 + E 1,0 -2E 1,1 E 2,0 -E 1,0 2 . ( 20
)
To state whether the lower bound character of this relation is ensured requires to know if (4) is actually an upper bound to the true ϕ. From systematic studies [3,4], we found this situation to occur in a large class of potentials, and to be slightly violated in the case of weakly bound systems. This violation is however very small, not exceeding a few percents, and cannot overbalance the inequality (19). Thus we conjecture that (20) is still a upper bound, although it cannot be proved mathematically.

It is interesting to note that (19) yields an upper bound to the transition probability and a lower bound to the life-time of the dipole state in terms of the dipole excitation energy and the energy of the second = 0 state. If we consider the decay of the dipole state by photon emission, the transition rate is given by [8]

R = 4 3 α 3 (E 1,1 -E 1,0 ) 3 he 4 | < 1, 1|z|1, o > | 2 , (21)
with α the fine structure constant. By using (19) re-establishing the h and m dependence, it leads to

R ≤ 2 3 α 3 (E 1,1 -E 1,0 ) 2 h me 4 ϕ (22)
and inversely

τ ≥ 3 2α 3 me 4 (E 1,1 -E 1,0 ) 2 h 1 ϕ . ( 23
)
The equality is reached for the harmonic oscillator, as it can be verified easily. The Coulomb case offers an obvious possibility of testing ( 23). Taking into account E 1,1 -E 1,0 = 3me 4 8h 2 , α = 1/137.0037 and h3 me 4 = 2.42 10 -17 sec, we get τ ≥ 0.910 -9 sec , where the actual value is 1.610 -9 sec [START_REF] Flügge | Practical Quantum Mechanics[END_REF].

Conclusions.

We have studied the relationships between the ground state ms radius and the lowest dipole excitation. The problem addresses one particle in a local potential or two particles interacting via a scalar interaction, in D = 3 space, assuming spherical symmetry. Starting from the Bertlmann-Martin inequality, we have discussed two kinds of correction factors transforming the inequality into approximate or exact expressions. The first one, due to Bertlmann and Martin, requires only the knowledge of the 2s state energy. It is very efficient for a wide class of potentials.

The second one is based on the sum rule rate exhausted by the lowest dipole excited state T . It is expressed as a power-law : ϕ = T ν . This relationship is found dependent on the potential. The results derived for the power-law potentials constitute a limit. In this respect, by measuring T and using this limit, a lower bound to the ms radius is obtained, efficient for the short range potentials.

The present results show also that the simultaneous measurements of the ground state ms radius, T and the energy of the 1p and 2s excited states contain a lot of information concerning the potential. They fix characteristics of the potential at least at a qualitative level.

Finally, the inequality (3) expresses a lower bound to the ground state ms radius in terms of the dipole transition matrix element squared. Inverting this relationship produces an upper limit to this square matrix element. Combined with the perturbative treatment of the one photon emission, it produces bounds to the transition probability and the lifetime of the lowest dipole state. To the extend that the Eq. ( 4) yields a sound approximation to ϕ, these bounds require only the knowledge of spectroscopic data.

The present work concerns one-and two-body systems assuming a local potential and spherical symmetry. Its extension to non-local potentials, three-body systems and/or relaxing spherical symmetry deserves attention. However, for this extended cases, the exact solutions are scarce, which makes the investigations much more tedious. Nevertheless they are under consideration for future work.

[11] I.S. Gradshteyn and I.M. Ryzhik, Table of integrals series and products, Academic Press, New York and London, 4th edition (1965).

Captions for figures Fig 1

The correction factor ϕ (dashed line) and T , Eq(6) (solid line) are plotted against α for the power-law potentials (Eq (15)).

Fig 2

The exponent ν is plotted versus T (see Eq (12)) for the sample of potentials studied in the present work. The solid line is the power-law potentials. Besides potentials quoted in the insert with their symbols, the full triangles are the polynomial, and the full circles the Gaussian potentials, respectively. Note that for T ≥ 0.5, the Hulthén and Yukawa results coincide with the power-law values.

Acknowledgments

This work is part of the Franco-Algerian cooperation project CMEP 02 MDU 551.


Florent Altché

Arnaud De

La Fortelle

Analysis of Optimal Solutions to Robot Coordination Problems to Improve Autonomous Intersection Management Policies

The deployment of Cooperative Intelligent Transportation Systems (C-ITS) raises the question of future traffic management systems, which will be operating with an increasing amount of information and control over the infrastructure and the vehicles. This topic of research shares some similarities with robot coordination problems, inspiring our research on autonomous intersection management. In this article, we use a mixed-integer linear programming formulation for timeoptimal robots coordination along specified paths and apply it to intersection management for autonomous vehicles. Our formulation allows to simultaneously solve a discrete optimal vehicle ordering problem, and a (discretized) continuous optimal velocity planning problem taking into account kinodynamics constraints. This allows faster pruning of the decision tree for the discrete problem, thus reducing computation time. A possible application for ITS is to evaluate the efficiency loss from a given vehicle ordering policy, or dynamically adapt policies to improve their efficiency. Moreover, any intermediary solution found by the solver can be used as a heuristically good policy, with proved bounds on sub-optimality.

I. INTRODUCTION

Cooperative Intelligent Transportation Systems (C-ITS) and particularly autonomous vehicles will likely have important implications for future traffic management systems. With advanced communication systems and increased control on both the vehicles and the infrastructure, it becomes imaginable to "perfectly" coordinate vehicular traffic flows. Although this is certainly overoptimistic, as there are many hypotheses to refine and a lot of engineering to do, this thought gives a direction for our research. Provided we have reasonably good information sharing and control, how can we design the best intersection management system for autonomous vehicles? Should we keep traffic lights? Should we only allow vehicles to enter the intersection at full speed to minimize occupation time? Several methods have been proposed [START_REF] Dresner | Multiagent traffic management: a reservation-based intersection control mechanism[END_REF]- [START_REF] Colombo | Least restrictive supervisors for intersection collision avoidance: A scheduling approach[END_REF], as well as several metrics to evaluate such systems [START_REF] De La Fortelle | Autonomous driving at intersections: combining theoretical and practical approaches[END_REF]. However, finding the optimal solution is ultimately NP-hard so that it is difficult to evaluate methods.

In this paper, we focus on the analysis of such optimal solutions. We measure optimality in terms of traffic flow, by minimizing the average sojourn time of the vehicles in the intersection, which is equivalent to maximizing intersection throughput. We present a simplified model combining two important and difficult properties: the optimality of the 1 MINES ParisTech, PSL Research University, Centre for robotics, 60 Bd St Michel 75006 Paris, France [florent.altche, arnaud.de_la_fortelle] @mines-paristech.fr discrete choice of a vehicle ordering (with relative priorities between conflicting vehicles, which encodes a homotopy class of trajectories), and the optimality of individual vehicle trajectories respecting this ordering, which is a continuous problem. Most previous work mostly considered the optimization of the continuous component, leaving the discrete optimization up to exhaustive enumeration. Our objective is to show that treating both problems simultaneously allows to find good heuristics for near-optimal priority-assignment policies.

This paper is articulated as follows. In Section II, we present a review of related work in the field of autonomous intersection management and multi-robot coordination. In Section III, we briefly describe the mixed-integer linear programming framework allowing to solve the time-optimal coordination problem for multiple moving robots with kinodynamic constraints. In Section IV, we illustrate on a simple example how this framework can help evaluate intersection management policies. In Section V, we show how nearoptimal policies could be designed for real-time applications. Finally, Section VI gives some perspectives.

II. RELATED WORK

Traffic lights are the most widespread method used to regulate heavy-traffic intersection; in this setting, right-ofway is alternated between lanes to avoid possible collisions between vehicles. A good timing of the phases is critical to the performance of the regulated intersection, and a lot of work has focused on finding good or optimal timings [6]- [START_REF] Yin | Robust optimal traffic signal timing[END_REF] off-line. However, such systems are suboptimal by nature, as they do not adapt to real-time changes in the volume of traffic. More recently, several authors have proposed using sensors [START_REF] Mckenney | Distributed and adaptive traffic signal control within a realistic traffic simulation[END_REF], [START_REF] Faye | A distributed algorithm for adaptive traffic lights control[END_REF] or communication [START_REF] -N. Pandit | Adaptive Traffic Signal ControlWith Vehicular Ad Hoc Networks[END_REF] to adapt traffic lights phases and timing in real-time. Other authors [START_REF] Gregoire | Back-pressure traffic signal control with unknown routing rates[END_REF], [START_REF] Gregoire | Capacity-aware backpressure traffic signal control[END_REF] have proposed methods based on queuing theory to demonstrate asymptotic properties for a certain class of traffic control algorithms. Although all of these systems are shown to improve intersections throughput by allowing a finer control of traffic, the need of using phases during which vehicles in a lane are indiscriminately allowed to pass does not guarantee optimality, and it is difficult to measure this loss of efficiency.

Naumann et al. [START_REF] Naumann | Managing autonomous vehicles at intersections[END_REF] are the first to propose using the foreseen abilities of autonomous vehicles to improve intersection management a step further. They introduce a distributed technique based on a reservation mechanism, where vehicles reserve an area of the intersection for a certain duration, and only the vehicles having a valid reservation can enter at a given time. This approach has been perfected by Dresner et al. [START_REF] Dresner | Multiagent traffic management: a reservation-based intersection control mechanism[END_REF] to accommodate a larger number of vehicles. In such systems, the order in which reservations are granted can obviously have a huge influence on the efficiency of the autonomous intersection. In [START_REF] Dresner | Multiagent traffic management: a reservation-based intersection control mechanism[END_REF], reservations are granted in a simple first-come, first-served (FCFS) policy, which is known [START_REF] Vasirani | Evaluating Policies for Reservationbased Intersection Control[END_REF] to be less efficient than pre-timed traffic lights for higher levels of traffic.

In [START_REF] De | Analysis of reservation algorithms for cooperative planning at intersections[END_REF], a reservation-based intersection management method is proposed, along with a heuristically efficient reservation-granting method. Reference [START_REF] Bill | A Next-Generation Intersection Control Algorithm for Autonomous Vehicles[END_REF] further studies improvements to this method. These techniques, however, are based on a set of general rules, such as having the fastest vehicle go first, and give no guarantee regarding optimality.

We argue that the problem of autonomous intersection management is closely related to that of (optimal) multirobot coordination, where the robots are the autonomous vehicles, and that the study of the latter can help find improved methods for the former. Peng et al. [START_REF] Peng | Coordinating Multiple Robots with Kinodynamic Constraints Along Specified Paths[END_REF] are among the first to apply optimization techniques to solve a multirobot coordination problem, taking into account kinodynamic constraints. In their approach, robots evolve along fixed paths divided in segments of predetermined length, and collision avoidance is ensured by requiring that two robots cannot simultaneously be in conflicting segments. The coordination problem is formulated as a mixed-integer nonlinear programming problem, which is very difficult to solve. Using approximations of the vehicle dynamics, the exact problem is translated into two easier mixed-integer linear problems. Digani et al. [START_REF] Digani | A Quadratic Programming approach for coordinating multi-AGV systems[END_REF] use a minimum-velocity hypothesis to formulate a (non-convex) quadratic programming problem for multi-robot coordination with kinematic constraints. A shared issue with those two approaches is the difficulty of handling the dynamic constraints, which need to be either approximated or ignored. The foundation of this paper is a new formulation for the multi-robot coordination problem along fixed paths using the concept of priorities [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF], which has already been used for intersection management [START_REF] Qian | Decentralized model predictive control for smooth coordination of automated vehicles at intersection[END_REF]. Using time discretization, our method allows a resolutionexact modeling of second-order dynamics, presented in [START_REF] Altché | Time-optimal Coordination of Mobile Robots along Specified Paths[END_REF].

III. MILP-BASED OPTIMAL ROBOTS COORDINATION

In this section, we briefly present our mixed-integer linear programming (MILP) formulation for multi-robots coordination under kinodynamic constraints. A more detailed presentation can be found in [START_REF] Altché | Time-optimal Coordination of Mobile Robots along Specified Paths[END_REF].

We consider a set of N robots evolving on predetermined paths inside a bounded coordination region. For automated driving where coordination is mostly needed in the middle of road intersections, the coordination region would be chosen as the center of the intersection including a portion of the roads leading to, and exiting from this center, as illustrated in Fig. 1.

We assume that each robot i follows a predetermined path γ i inside the coordination region, and we only consider the longitudinal behavior of the robots. The dynamics of robot

Coordination region

Fig. 1. Example of paths inside and outside the coordination region for a three-lanes road intersection.

i are described as a double integrator:

ṡi =v i , (1)
vi =a i , (2)
where s i the curvilinear position of robot i along its path, v i its longitudinal velocity and a i its acceleration. This fixedpaths hypothesis is well-suited for road intersections, as vehicles generally stay inside circulation lanes. The origin of s i is chosen so that s i = 0 when the front of the robot enters the coordination region, and s i = s out i > 0 when it fully exits. The velocity is non-negative bounded such that v i ∈ [0, v i ], and we assume that the longitudinal acceleration a i of robot i is bounded to an interval a i , a i with a i < 0 < a i .

We assume that robot i enters the coordination region at time t in i with speed

v in i ∈ [0, v i ],
and is required to leave the coordination region with speed v out i . We let t out i denote the corresponding exit time. Note that v out i should be properly chosen to avoid collisions outside of the coordination region.

For a pair of distinct robots i and j, we call collision set between i and j, noted C ij ⊂ [0, S out i ] × 0, S out j the set of positions (s i , s j ) inside the coordination region where i and j would collide; we say that i and j are conflicting if C ij = ∅. To simplify this presentation, we assume that collision sets are connected, which amounts to assuming that two robots can cross paths at most once. Moreover, we approximate the exact collision sets with minimally bounding polygons, with edges either parallel to the coordinate axes or to the first diagonal s i = s j .

A key remark to be made is that, for two conflicting robots i and j, one necessarily clears a conflict area before the other can enter. Previous work [START_REF] Gregoire | Priority-based coordination of mobile robots[END_REF] has shown that the relative order of passage between conflicting robots can be interpreted as priorities, and that the given of priority relations for all conflicting robots define discrete classes of solutions to the coordination problem. In general, there can exist up to 2 N (N -1)/2 such classes.

Building from these results, we show in [START_REF] Altché | Time-optimal Coordination of Mobile Robots along Specified Paths[END_REF] that comply-ing with assigned priorities is equivalent to respecting a set of linear constraints involving binary variables. Indeed, under our hypotheses, four possible types of conflicts can exist: robots can follow one another throughout the coordination region (a), or have crossing (b), merging (c) or diverging (d) paths, as illustrated in Fig. 2. In these diagrams, priority "i goes before j" (noted i j) corresponds to cases where the trajectory lies on the right-hand side of the collision set C ij . Note that the additional integer variables are needed to deal with the fact that the boundary of these sets are piecewiselinear rather than linear. Moreover, we show that using temporal discretization allows an exact linear modeling of second-order kinodynamic constraints for constant accelerations, unlike the previous approaches relying on spatial discretization which require approximating robots dynamics. More specifically, with a time step duration of τ and noting s k i and v k i the position and speed of vehicle i at step k, constant acceleration secondorder dynamics are enforced as:

s i s j C ij (a) Following paths s i s j C ij (b) Crossing paths s i s j C ij (c) Merging paths s i s j C ij (
s k+1 i -s k i - 1 2 v k+1 i + v k i = 0 (3)
v k+1 i -v k i ≤ a i τ (4)
v k+1 i -v k i ≥ a i τ (5)
Since time is not directly a variable in this formulation, we also introduce a binary variable at each time step and for each vehicle, σ k i , which indicates if robot i has exited the coordination region at step k. To minimize the average sojourn time of robots, we simply can maximize the sum i,k σ k i . Interestingly, this is equivalent to maximizing the throughput of the intersection.

This formulation, which can be handled by many widely available solvers, allows to simultaneously solve two problems: the discrete problem of choosing the best priority assignment, which is the aim of an automated intersection management policy, and the continuous (but discretized) problem of finding optimal trajectories respecting those priorities under kinodynamic constraints.

IV. POLICY EVALUATION

In this section, we show how this mixed-integer linear programming formulation can be used to evaluate an intersection management policy on the example of autonomous vehicles at an intersection. The simulation is based on the free traffic modeling tool SUMO [START_REF] Krajzewicz | Recent development and applications of SUMO -Simulation of Urban MObility[END_REF] and uses its path generation algorithm to compute collisions sets. Vehicles are added using Poisson arrival times and normally-distributed entry speeds (truncated to a minimum speed, and a maximum speed). The optimization problem is solved using the commercial MILP solver Gurobi [START_REF] Gurobi Optimization | Gurobi optimizer reference manual[END_REF], using its Python interface to generate the constraints. Lastly, if the problem is feasible, the solution trajectories are simulated in SUMO using the TraCI interface to verify that they do not generate collisions. Vehicles are modeled as rectangles of 5 m length by 2 m width, with [a i , a i ] = [-3 m s -2 , 4 m s -2 ]. The exit speed for all vehicle i is set as v out i = v = 15 m s -1 and v in i is normally distributed with average 12 m s -1 and standard deviation 3 m s -1 , truncated to [10 m s -1 , 15 m s -1 ].

Simulations are performed on a personal computer running on an Intel Core i7-4790 CPU (clocked at 3.60 GHz) with 16 GB of RAM.

We present the results from three distinct sets of simulations; our objective is to evaluate the performance loss caused by the choice of sub-optimal priority assignments. In order to do so, we compare the optimal assignment found using the problem of Section III with a first-come, first-served (FCFS) ordering of the vehicles as proposed, e.g., in [START_REF] Dresner | Multiagent traffic management: a reservation-based intersection control mechanism[END_REF] which allows simple computation of optimal trajectories.

A. Microscopic simulation

The first simulation compares the optimal trajectories with those resulting from first-come, first-served (FCFS) ordering. Figure 3 shows these trajectories, for a given set of 57 vehicles on a single-lane, four-roads intersection. For illustration purposes, vehicles only travel from south to north, and from east to west in this example. It can be seen that FCFS ordering causes a higher average sojourn time for vehicles, and in fact saturates the intersection, as can be seen in the increasing time vehicles remain stopped. By contrast, an optimal ordering allows the same vehicles to pass without stopping. Note that this does not guarantee that every vehicle is better served by an optimal ordering. For instance, the vehicle with trajectory highlighted in red exits the intersection later in the optimal ordering than in the FCFS ordering.

The attached video submission 1 shows two simulations run on the intersection presented in Fig. 1, with optimal and first-come, first-served ordering.

B. Number of vehicles

It is known that first-come, first-served policies are often efficient for lower volumes of traffic, whereas they lose efficiency when the number of vehicles increases [START_REF] Vasirani | Evaluating Policies for Reservationbased Intersection Control[END_REF]. Fig. 4 1 Also available at https://youtu.be/U7I7x09Hkeo shows the influence of the number of vehicles N on the average relative time loss for the optimal priority assignment and FCFS ordering. This indicator is computed as

t policy i -t alone i t alone i ,
where t alone i is the time vehicle i would have needed to exit the coordination region if there were no other vehicle, and t policy i the sojourn time of the vehicle in the given policy. These data have been computed over a set of 20 initial configurations on the intersection presented in Fig. 1.

C. Influence of time step duration

The discretization time step has a double influence on the solution: first, the safety constraints force vehicles to have passed the conflict area one time step before another vehicle is admitted. Second, we assume constant acceleration during one time step. As the duration of a time step decreases, the average normalized delay is expected to decrease for the same initial configuration.

To measure the effect of the choice of time step duration on solution quality, we show the average relative time loss for various values of the time step, for both optimal and first-come, first-served (FCFS) ordering in Fig. 5. These data have been computed over a set of 20 initial configurations comprising 15 vehicles. The shaded areas correspond to the [50%, 75%] percentiles. Note that, due to the relatively small number of vehicles in the simulations, the first quartile (corresponding to the first 25% of the vehicles) have a nearzero time loss for both optimal and FCFS orderings. It can be seen that the average delay in the optimal ordering is always below that of FCFS ordering. This result emphasizes the fact that optimizing vehicles trajectory for a given ordering is less effective than choosing a good ordering with lowerquality trajectories. In the remainder of this article, a time step duration of 1 s is used.

V. DESIGN OF NEAR-OPTIMAL POLICIES

In this prospective section, we describe how using the presented formulation could allow to design near-optimal policies for real-time applications. A first possibility for reducing computation time with bounded optimality loss is to use less strict termination criteria for the MILP solver. Indeed, branch-and-bound algorithms used in most solvers keep exploring the decision tree until the current solution is very close to the best known bound for the objective function, usually with tolerated relative difference (or gap) of less than 10 -4 . In our particular formulation, we observe that the solver generally finds relatively quickly a reasonably good solution, then spends a lot of time improving this solution to gain a few percents of optimality. For real-time applications, a possible approach to decrease computation time is to increase the maximum tolerated gap, leading to a potentially sub-optimal solution in a much shorter time while ensuring an upper bound on the loss of optimality. In Fig. 6 (red and blue curves), we compare the average computation times for a gap tolerance of 10 -4 and 10 -1 over 10 instances for a time step duration of 1 s; in our test instances, the actual average optimality loss from increasing gap tolerance to 10% is approximately 5%. Using this method, up to 14 vehicles can be treated in less than a second, whereas only 10 vehicles can be treated in the same time span with a gap tolerance of 10 -4 .

Another possibility for real-time applications is to use a receding time horizon approach. At each time step, the optimization problem can be solved only considering the vehicles currently inside the coordination region, taking their position and speed at the beginning of the time step as initial conditions. Vehicles inside the coordination region are tasked to adjust their control to comply with the optimal trajectory. This process can then be repeated at the next time step. Note that in this approach, priorities assigned at a given time step can still be modified in the following steps, which should have a very limited impact on optimality. The green curve in Fig. 6 shows the average computation time needed to converge (with a 10% gap tolerance and 1 s time steps), as a function of the number of vehicles simultaneously present inside the coordination region. Using this approach, up to 24 simultaneous vehicles can be treated in less than a second, which is reasonable for real world applications. For an even higher number of vehicles, choosing a longer time step duration could also be a possibility, since the proposed formulation ensures that such solutions remain feasible.

Note that vehicles are not required to follow the resolutionoptimal trajectory found by the solver, as long as they respect their assigned relative priorities. For reasonable control laws with time as optimization criterion, the above-computed near-optimal policies will likely remain efficient although no optimality guarantee could be given. Notably, [START_REF] Qian | Decentralized model predictive control for smooth coordination of automated vehicles at intersection[END_REF] has demonstrated a decentralized model predictive control scheme allowing vehicles to cross an intersection under assigned priorities.

VI. CONCLUSION AND PERSPECTIVES

In this paper, we have used a simplified autonomous intersection model: there are no unexpected events, no pedestrian, no human-driven vehicles, perfect communication, etc. How-ever, this simplified model allows to reduce the complexity -although the problem remains NP-hard -and we can get optimal solutions for small instances of traffic flows. This is a first lesson: very clearly, such a time-optimal algorithm can be used to manage simple intersections where the traffic flow is limited. Using the rather standard average normalized delays criteria, this allows a significant decrease in the waiting times (from up to 40% to around 10%) comparing with a First-Come First-Served policy where all vehicles also optimize their speed profile. We must emphasize this is true for fluid traffic patterns. However, it shows that theoretically, one can compute an optimal scheduling in some cases. Moreover, we know from previous work that our priority framework allows easily to take into account unexpected events [START_REF] Gregoire | Priority-based intersection management with kinodynamic constraints[END_REF] or even distributed schemes [START_REF] Qian | Decentralized model predictive control for smooth coordination of automated vehicles at intersection[END_REF], which could allow to relax some hypotheses.

What we need to do now it to get a better idea of the optimal patterns. We know from previous studies that FCFS is not optimal for dense traffic situations, where gated or exhaustive policies (in the sense of queuing theory, similar to adaptive traffic lights) are better, and may be optimal. This study is another confirmation FCFS is not optimal, though for very small traffic it is closer to optimality. Very clearly, we must change our methods for intersections subjected to large traffic flows since real-time computation of an optimal solution, even for our simplified model, is out of reach.

Moreover, we also showed that the branch-and-bound algorithm applied to this particular problem converges much faster than it terminates, meaning we can get optimal (or slightly suboptimal) solutions long before we have the proof of its optimality. This may lead to some approximated optimization that could still be operational for large flows. This is clearly an interesting avenue of research, as very few results exist on what a very good policy would be like for an intersection with dense to saturated traffic.

Fig. 2 .

 2 Fig. 2. Possible bounding polygons for the collision set between two conflicting robots i and j. The edges of C ij are straight lines with equations s i = A, s j = B or |s i -s j | = C where A, B, C are constants.

Fig. 3 .

 3 Fig.3. Trajectories of the vehicles in a single-lane intersection for FCFS and optimal ordering of the vehicles. The gray areas correspond to the center of the intersection. Notice that trajectories do not intersect in the gray areas.

Fig. 4 .

 4 Fig.4. Comparison of average normalized delays for optimal and FCFS ordering for an increasing number of robots. The shaded areas correspond to the [50%, 75%] percentiles of vehicles for each policy.

Fig. 5 .

 5 Fig.5. Influence of the time step duration on the average normalized delay for optimal and first-come, first-served ordering. The shaded areas correspond to the [50%, 75%] percentiles of vehicles for each policy.

Fig. 6 .

 6 Fig.6. Computation time, depending on the number of vehicles, for a gap tolerance value of 10 -4 and 10 -1 , and for the time-receding approach with tolerance value of 10 -1 .

École des Ponts ParisTech, Cité Descartes, 6-8 Av Blaise Pascal, 77455 Champs-sur-Marne, FranceThis research was supported by the European FP7 project AutoNet2030 (Grant Agreement NO. 610542).


Konstantin V Avramov

Bifurcations at combination resonance and quasiperiodic vibrations of flexible beams

Keywords: vibrations, flexible beam, Bubnov-Galerkin method, multiple-scale method, bifurcations, quasiperiodic vibrations Introduction

where

F x t F x l t ( , ) ( / ) ( ) = - 0 3 δ
cos Ω is the periodic transverse force, ρ is the density of the material, E is the elastic modulus, A and J are the area and moment of inertia of the cross section, and δ( ) is the delta function.

Introduce the following dimensionless parameters and variables:

ε ρ β ρ µε W W r x x l t EJ Al t l Ar E * * * , , , = = = = 4 2 2 , ε δ δ f l F EJr x l x l 0 3 0 1 3 3 = -= - , ( ) ( ), * (2)
where ε << 1and r is the radius of inertia of the cross section.

Omitting the asterisks in (2), we rearrange [START_REF] Avramov | Forced nonlinear vibrations of flexible beams at combination resonance[END_REF] in terms of the new variables:

W W W W dx W f x t tt xxxx xx x t + = -         + ∫ ε µ 1 2 2 2 0 1 ( , ), (3)
where

f x t f x t ( , ) ( / ) ( ). = - 0 1 3 δ cos Ω Represent the vibrations in the form W t nx n n = ∑ 2 η π ( ) (
). sin Applying the Bubnov-Galerkin method to (3), we obtain the following system of ordinary differential equations:

&& & η ω η ε π η η µη k k k k i i k k k i h t + = - -         + ∑ 2 4 2 2 2 2 2 2 c o sΩ , ω π π k k k h f k k = =       = 2 2 0 2 3 1 2 , , , ,... sin . ( 4
)
Consider a combination resonance: 2 4), we will apply the multiple-scale method [START_REF] Nayfeh | Perturbation Methods[END_REF] based on the following change of variables:

1 2 Ω = + + ω ω εσ. To analyze Eqs. (
η π β ε k k k k a t k t t t O = + + + ( ) ( ( )) ( ) cos cos 2 2 2Λ Ω , (5)
where

Λ k k h k = - 4 4 25 4 4 π ( ) .
As a result, we obtain the system of modulation equations

γ σ π π χ π χ π ′-+ + + - +     5 16 5 4 6 4 2 1 2 2 2 2 0 2 6 0 2 6 2 1 1 2 a a f f a a a a     = cos γ 0, ′ + - = ′ + - = a a f a a a f a 1 1 0 2 1 6 2 2 2 0 2 1 6 1 4 0 0 µ χ π γ µ χ π γ sin sin , , (6)
where . . m m [START_REF] Agalovyan | Solution asymptotics of classical and nonclassical static and dynamic boundary-value problems for thin bodies[END_REF] The transition from (4) to [START_REF] Nayfeh | Perturbation Methods[END_REF] through the change of variables ( 5) is based on rigorous mathematical considerations. This procedure is fully described in [START_REF] Nayfeh | Perturbation Methods[END_REF]. The stability of the periodic motions of (6) corresponds to the stability of the periodic motions of (4). System (6) can be written in terms of the variables ( , , ) ( , , )

γ σ β β χ = -- = ⋅ - T 1 1
x y z a a a = 2 2 1 cos sin γ γ : ′ = ---         + + + - x x f y z x y y f µ σ χ π π χ 6 5 16 4 4 4 0 2 6 2 2 2 2 0 2 ( ) 1 6 π xy z , ′ = -         -+ - + + y f x y f z z x y σ χ π µ χ π π 6 5 16 4 4 0 2 6 0 2 1 6 2 2 2 2 ( )x f x z + 4 0 2 1 6 2 χ π , ′ = -+ z z f y µ χ π 4 0 2 1 6 . ( 8
)
To analyze the nonlinear vibrations of flexible beams, we will represent the solution of Eq. ( 3) in the form

W a t x a t x = + + + + 2 2 4 2 1 2 1 2 2 2 cos sin cos sin ( ) ( ) ( ) ( ) ( π β π π β π ψ x t O , ) ( ) + ε , ψ π π ( , ) ( ) [( ) ] ( ) x t f t m x m m m = - + + - = 4 3 1 3 1 4 3 1 25 0 4 4 0 cos sin Ω ∞ = ∞ ∑ ∑     + - + + -     ( ) [( ) ] ( ) . 1 3 2 4 3 2 25 4 0 m m m x m sin π (9)
2. Bifurcations of Fixed Points. It is easy to verify that the dynamic system (6) has a fixed point ( , ) ( , ) a a

1 2
0 0 = and fixed points that satisfy the system of equations

a a f a f 1 2 2 0 2 6 2 1 2 2 0 4 1 2 12 2 4 6 5 8 16 = + -+         = , . µ χ π σ π χ π (10)
The solutions of Eqs. [START_REF] Evensen | Nonlinear vibrations of beams with various boundary conditions[END_REF] have the form

a p A B 1 1 2 2 5 ( , ) = ± π σ , ( 11
)
where p f Equality [START_REF] Nayfeh | Nonlinear Oscillations[END_REF] defines two groups of fixed points, denoted by the superscripts A and B. These points correspond to two groups of equilibrium states of the dynamic system (8):

r x y z r x y z A A A A B B B B = = ( , , ), ( , , ),
where

δ σ χ π = - 6 0 2 6 f , x p f p y f p A B A B , , , = ± ± = ± π χ δ π µ χ δ 5 0 2 1 5 0 2 1 2 10 10 , z p A B , = ± 4 10π δ . ( 12
)
Let the above fixed points be associated with an amplitude surface (see Fig. 1) that represents the dependence a a f To analyze the fixed points for stability, let us determine the eigenvalues λ i of the Jacobian matrix of the vector field [START_REF] Avramov | Bifurcation analysis of a vibropercussion system by the method of amplitude surfaces[END_REF]. The stability of the fixed points r A B , is defined by quantities λ i A B ( , ) as follows:

Fig. 1 S a 1 σ C S A u B F L G 6χ G 1 4 1 χ f 0 2 6 / π f 0 1 2 6
( ) / π C 1 L 1 λ µ λ µ λ µ 1 2 3 1 2 3 2 2 ( , ) , ( ) , ( )
, , ,

A B A B R R = - = -± = -± (13)
where R ap

1 2 2 2 1 2 5 8 , = µ π m .
To describe the stability of the fixed point (a 1 , a 2 ) = (0, 0), we will use the quantities

λ µ λ µ χ π σ χ π 1 2 3 0 4 1 2 12 0 2 6 2 2 4 1 4 6 = - = -± - -         , . , f f ( 14
)
Note that the nonhyperbolic fixed points of the sheet C satisfy the equation λ 2 = 0. Therefore, we have

σ χ π χ π µ -         = - 6 16 4 0 2 6 2 0 4 1 2 12 2 f f . ( 15
)
To analyze Eq. ( 15), we introduce new variables ( ~, ~/ ) σ π f 0 2 6 . Then this equation rearranges to

- +         = ~σ µ π χ µ 2 2 0 2 6 2 1 2 2 4 4 f 1. ( 16
)
Equation ( 16) describes the hyperbole L in Fig. 1. The domains of stable and unstable states of the dynamic system (8) on the amplitude surface are labeled by S and U, respectively. The line L represents supercritical pitchfork bifurcation on the section GC 1 and subcritical bifurcation on the section C 1 B. The codimension-two bifurcation point C 1 separates the lines of supercritical and subcritical bifurcation. The amplitude surface contains a saddle-node bifurcation line L 1 , which connects the fixed points of the sheets A and B.

To analyze the dynamics of system (8) at C 1 , we will take advantage of the central-manifold method [START_REF] Bryuno | Analytical form of differential equations[END_REF][START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF], which allows reducing the dimension of a dynamic system. Note that it is central manifolds on which bifurcation phenomena occur. In the neighborhood of the fixed point x = y = z = 0, Eqs. ( 10) can be represented as

′ ′ ′           = - - + + -        x y z 2 0 0 0 2 0 4 2 1 1 1 1 µ µ ν χ µ χ ν µ µ                + - + +         + x y z y z z δ µ χ ν µ π
) ( ) ( ) δ χ ν µ π              , ( 17
)
where ν π µ χ

1 0 2 6 1 2 = - f .
Note that the bifurcation point C 1 is observed when ν δ 1 = = 0. Write system (17) in terms of the new variables ( , , ) ( , / , / ): ,,),( ,,, ), ( χ ν δ

u v w x y z z y = + - 2 2 ′ = - + ′ = + ′ = - + u u r u v w v v r u v w w 2 2 2 1 1 1 2 µ δ νχ δ µ ( ,
1 1 2 ) ( , , , ), w r u v w - (18)
where

[ ] r w v v v w v w u 1 1 1 1 1 1 2 2 2 2 2 2 5 16 2 = -- + + - + + δ χ ν µ χ ν µ π ( ) ( ) ( )( ) , r v w u v w u 2 2 1 1 1 2 2 2 2 8 2 5 8 2 = - + + - + + δ χ ν µ π ( ) ( ) ( ) . ( 19
) λ δ χ ν µ ν 1 2 1 2 , ( ). = ± +O (31)
The equation ν = 0describes a saddle-node bifurcation line. The behavior of trajectories near the bifurcation point C 1 is shown in Fig. 2 as a bifurcational diagram [START_REF] Nayfeh | Nonlinear Oscillations[END_REF].

3. Quasiperiodic Vibrations of a Flexible Beam. The steady vibrations of a flexible beam that correspond to the fixed points of the dynamic system (8) are defined by

W a t C x a t C x = - + + - + 2 2 2 1 1 1 2 0 cos ( )sin ( ) cos ( )sin ( ) * * ν π ν γ π Ψ ( , ) ( ), x t O + ε (32)
where γ µ π χ

0 6 0 2 1 2 = arcsin ( / ( )) f .
Thus, the vibrations of the beam are represented by the sum of three functions of time with frequencies ν 1 , ν 2 , and Ω , As follows from (34), vibrations (32) occur at two frequencies ν 1 and ν 2 and their combination Ω . It is well known that if the frequencies ν 1 and ν 2 are not related as integers, then the vibrations are quasiperiodic. Minor changes in ν 1 or ν 2 may turn quasiperiodic vibrations into periodic. This phenomenon is called synchronization. Note that based on the results presented in this paper we may demonstrate that vibrations (32) are always quasiperiodic.

  This surface consists of joined sheets marked by A, B, and C in Fig.1. The sheet C corresponds to the fixed point (a 1 , a 2 ) = (0, 0); and the sheets A and B represent fixed points marked by the same letters.

Fig. 2 C

 2 Fig. 2

The central manifold of the dynamic system (18) can be represented as power series with undetermined coefficients. After determining these coefficients by the method described in [START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF], we represent the central manifold in the form

where O(4) denotes terms of orders v v v v The motions of the dynamic system (8) on the central manifold are described by the equation

where a O a O µχ ν δ = . This equation approximates the curve L at C 1 (Fig. 1). The equation

describes the motions of the dynamic system (8) on the central manifolds of the fixed points of the line L. This manifold is stable when δ < 0 and unstable when δ > 0. The motion of system (8) on the central manifold of the point C 1 is described by the equation

It follows from this equation that the point C 1 is stable. Consider saddle-node bifurcation L 1 (Fig. 1) at C 1 . To this end, let us introduce the variables

where 10 2 1 π δ δ = .

Rearranged in terms of the new variables, system (8) has the form

We do not present the expressions of the functions f x y z

( , , ), and f x y z ,),( ,,),( ,, ). ( 27)

The functions F 1 , F 2 , and F 3 are also unwieldy and, thus, are not presented. The central manifold is

where

, ν, ν , and ν 3 . The motion of system (8) on the central manifold is described by the differential equation The stability of the fixed points v 1 2 , is determined by the characteristic numbers λ 1 2 ,