Enhancing sampling-based kinodynamic motion planning for quadrotors

Alexandre Boeuf, Juan Cortés, Rachid Alami and Thierry Siméon Abstract-The overall performance of sampling-based motion planning algorithms strongly depends on the use of suitable sampling and connection strategies, as well as on the accuracy of the distance metric considered to select neighbor states. Defining appropriate strategies and metrics is particularly hard when considering robot dynamics, which is required to treat constrained motion planning problems for quadrotors. This paper presents an accurate but computationally fast quasimetric to determine the proximity of dynamic states of a quadrotor, and an incremental state-space sampling technique to avoid generating local trajectories that violate kinodynamic constraints. Results show that the integration of the proposed techniques in RRT-based and PRM-based algorithms can drastically decrease computing time, up to two orders of magnitude.

I. INTRODUCTION

Sampling-based motion planning algorithms [START_REF] Lavalle | Planning Algorithms[END_REF] explore the connectivity of the configuration space (or C-space) by sampling configurations and attempting to connect them to build an underlying graph whose nodes and edges are feasible configurations and local paths, respectively. The sampling method and the steering method (or interpolation method) used to connect samples are therefore essential components of this type of algorithms. Most of these algorithms also rely on a notion of distance on the C-space to define a connection policy, which is also determinant for the efficiency of the exploration. In the Rapidly-exploring Random Tree (RRT) algorithm [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF] for instance, expansion bias towards the largest Voronoi regions is a direct consequence of the nearest node connection strategy. In Probabilistic Roadmap (PRM) based algorithms [START_REF] Kavraki | Probabilistic roadmaps for path planning in high-dimensional configuration spaces[END_REF], significant efficiency improvements are obtained by prioritizing connections to nodes according to their distance to a sampled configuration. Since the topology of the space being explored is induced by the steering method, a suitable distance metric in the connection policy should be correlated to the actual length of the local paths, rather than using a simple Euclidean distance. Note that discussions on the importance of the distance metric have been recurrent since early work on sampling-based motion planning algorithms [START_REF] Amato | Choosing good distance metrics and local planners for probabilistic roadmap methods[END_REF]. Indeed, this is an open problem that still motivates the development of new approaches [START_REF] Bharatheesha | Distance metric approximation for state-space RRTs using supervised learning[END_REF].

The motion planning problem becomes even harder when dealing with under-actuated systems. The problem difficulty if often circumvented by using a decoupled approach. For instance, motion planning for quadrotors can be treated this way: In a first stage, a collision-free path is planned for the robot's bounding sphere. This path is then transformed into a trajectory that meets both continuity requirements and the differential constraints [START_REF] Richter | Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments[END_REF]. However, for very constrained problems, this approach is insufficient since there is no solution path for the bounding sphere. Several approaches have been proposed to generate local trajectories of quadrotors [START_REF] Mueller | A computationally efficient algorithm for state-to-state quadrocopter trajectory generation and feasibility verification[END_REF], [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF], but few works tackle the global kinodynamic motion planning problem [START_REF] Lavalle | Randomized kinodynamic planning[END_REF], which is mandatory in those cases.

Defining appropriate sampling and connection strategies is indeed more difficult when considering the robot dynamics.

Treating the kinodynamic problem implies to explore the state space (or X -space) rather than the C-space, which means also considering derivatives of the degrees of freedom. Finding the ideal metric in this case may be as hard as solving an optimal control problem.

This paper builds on our recent work on quadrotor motion planning [START_REF] Boeuf | Planning agile motions for quadrotors in constrained environments[END_REF], which presented a steering method to connect kinodynamic states of this type of robots, together with first results on its integration inside a simple RRT algorithm. The principle of the steering method is recalled in Section II. Here, we propose a (quasi-)distance metric and a state-space sampling strategy aimed to significantly improve the performance of motion planning algorithms making use of this type of steering method. Sections III and IV present these two contributions together with some results that validate their relevance. Then, results showing the influence of both the metric and the sampling strategy on the overall performance of RRT-based and PRM-based planners are presented in Section V.

II. STEERING METHOD

This section briefly recalls some notions about the state space of a quadrotor and gives an overview of the steering method presented in [START_REF] Boeuf | Planning agile motions for quadrotors in constrained environments[END_REF] 1 . The steering method uses fourthorder splines to generate smooth trajectories between two arbitrary states of a quadrotor, ensuring continuity up to the jerk (and therefore, up to the angular velocity), which is important to facilitate controllability along the planned trajectories. These trajectories minimize flying time (with respect to the proposed closed-form solution) while respecting kinodynamic constraints.

A. State space of a quadrotor

A quadrotor is a free-flying object in R 3 . It has three degrees of freedom in translation (position of the center of mass: [x, y, z] ∈ R 3 ) and three in rotation (roll, pitch and yaw angles: [θ, φ, ψ]). Thus, its configuration space is C = SE(3). Considering a standard dynamic model, the system is differentially flat for the outputs [x, y, z, ψ] (see [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF]). This implies that any trajectory in the state of the flat outputs is flyable if their derivatives are correctly bounded. We define a kynodynamic state as a vector x = [x, y, z, ψ, ẋ, ẏ, ż, ψ, ẍ, ÿ, z, ψ] ∈ X . With X the state space:

X = P × V × A
where P, V and A are intervals of R 4 . Note that although there are no actual limits (if friction is neglected) on linear velocity of the center of mass of a quadrotor, for safety reasons, it might be interesting to take such limits into consideration for trajectory planning.

B. Steering method

The steering method we consider provides a solution (S, T ) to the problem:

               [S(0) Ṡ(0) S(0)] = x 0 ∈ X [S(T ) Ṡ(T ) S(T )] = x T ∈ X ∀t ∈ [0, T ]        Ṡ(t) ∈ V S(t) ∈ A ... S(t) ∈ J .... S (t) ∈ S (1)
where S is the trajectory in the state of the flat outputs, T is the total time of the motion and V, A, J and S are zero centered intervals of R 4 . Obtaining the time-optimal trajectory is difficult and computationally expensive. This is not well suited for integration in sampling-based motion planners since they make extensive use of the steering method (usually thousands of calls to solve constrained planning problems). Therefore, we propose a method to compute a trajectory that approaches the optimal one with a simpler (imposed) shape enabling a rapid, analytical solution. The method works in two stages. First a solution is computed independently for each output x, y, z and ψ. Each solution is a fourth order spline s(t) with a bang-null snap profile (i.e. a trapezoidal jerk profile) such that [s(0) ṡ(0) s( 0 for each output of index i). These splines are divided into seven main temporal phases (see Fig. 1). Four of them (noted A, C, E and H) correspond to acceleration variations. Phases B and G are constant acceleration phases. The goal is to maximize acceleration variations in order to reach full speed as fast as possible and maintain it as long as possible. Solutions are then synchronized to form the local trajectory whose total duration is T = max i=1..4

)] = [x 0i v 0i a 0i ] and [s(T i ) ṡ(T i ) s(T i )] = [x T i v T i a T i ] with T i ∈ R + (different
TIME (s) POSITION VELOCITY ACCELERATION A B C D E G H a 0 v 0 x 0 a F v F x F a B a G v D
T i . Note that this method can be generalized to various types of systems. For a serial manipulator, for instance, each degree of freedom would be treated as a flat output.

III. QUASI-METRIC IN THE STATE-SPACE

A. Motivations

The efficiency of the state-space exploration using randomized kinodynamic motion planning algorithms relies on a good distance metric. However, as discussed in [START_REF] Lavalle | Randomized kinodynamic planning[END_REF], computing the actual distance between two states is as hard (and thus as costly) as solving the corresponding optimal control problem. Our steering method provides a deterministic suboptimal solution to such a control problem. Therefore, it defines a quasi-metric2 M * SM : (x 0 , x T ) → T on the state space. Because of the dynamics of the system, a trajectory in the state space from x 0 to x T is indeed necessarily different from a trajectory from x T to x 0 , and thus M * SM is not symmetric. Although this steering method is computationally fast, it is still too costly to be used for nearest neighbor search inside a sampling-based planner. This section presents a method to approximate the quasi-metric M * SM at a very low computational cost, and presents results that show its relevance.

B. Approximate quasi-metric

The complexity of the problem (1) defined in Section II is mainly due to its order (four) and to the inequality constraints on the derivatives. We propose to solve a simpler time optimal control problem for the third order (i.e. by considering the jerk as the control input), in one dimension and without constraints other than the bounds on the control input. The problem is then to find for each output of index i the couple (S i , T i ) such that:

       min T i ∈ R + s.t. [S i (0) Ṡi (0) Si (0)] = [x 0 v 0 a 0 ] ∈ R 3 [S i (T i ) Ṡi (T i ) Si (T i )] = [x Ti v Ti a Ti ] ∈ R 3 ∀t ∈ [0, T i ], | ... Si(t)| ≤ j max ∈ R + (2)
Pontryagin maximum principle (see for example [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]) says that the optimal control is necessarily saturated, i.e.:

∀t ∈ [0, T i ], ... Si(t) ∈ {-j max , j max }
with at most two control commutations. Solving (2) implies to find T i and these (at most) two commutation times, which requires to solve polynomial equations of maximum degree four. Further details can be found at http://homepages.laas.fr/aboeuf/iros15annex.pdf. The proposed quasi-metric is then defined as:

M SM : (x 0 , x T ) → max i=1..4 T i

C. Results

Here we present results of an experimental test to validate the proposed approximate quasi-metric. 10 4 pairs of kinodynamic states were randomly sampled in X = [-5, 5] 3 × [-5, 5] 3 × [-10, 10] 3 , considering J = [-20, 20] 3 and S = [-50, 50] 3 . Note that, without loss of generality and for simplification purposes, we consider here a constant yaw. For each pair (x 1 , x 2 ), we computed the value M * SM (x 1 , x 2 ) of the quasi-metric induced by our steering method, the value M SM (x 1 , x 2 ) given by the proposed approximation, and the value ED(x 1 , x 2 ) of the euclidean distance in R 3 considering only the position of the center of mass. We study the distribution of the relative error between M * SM and M SM , i.e. the quantity:

RE M SM (x 1 , x 2 ) = 1 - M SM (x 1 , x 2 ) M * SM (x 1 , x 2 )
For comparison, we also provide the relative error RE ED between M * SM and ED. Fig. 2 shows histograms of the distributions of these errors, Tab. I shows key statistical values of these distributions and Tab. II gives mean CPU times in milliseconds for a single core of an Intel Xeon W3520 processor at 2.67GHz.

The low standard deviation of the distribution of the relative error for the proposed quasi-metric is a measure of the quality of the approximation. These results also provide empirical evidence that M SM and M * SM are equivalent since for all pairs (x 1 , x 2 ),

0.16396 ≤ RE M SM (x 1 , x 2 ) ≤ 0.85540 which implies 1 10 .M * SM (x 1 , x 2 ) < M SM (x 1 , x 2 ) < 10.M * SM (x 1 , x 2 )
This means that M SM and M * SM are inducing the same topology on X which means that the cost-to-go defined by our steering method is correctly evaluated by M SM . This is clearly not the case for ED.

IV. SAMPLING STRATEGY

This section presents an incremental state-space sampling technique that increases the probability of generating connectible states. The definition of such states is first presented together with our motivations. Then the different steps of the method are presented. Finally some results are provided.

A. Connectible states and motivations

Trajectories generated by the steering method presented in Section II do not guarantee to respect bounds on the position of the robot. Such a constraint is typically violated when samples are close to the boundary of the workspace and the velocity is high, so that it is not possible to decelerate to avoid crossing this positional limit. In a similar way, bounds on velocity can also be violated. If acceleration is too high and velocity is close to the limit, produced trajectories will be invalid because velocity can not be reduced in time to meet the constraints. Note however that the imposed shape for the trajectories produced by our steering method guarantees that bounds on acceleration are respected.

We say that a kinodynamic state x ∈ X is forwardconnectible (respectively backward-connectible) if and only if there is a state x ∈ X such that the local path (x, x) (respectively (x, x)) produced by the steering method lies entirely in X (i.e. respects bounds on position, velocity and acceleration). A state that is forward-connectible and backward-connectible is said to be connectible (nonconnectible otherwise). This is illustrated on Fig. 3.

In case of uniform sampling of the state space, nonconnectible states can be generated, and thus, local paths computed to connect those states have to be discarded a posteriori by the planner. This is rather inefficient since generating and testing a local path for validity is a costly operation. The goal of the sampling technique proposed below is to notably reduce the probability of generating nonconnectible states, and hence to improve the performance of planning algorithms.

The sampling technique proceeds in a decoupled and incremental way. First, acceleration is uniformly sampled. The idea is then to compute a set of velocity values for which the state is known to be non-connectible. Velocity is then uniformly sampled outside this set. Finally, given this couple (velocity, acceleration) a set of position values for which the state is known to be non-connectible is computed, and the position is then uniformly sampled outside this set. 4. Acceleration and velocity (blue curves) around t = 0. Saturated acceleration variation is applied in order to determine the limits on initial velocity (red triangles). Red curves are velocities for these limits.

B. Sampling velocity

This subsection explains how to compute the set of velocity values for which the state is known to be nonconnectible, given a uniformly sampled acceleration value a s . Explanations are given for one output. Fig. 4 illustrates this explanation. Let us denote v max and a max as the bounds on the absolute value of velocity and acceleration respectively. We study acceleration a(t) and velocity v 0 (t) on a neighborhood around t = 0 for a(0) = a s and v 0 (0) = 0. The idea is to apply a saturated acceleration variation and determine the extrema of v 0 (t) in this neighborhood. Using them, we can compute the limits on velocity v s such that v(t) = v 0 (t) + v s lies in [-v max , v max ]. We use notations defined in section II. For t > 0, phase A of our steering method is applied. Phase H is applied for t < 0. The sampled value a s locally imposes a direction of variation of v 0 (t) on phases A and H. We want to reverse this direction of variation in minimum time. This is equivalent to driving a(t) to zero in minimum time. For that, we set a B = a G = -sign(a s ).a max . This corresponds to the highest acceleration variation achievable by our steering method. Note that, by construction, acceleration is symmetric during phases A and H (i.e a(-t) = a(t)) and v 0 (t) is anti-symmetric (i.e v 0 (-t) = -v 0 (t)). Since a(t) is a second order spline strictly monotonic on phase A, it is straightforward to compute the unique t 0 > 0 such that a(t 0 ) = 0. The value v bound = v max -|v 0 (t 0 )| is then the upper bound on the absolute value of v s . This means that if |v s | > v bound then v(t) = v 0 (t) + v s will violate the constraints on velocity. A velocity value v s is then uniformly sampled in [-v bound , v bound ].

C. Sampling position

Given a couple (v s , a s ) for one output, this subsection explains how to compute the set of position values for which the state is known to be non-connectible. Fig. 5 illustrates this explanation. The principle is similar to the one in the previous subsection. Velocity v(t) and position x 0 (t) are studied around t = 0 for a(0) = a s , v(0) = v s and x 0 (t) = 0. We apply a saturated velocity variation and determine the extrema of x 0 (t) in this neighborhood. Using them, we can compute the limits on position x s such that x(t) = x 0 (t)+x s lies in [-x max , x max ] (bounds on position). For t > 0, phases A to C of our steering method are applied. Phases E to H are applied for t < 0. We want to reverse the direction of variation of the position imposed by v s as fast as possible. This is equivalent to driving v(t) to zero in minimum time. For that, we set v D = -sign(v s ).v max for both phases A to C and E to H. This corresponds to the highest velocity variation achievable by our steering method. The only difference here is that neither v(t) nor x 0 (t) have symmetry proprieties. We compute t + > 0 such that v(t + ) = 0 and t -< 0 such that v(t -) = 0. If v s ≥ 0 then x + = x max -x 0 (t + ) and x -= -x max -x 0 (t -) else x + = x max -x 0 (t -) and x -= -x max -x 0 (t + ). A position value x s is then uniformly sampled in [x -, x + ].

PHASE E PHASE H PHASE

A PHASE C TIME (s) POSITION VELOCITY vs v ( ) t t + t + t - t - vmax - vmax xmax xmax v D v D x +
x -x 0 ( ) t -Fig. 5. Velocity and position (blue curves) around t = 0. Saturated velocity variation is applied in order to determine the limits on initial positions (red triangles). Red curves are positions for these limits.

D. Results

We provide here some results concerning the sampling strategy. The conducted experiment consisted in testing the validity of local paths computed between uniformly sampled pairs of states in X = [-5, 5] 3 × [-5, 5] 3 × [-10, 10] 3 , with J = [-20, 20] 3 and S = [-50, 50] 3 . Yaw was kept constant. We measured the percentage of valid paths (i.e. which lies entirely in X ) over 10 4 calls. We then repeated this operation using our sampling technique. With uniform sampling only 11.53% of the produced paths were valid. With our sampling technique 95.58% of the paths were valid. Moreover, we can test for a given state if each output is respecting the constraints defined by our sampling technique, i.e. if velocity lies in [-v bound , v bound ] and if position lies in [x -, x + ]. When this is not the case we consider that the state is not connectible. Our sampling technique is obviously generating connectible states every time whereas with uniform sampling about 90% of the generated states are not connectible (with respect to this criteria).

V. INFLUENCE ON PLANNING ALGORITHMS

This section describes variants of RRT and PRM algorithms, and provides some results on the influence of both the proposed quasi-metric and the sampling strategy on them. These algorithms are well known, but had to be adapted to the non-symmetry of the steering method and the associated quasi-metric. Because of this non-symmetry, the underlying graph is directed.

A. Directed bi-RRT

In the well known undirected version of the bi-RRT algorithm [START_REF] Lavalle | Rapidly-exploring random trees: progress and prospects[END_REF], two trees T S and T G are constructed in parallel. T S grows from the start configuration and T G from the goal configuration. Each iteration for one of the trees consists of sampling a configuration q rand , finding its nearest neighbor q near in the tree (according to a defined metric), and extending it toward q rand (using a steering method) to create a new configuration q new . Each time an expansion is successful for one of the trees, a direct connection is attempted between q new and its nearest neighbor in the other tree. The algorithm ends if this local path is valid (i.e. when the trees are connected). In our directed version, both the steering method and the quasi-metric M SM are nonsymmetric, and thus have to be called taking care of the order of the two states. The nearest neighbors N S (x) and N G (x) of a state x in T S and T G respectively are defined as such:

   N S (x) = arg min x∈T S M SM (x, x) N G (x) = arg min x∈T G M SM (x, x)
For an expansion of T S , we test the local path N S (x rand ), x new for validity. In case of success, the algorithm ends if the local path x new , N G (x new ) is valid. For an expansion of T G , the local path x new , N G (x rand ) is tested for validity, and the algorithm ends in case of validity of the local path N S (x new ), x new .

B. Directed PRM

At each iteration of the undirected version of the PRM algorithm [START_REF] Kavraki | Probabilistic roadmaps for path planning in high-dimensional configuration spaces[END_REF], a collision free configuration q is sampled and added to the graph G. For every connected component G i of G, connections are attempted between q and each node of G i in increasing order of distance from q until one is successful. A threshold on this distance can be considered with the aim to reduce computational cost. In our directed version, we consider the strongly connected components G i of G. Moreover, we maintain during the execution the adjacency matrix A G of the transitive closure of the graph of the strongly connected components of G. This square matrix, whose dimension is the number of strongly connected components, is defined by

A G [i][j] = 1 if a path in G exists from every node of G i to every node of G j and A G [i][j] = 0 otherwise. If A G [i][j] = 1 we say that G i is connected to G j . Note that A G [i][j] = A G [j][i] = 1 if and only if i = j.
At each iteration, a valid state x is sampled and added to G (which has n strongly connected components). Its strongly connected component G n+1 = {x} is added to the matrix A G . For every connected component G i of G (i = 1..n), if G i is not connected to G n+1 , connections from every node x j of G i to x are attempted in increasing order of M SM (x j , x) until one is valid. As for the undirected version, a threshold on the value of M SM can be considered here. A G is updated if neecessary. Then, if G n+1 is not connected to G i , connections from x to every node x j of G i are attempted in increasing order of M SM (x, x j ) until one is valid. If used in single query mode, the algorithm ends when the strongly connected component of the initial state is connected to the strongly connected component of the goal state.

C. Results

Results presented below show the influence of the quasimetric and the sampling technique on the two previously presented motion planners. Experiments have been conducted on two different environments shown in Fig. 6 and for the same quadrotor whose diameter is equal to 0.54 meters. We consider V = [-5, 5] 3 , A = [-10, 10] 3 , J = [-20, 20] 3 , S = [-50, 50] 3 (using SI base units). Yaw is kept constant. The first environment, referred to as boxes, is a cube with side length of 10 meters filled with box shaped-obstacles of different sizes. The second environment, referred to as slots, is also a cube with side length of 10 meters but divided in ). Each experiment is designated by an acronym whose meaning is explained in the caption. Results are averaged over 100 runs and are for an implementation in C, integrated in our motion planning software Move3D [START_REF] Siméon | Move3D: a generic platform for path planning[END_REF], and run on a single core of an Intel Xeon W3520 processor at 2.67GHz. Results show a significant improvement of the performance of both algorithms thanks to the integrations of the proposed techniques. However, one can clearly see that the metric and the sampling technique have a more notable effect on one or the other planner. Results for the PRM algorithm shows that the sampling method has a great influence on its performance. Its integration indeed improves CPU time by two orders of magnitude for both environments. On the other hand, one can see that for the slots environment CPU times are slighlty worse with the use of the metric. This can be explained by the difference of computing time between our quasi-metric and the euclidean distance. This is also observed in one case for the RRT algorithm (SRMU vs. SREU). For the RRT algorithm, results show that the influence of the metric is more important. This was to be expected since RRT-based algorithms are known to be very sensitive to the metric. One can see that the number of iterations is significantly reduced, meaning that the search is better guided. The improvement produced by the sampling technique is also very significant for the slots environment but less noticeable for the boxes environment. This can be explained by the fact that, in RRT-based algorithms, the sampled states are not tested for connections but used to define a direction for extension. A new state is then generated according to that direction. One can see that, for the boxes environment, about half of these states are not connectible regardless of the sampling method. Finally note that flying times are given for the raw, non-smoothed trajectories, which explains the rather large difference of path quality between PRM and RRT results.

VI. CONCLUSION

We have presented two techniques to enhance the performance of kinodynamic motion planning algorithms applied to quadrotors: 1) a quasi-metric in the state space that accurately estimates the actual quasi-distance induced by a steering method, but which is two orders of magnitude less computationally expensive; 2) An incremental state-space sampling technique that notably increases the probability of generating connectible states. Results have show that the integration of these techniques can drastically change the performances of PRM-based and RRT-based planners, particularly when solving very constrained problems.

We are currently working on the comparison of the trajectories generated by the proposed steering method and those obtained by a numerical optimal control method. First results tend to show the good quality of our approximate solution, which is orders of magnitude faster to compute. We are also working on the integration of the planning methods in an indoor testbed for experimental validation.

Fig. 1 .

 1 Fig. 1. Illustration of the seven main temporal phases of the spline for one output. Red circled dots are the values to be interpolated. Pink dashed lines are the boundary values for acceleration and velocity.

Fig. 3 .

 3 Fig. 3. Examples of non-connectible states in two dimensions. Red squares are positions and red arrows are velocity vectors. Blue curves are examples of trajectories. Bounds on position are represented in black. The state on the left is not backward-connectible. The state on the right is not forwardconnectible.

Fig. 2 .

 2 Fig. 2. Histograms of the distributions of the relative errors

Fig. 6 .

 6 Fig. 6. Testing environments (a) boxes (b) slots

TABLE III B

 III : BOXES, S: SLOTS, P: PRM, R: RRT, M: PROPOSED METRIC, E: EUCLIDEAN METRIC, I: INCREMENTAL SAMPLING, U: UNIFORM SAMPLING

	Experiment	B.P.M.I.	B.P.E.I.	B.P.M.U.	B.P.E.U.	S.P.M.I.	S.P.E.I.	S.P.M.U.	S.P.E.U.
	CPU time (s)	0.05648	0.07884	3.284	4.409	1.505	1.220	578.5	444.2
	Flying time (s)	8.180	8.772	8.000	8.126	9.074	8.979	8.615	8.387
	Number of nodes	12.11	13.77	78.64	88.33	71.93	61.59	767.9	725.8
	% of not connectible nodes	0	0	82.53	84.38	0	0	89.11	89.19
	Experiment	B.R.M.I. B.R.E.I. B.R.M.U. B.R.E.U. S.R.M.I. S.R.E.I. S.R.M.U.	S.R.E.U.
	CPU time (s)	0.02780	0.04088	0.04144	0.05612	2.165	2.466	558.8	512.9
	Flying time (s)	9.674	10.84	9.365	10.09	25.42	34.72	33.96	55.98
	Number of nodes	8.79	8.84	9.18	10.77	334.5	565.6	4429.4	8813.0
	Number of iterations	26.45	45.04	45.58	65.02	982.9	2502.9	30253.7	196233.3
	% of not connectible nodes	50.34	54.94	51.51	59.85	25.60	34.86	79.41	82.48
	two halves by a series of aligned obstacles separated by 0.40
	meters (hence smaller than the robot diameter). This problem
	is particularly challenging since going across these obstacles
	requires to find a path in a very narrow passage in the state-
	space. Every combination of environment, algorithm, metric
	and sampling strategy has been tested. Results are provided
	in Tab. III for CPU and flying times in seconds, number of
	nodes (and iterations for the RRT) and percentage of not
	connectible nodes (with respect to the criteria defined in
	section IV.D.

The code (C and MATLAB versions) is available upon request.

A quasi-metric has all the properties of a metric, symmetry excepted.

This work has been partially supported by the European Community under Contract ICT 287617 "ARCAS".


Student Member, IEEE Marko Pavlov

Anne Migan

Vincent Bourdin

Michel Pons

Martial Haeffelin

Jordi Badosa

Experimental and Numerical Study of the Influence of String Mismatch on the Yield of PV Modules Augmented By Static Planar Reflectors

Keywords: photovoltaic modules, low concentration, planar reflectors, increase of production, multiphysics model, simulation, optimization

Photovoltaic (PV) modules are generally installed by the application of empirical rules aimed at reducing shadows during the periods of high solar irradiation. A traditional installation on a horizontal surface results in largely spaced rows of modules with a relatively low tilt angle. The addition of inter-row reflectors results in more direct and diffuse flux transmitted to the cells. The "Aleph" (Amélioration de l'Efficacité Photovoltaïque) project aims to define clear rules for optimal settings of systems of PV module rows with fixed inter-row planar reflectors in a given location and under a given climate. Two PV technologies are tested for performance with this type of system: amorphous silicon (a-Si) and polycrystalline silicon (p-Si). This work combines experiments on panel behavior in an outdoor environment on the SIRTA (Site Instrumental de Recherche par Télédétection Atmosphérique) meteorology platform and a multiphysics numerical model used to couple all the important physical phenomena and accurately describe the system behavior. The model includes a ray tracing radiation/optics module based on the Monte-Carlo method, as well as an electrical module simulated in SPICE. This work presents the influence of the string mismatch losses, present at periods of heterogeneous illumination, on the yield of PV modules augmented by static planar reflectors.

I. INTRODUCTION

T HE Aleph experiment is rooted in the recognition of the constraints on the design of conventional photovoltaic (PV) installations on a horizontal plane, such as rooftop installations and large-scale solar farms installed in the field. In a typical installation, system density, defined here as the ratio between the cell surface and the surface of the field, is limited by the need to minimize shading periods which occur in the morning and evening periods. The resulting system is of a relatively low angle of inclination "s" of the PV modules and an inter-row spacing "D" equal to 4-5 times the module height "L" (see Fig. 1). This arrangement results in the interrow space being illuminated during the period of highest solar M. Pavlov, V. Bourdin and M. Pons are with the Laboratoire d'Informatique et de Mécanique pour les Sciences de l'Ingénieur (LIMSI), Orsay Cedex, France (e-mail: marko.pavlov@limsi.fr).

M. Pavlov and A. Migan are with the Laboratoire de Génie électrique et électronique de Paris (GeePs), Gif-sur-Yvette Cedex, France (e-mail: anne.migan-dubois@geeps.centralesupelec.fr).

M. Haeffelin and J. Badosa are with the Laboratoire de Météorologie Dynamique (LMD), Palaiseau Cedex, France (e-mail: jordi.badosa@lmd.polytechnique.fr).

Manuscript received June, 2015.

irradiation (around solar noon). The Aleph project aims to explore the benefit of adding inter-row planar reflectors to increase the system yield. The use of planar reflectors is an already established technique for increasing the solar flux incident on a collector, and various studies have been performed for estimating the increase in irradiation from a given collector-reflector system geometry [START_REF] Ahmad | Comparative study of pv modules with and without a tilted plane reflector[END_REF]- [START_REF] Winston | Nonimaging solar concentrator with uniform irradiance[END_REF]. Ronnelid et al. [START_REF] Ronnelid | Booster reflectors for pv modules in sweden[END_REF] performed a theoretical study of the planar reflector system. For a latitude of 60°N and a specular reflector with a reflectivity of 0.8, a gain in annual output from a-Si PV modules of up to 25% was predicted. Matsushima et al. [START_REF] Matsushima | Concentrating solar module with horizontal reflectors[END_REF] explored the benefit of adding booster reflectors by placing the reflector at a 90°angle with a c-Si PV module and varying the inclination of the assembly. With a latitude of about 35°N (Tokyo, Japan) and a specular reflector with a reflectivity of 0.95, a gain of 50% in produced electricity was estimated. Andrews et al. [START_REF]Photovoltaic System Performance Enhancement With Non-Tracking Planar Concentrators: Experimental Results and BDRF Based Modelling[END_REF] experimentally tested the PV-booster system over a year. For the latitude of 44°N (Kingston, Canada), a specular reflector with a reflectivity of 0.90 and a c-Si module, a gain of 45% of produced energy was observed.

This work presents the influence of the string mismatch losses, present at periods of heterogeneous illumination, on the yield of PV modules equipped by static planar reflectors.

II. EXPERIMENTAL SETUP

The experiment is installed in the SIRTA [START_REF] Haeffelin | Sirta, a ground-based atmospheric observatory for cloud and aerosol research[END_REF] meteorological station in Palaiseau, France, located at 48.71°N, 2.21°E. It currently consists of three rows of six PV modules. Each row is comprised of modules of two technologies: amorphous silicon (a-Si) and polycrystalline silicon (p-Si).

The modules are inclined with an angle of s = 49°(roughly equal to the latitude) and the rows are spaced apart with a distance D/L = 3/1 (see Fig. 1), avoiding most shading in the winter. The two rear rows are equipped with planar reflectors consisting of low-cost mirrors of 3mm thickness. On September 26, 2014, aluminum foil was installed as reflector on row 3. The configuration is adjustable and provides space for an additional row of modules and mirrors. The PV modules are orientated as portrait or landscape to allow for a study of the impact of the module orientation on the degree of string mismatch.

The set is instrumented to acquire data continuously from sunrise to sunset. The I-V characteristic is measured using Access to all atmospheric variables measured by SIRTA has been granted, including global, diffuse, and direct irradiances, local albedo, air temperature and humidity, speed and wind direction at various altitudes, and satellite images.

III. EXPERIMENTAL RESULTS

The collected data demonstrate the interest of this system. Fig. 2 presents the gain in the daily and monthly produced energy achieved by p-Si and a-Si PV modules under concentration compared to the control modules. The modules under concentration represented here are situated in the middle of rows 2 and 3, and are therefore not subject to severe edge effects. The monthly gain is highest at spring and fall, and reaches 26% for a-Si modules and 18% for p-Si modules. Comparable gains were observed for two reflector types, ordinary mirror and aluminum foil. The daily gain reaches 35% at certain periods of the year. Note that higher gains are achieved for a-Si compared to p-Si modules, owing to different layouts of cells.

As all the modules are measured with the same sourcemeter, simultaneous continuous measurements are not feasible and a time lag exists between the determinations of the P mpp of each module. Therefore, measurements made during periods of highly varying solar irradiation (e.g. during the presence of clouds) may suggest unfavourable ratios of energy production.

Unfavourable monthly ratios of energy production in the winter months, as suggested by Fig. 2, can be largely attributed to the above-mentioned data acquisition limitation when taking measurements during periods of intermittent solar resource. This becomes clear when one notes a rather constant gain in daily produced energy during overcast days (days with high daily average DHI/GHI, shown in grey), and a seemingly variable gain during intermittent days (shown in yellow). A secondary factor contributing to unfavourable ratios of energy production in the winter period is the presence of shadows on reflector-equipped modules during early morning and late evening periods of winter days. This shadow does not appear on the frontal modules used as control, resulting in an apparent loss of energy. However, the shadow effect would be identical for the two systems (PV+Reflector and classical installation) constructed with the same inter-row spacing. For these reasons and contrary to what the data suggests, no losses of energy should be expected with the addition of inter-row planar reflectors.

It is interesting to note that the shadow has a different effect on the two technologies, having a greater impact on the performance of p-Si modules. This is because the resulting heterogeneous illumination produces a more severe string mismatch effect for the p-Si compared to the a-Si modules.

IV. NUMERICAL MODELLING

To correctly describe the electricity production of the planar concentrator PV system, it is important to build a multiphysics model integrating all the relevant phenomena, such as the solar resource, the photovoltaic conversion process, as well as the applicable loss mechanisms. Andrews et al. [START_REF] Andrews | Model of Loss Mechanisms for Low Optical Concentration on Solar Photovoltaic Arrays with Planar Reflectors[END_REF] give an excellent overview of the loss mechanisms present in this kind of system, namely angle-of-incidence losses, reflection losses, string mismatch, temperature effects and spectral mismatch.

A. Ray Tracing Module

To describe the total flux transmitted to the cells, which includes the contribution of the reflector and the optical losses, a ray tracing model was developed in the EDStar environment [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF]. The geometry was modelled in 3D using the PBRT framework, commonly used for rendering 3D models. It was described as an inclined reflector surface and an oppositely inclined receiver surface, representing the imaginary PV-cell surface in the plane of the modules encompassing all the modules in the experiment. Reflectance of the reflector as well as the transmittance and reflectance of the PV module front glass were described using the procedure outlined by Hoang [START_REF] Hoang | Coupling optical and thermal models to accurately predict pv panel electricity production[END_REF], applying the refractive indices of materials found in the literature [START_REF] Polyanskiy | Refractive index database[END_REF]. The reflectivity was adjusted to account for the spectral mismatch between the reflected wavelengths and the spectral response of the modules.

To precisely describe the position of the sun relative to the system at any moment, several astronomical variables had to be calculated. Solar declination and the equation of time were defined by the means of empirical formulas as outlined by Bourdin [START_REF] Bourdin | Master physique et ingénierie de l'energie[END_REF], based on the data provided by IMCCE [START_REF] Imcce | Ephémérides générales de position des corps du système solaire[END_REF]. The corresponding maximum error resulting from using the empirical formulas is ±0.20°and 16 seconds, respectively. The local hour angle, solar altitude and solar azimuth, were found using the equations derived by Duffie and Beckman [START_REF] Duffie | Solar engineering of thermal processes[END_REF].

The direct flux was modelled using the "pill-box" sun model, meaning that the surface of the sun, as seen from Earth, was sampled with a uniform probability during the Monte-Carlo simulations. The diffuse flux was modelled using the isotropic sky model [START_REF] Hottel | Performance of flat-plate solar-heat collectors[END_REF].

Monte-Carlo simulations were run with 10,000,000 histories. This yielded a 2D distribution of the factor of concentration or shading, present on the imaginary PV-cell surface (see Fig. 3a). Since the angle-dependent transmission through the module front layers is calculated at each history, this factor includes the corresponding losses. The resolution of this distribution was set at 48x288 "pixels", which is high enough to study the string mismatch effect while resulting in an acceptable standard deviation.

MATLAB was used to transform this 2D distribution into an average concentration by cell, taking as an input the physical arrangement of the modules (see Fig. 3b and3c). values transmitted to the cells. The electrical behaviour of the two kinds of modules were modelled in SPICE, starting with developing an equivalent circuit model of a cell of each PV technology.

B. Electrical Module

A p-Si cell was modelled using a modified 2-diode equivalent circuit model where the parallel resistance was replaced with a voltage-dependent term describing the avalanche breakdown during reverse bias, as described by Bishop [START_REF] Bishop | Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits[END_REF] and Silvestre et al. [START_REF] Silvestre | Study of bypass diodes configuration on pv modules[END_REF]. Figure 4 shows the schematic of the electrical model of a p-Si cell.

Figure 5 shows the current in each component of the p-Si cell model as a function of the applied voltage bias. The cell was parametrized using the values shown in Table I, and subjected to STC conditions (1000W/m 2 , 25°C). The breakdown voltage was set to -15 V, resulting in a rapid increase in the shunt current as the voltage bias approaches -15 V.

An a-Si cell was modelled using a 1-diode model with two additional terms: one to describe the avalanche breakdown (as done for the p-Si cell), and a second to describe the recombination in the i-layer, as outlined by Merten et al. [START_REF] Merten | Improved equivalent circuit and analytical model for amorphous silicon solar cells and modules[END_REF]. Figure 6 shows the schematic of the electrical model of a a-Si cell.

Figure 7 shows the current in each component of the a-Si cell model. The cell was parametrized using the values shown in Table I, and subjected to STC conditions. Similarly to the p-Si cell, the shunt current increases rapidly as the voltage bias approaches the set breakdown voltage, -20 V. Note the

D1 V bias ↑ ↑ I shunt (V j ) I phot ↑ I rec (I phot ,V j ) R s V j V out
I dark I shunt I rec I out Fig. 6. Schematic of the electrical model of a a-Si cell, consisting of a current source, a diode, a series resistance, a recombination current term (as described by Merten et al. [START_REF] Merten | Improved equivalent circuit and analytical model for amorphous silicon solar cells and modules[END_REF]), and a shunt current term (as described by Bishop [START_REF] Bishop | Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits[END_REF]. effect of the i-layer recombination term. Compared to the diode equation, it rises more gradually with voltage.

For each technology, the PV modules were modelled by electrically connecting individual cells in series. In the case of p-Si, two bypass diodes were added, one across each string. These diodes were parametrized as follows: the saturation current and the emission coefficient were set to 1 × 10 -7 A/cm 2 and 1.00, respectively. To illustrate the influence of the bypass diodes, the I-V characteristic was simulated for an east edge p-Si module at various points throughout the morning period, characterized by different distributions of the concentration-shading-transmission factor (see Fig. 8). To emphasize the mismatch effect, the beam normal and diffuse horizontal irradiances were set to 800 and 0 W/m 2 , respectively. The operating temperature was fixed at 70 °C. The simulation was run for August 16, 2014, at: a) 07:55, b) 09:15, c) 10:47, d) 11:08, and e) 11:55 UTC. It can be seen that in forward bias conditions, the influence of the bypass diodes depends on the degree of mismatch, while in reverse bias (around -0.7 V) they result in a rapid increase in current.

The electrical model was calibrated using experimental measurements. The electrical parameters were adjusted through trial and error until the I-V curves measured at several times of the day were reproduced by the model (see Fig. 9 for an example). The obtained electrical parameters are listed in Table I. Figure 9 shows the model behaviour after calibration for an a-Si module under concentration, on May 17, 2014 at 10:39 UTC.

C. Integrated Model Validation

The integrated model was run using the measured module temperatures and atmospheric data of May 17, 2014, a relatively clear day, and the results were compared with module power output measurements. As seen in Fig. 10, the power output of 6 modules was simulated: a central, an edge and a control module of both a-Si and p-Si technologies. The string mismatch effect observable in the performance of edge modules was successfully reproduced, validating the model. On this day, the following gains in produced energy were achieved: 19% and 16% for central and edge a-Si modules, and 10% and 6% for the central and edge p-Si modules. This shows the advantage of the amorphous module topology in portrait orientation, with cells extending along the total height of the module.

D. Example Application

The model was used to evaluate the influence of the distance of a module from the edge of the row on the yield (see Fig. 11). To emphasize the edge effect, the annual production of east edge modules was considered during morning periods (a similar behaviour exists during evening periods for the west edge modules). During these periods, the influence of the reflector varies based on the distance from the edge, favouring the central modules. Furthermore, as the edge effect and the string mismatch are the most influential during clear sky periods, a clear-sky model [START_REF] Lefèvre | Mcclear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF] was used to set irradiance values.

The result of such a study can be used to decide on the optimal length of a reflector overhang to minimise the edge effect and in turn maximise production (in this case, 2-3 module widths). Note that in temperate climates, where the diffuse component presents a significant portion of the total irradiation, the edge effect will be less influential on the yield.

V. CONCLUSION

The experimental data show that the use of fixed planar concentrators is a low-cost technique to increase system yield, resulting in gains in daily produced energy of up to 35% during clear sky days of certain periods of the year. Gains in monthly produced energy of up to 26% and 18% were measured for a-Si and p-Si modules, respectively. A multiphysics model was developed to study the performance of a low-concentration PV system, introducing a new modelling methodology to evaluate the impact of the string mismatch on the performance of a-Si and p-Si modules under various configurations. The string mismatch can have a considerable impact on the performance of reflector-equipped PV modules, and was found to affect p-Si modules more severely than a-Si modules due to the difference in cell layout. The model can be used to find the system setting that will maximize the system yield.

Fig. 1 .

 1 Fig. 1. Schematic of the experiment. PV module rows 2 and 3 are equipped with reflectors, row 1 serves as the control.

Fig. 2 .

 2 Fig. 2. Experimental data: gain in daily and monthly produced energy resulting from concentration. After September 26, 2014, the third row was equipped with aluminum foil as reflector. Values below 1 are due to isolated spikes of high illumination during the measurement of a particular PV module in an otherwise very cloudy day. The color scale represents the daily ratio of diffuse to global irradiation.

Fig. 3 .

 3 Fig. 3. Concentration-shading-transmission factor, in the plane of the modules (a) and average by cell for PV modules in rows 3 and 2 (b),(c) -October 15, 2014, 13:50 UTC

Fig. 4 .

 4 Fig.4. Schematic of the electrical model of a p-Si cell, consisting of a current source, two diodes, a series resistance, and a shunt current term (as described by Bishop[START_REF] Bishop | Computer simulation of the effects of electrical mismatches in photovoltaic cell interconnection circuits[END_REF]).

Fig. 5 .

 5 Fig. 5. Current in each component of the p-Si cell model as a function of applied voltage. The cell was parametrized according to Table I and subjected to STC conditions.

Fig. 7 .

 7 Fig. 7. Current in each component of the a-Si cell model as a function of applied voltage. The cell was parametrized according to Table I and subjected to STC conditions.

Fig. 8 .

 8 Fig. 8. Simulated I-V characteristic of an edge p-Si module, equipped with two bypass diodes, resulting from different distributions of the concentration zone and hence degrees of string mismatch. The bypass diodes are active in mismatch conditions, as well as in reverse bias.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Calibration of the electrical model example -performed by comparing measured and modelled I-V curves. (a-Si module under concentration, May 17, 2014, 10:39 UTC)

Fig. 11 .

 11 Fig.11. Annual morning energy production of a-Si and p-Si modules as a function of distance from the east row edge, evaluated under clear sky conditions[START_REF] Lefèvre | Mcclear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF]. The equivalent energy productions for frontal a-Si and p-Si modules are 183 kWh and 282 kWh, respectively.

TABLE I ELECTRIC

 I PARAMETERS OF a-SI AND p-SI CELLS

	Parameter	p-Si	a-Si
	Area (cm 2 )	22.0	65.9
	J 0 (A/cm 2 )	5.00 × 10 -13	5.00 × 10 -19
	J 02 (A/cm 2 )	2.30 × 10 -7	-
	J SC,ST C (A/cm 2 )	3.80 × 10 -2	8.20 × 10 -2
	Rs (Ω.cm 2 )	0.55	26.36
	Rsh (Ω.cm 2 )	1650	6590
	α Isc (1/°C)	6.00 × 10 -4	1.20 × 10 -3
	n 1	1.00	1.60
	n 2	2.00	-
	a	2.00 × 10 -3	0.10
	m	3.70	3.70
	V br (V )	-15.0	-20.0
	d i (cm)	-	3.50 × 10 -5
	µτ ef f (cm 2 /V )	-	1.60 × 10 -8
	V bi (V )	-	1.75
	I(A)

ACKNOWLEDGMENT

The authors would like to thank the STARWEST team [29], SOLEMS company and LaSIPS for their kind support.


Wenmeng Xiong

José Picheral

Sylvie Marcos

PERFORMANCE ANALYSIS OF MUSIC IN THE PRESENCE OF MODELING ERRORS DUE TO THE SPATIAL DISTRIBUTIONS OF SOURCES

Keywords: array signal processing, distributed sources, model error, performance, MUSIC

In this paper, the direction of arrival (DOA) localization of spatially distributed sources impinging on a sensor array is considered. The performance of the well known MUSIC estimator is studied in the presence of model errors due to angular dispersion of sources. Taking into account the coherently distributed source model proposed in [1], we establish closed-form expressions of the DOA estimation error and mean square error (MSE) due to both the model errors and the effects of a finite number of snapshots. The analytical results are validated by numerical simulations and allow to analyze the performance of MUSIC for coherently distributed sources.

INTRODUCTION

The effects of model errors on DOA estimation has been widely studied in the literature with the sources assumed to be far-field point transmitters or reflectors [START_REF] Weiss | Effects od modeling errors on the resolution threshold of the MUSIC algorithm[END_REF]. In many applications, for instance, aero-acoustic sources imaging [START_REF] Chu | Robust Bayesian super-resolution approach via sparsity enforcing a priori for near-field aeroacoustic source imaging[END_REF], ocean waves parameter estimation [START_REF] Young | Cross-shore decay of cliff top ground motions driven by local ocean swell and infragravity waves[END_REF], bio-medical source localization [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF] or also mobile communication [START_REF] Pedersen | A stochastic model of the temporal and azimuthal dispersion seen at the base station in outdoor propagation environments[END_REF], the physical sources can no longer be considered as point sources and a spatially distributed model of the sources could be more appropriate.

The models for spatially distributed sources have been classified into two types, namely incoherently distributed (ID) sources and coherently distributed (CD) sources. On one hand, for ID sources, signals coming from different points of the same distributed source can be considered uncorrelated [START_REF] Zoubir | Efficient subspace-based estimator for localization of multiple incoherently distributed sources[END_REF] [START_REF] Hu | An ESPRITbased approach for 2-d localization of incoherently distibuted sources in massive MIMO SYSTEMS[END_REF]. On the other hand, in the scenario of CD sources, the received signal components are delayed and scaled replicas from different points of the same one [START_REF] Souden | A two-stage approach to estimate the angles of arrival and the angular spreads of locally scattered sources[END_REF] [START_REF] Zoubir | Efficient DSPE algorithm for estimating the angular parameters of coherently distributed sources[END_REF].

In the literature, different techniques propose to estimate the nominal DOA of distributed sources but they generally require that the shape of the dispersion is known which may not be true in practice. Even if it is the case, a mismatch on the angular spread dispersion between the real signals and the models will bring estimation errors on the DOA estimates.

In previous works, the model errors have been considered, in the case of punctual sources, as uncertainties on the array elements positions or patterns, or inter-element coupling, and are described as random variables [START_REF] Swindlehurst | A performance analysis of subspace method in the presence of model errors : Case of MUSIC[END_REF] [START_REF] Kangas | Finite sample and modeling error effects on ESPRIT and MUSIC diretion estimators[END_REF] [13] [START_REF] Marcos | Les méthodes à haute résolution, traitement d'antenne et analyse spectrale[END_REF]. For ID sources, the performance analyses have been investigated in the presence of the model error as a mismatch between the angular dispersion of the actual sources and the MUSIC estimator in [START_REF] Astély | The effects of local scattering on direction of arrival estimation with music[END_REF], and the ultimate asymptotic performance, that is to say, the lower bound of the estimation covariance has been investigated in [START_REF] Raich | Bearing estimation for a distributed source: Modeling, inherent accuracy limitations and algorithms[END_REF]. We here focus on the CD sources model proposed in [START_REF] Valaee | Parametric localization of distributed sources[END_REF], which is well adapted to the applications such as the aero-acoustic sources imaging [START_REF] Chu | Robust Bayesian super-resolution approach via sparsity enforcing a priori for near-field aeroacoustic source imaging[END_REF]. The MUSIC estimator has been extended for a model with an angular dispersion, and the model error refers to a mismatch between the angular dispersion of the actual sources and that of the extended-MUSIC estimator in our scenario.

In this paper, we first propose a new approach to obtain an accurate approximated analytical expression of the DOA estimation error and MSE. Then, the finite sample effects due to the estimation of the received data covariance matrix will be investigated. Finally, assuming that the shape of the angular distribution is known, an analytical expression of the DOA estimation error as a polynomial function of the angular spread dispersion error is proposed. This expression allows to see explicitly the influence of the model error on the DOA estimation performance, and could be useful in a future work to optimize the antenna parameters in order to reduce the DOA estimation error due to the angular spread of the sources.

The organization of this paper is as follows. The signal model and 1D-DSPE are given in section 2. In section 3, the sensitivity of the estimator is theoretically analyzed. Numerical simulations are presented in section 4 to validate the analytical expressions of the previous section. Finally, conclusions are given is section 5.

SIGNAL MODEL AND 1D ESTIMATOR FOR DISTRIBUTED SOURCES

Let us consider q spatially CD far-field sources impinging on an array of M sensors. The q sources and the M signals received at the array at moment t are denoted by s(t) = [s1(t), . . . , sq(t)] T and y(t) = [y1(t), . . . , yM (t)] T , respectively, with:

y(t) = C(θ)s(t) + n(t), (1)
where n(t) ∈ C M ×1 represents the complex Gaussian distributed additive noise, C(θ) = [c h 1 (θ1), . . . , c hq (θq)] ∈ C M ×q is the array steering matrix composed of q steering vectors c hi (θ) that can be written as proposed in [START_REF] Valaee | Parametric localization of distributed sources[END_REF] by:

c hi (θ) = π 2 -π 2 a(θ + φ)hi(φ)dφ, (2)
where i = 1 . . . q, and a(θ) is the steering vector for a point source which arrives from the DOA θ. In the most general case, the steering vector a(θ) is also a function of the array geometry, the sensor gains, the form of the wavefront, and other possible parameters which are supposed to be known. The function h(φ) is introduced to describe the angular spread distribution (for instance, Uniform and Gaussian distributions) and it can be parameterized by an angular dispersion ∆ which is omitted in the notation. The source signals and the additive noise are considered to be centered Gaussian independent random variables. Assuming that signals and noises are uncorrelated and the sources are uncorrelated with each other, the correlation matrix is given by:

Ry = E[yy H ] = CRsC H + σ 2 b I, (3)
where E[.] is the expectation operator, Rs and σ 2 b are the source covariance matrix and the noise variance, respectively.

The DSPE method proposed in [START_REF] Valaee | Parametric localization of distributed sources[END_REF] is a 2D MUSIC-like estimator allowing to jointly estimate the two parameters θ and ∆ by maximizing the criterion

1 c H h (θ,∆)V b 2
, where V b is the matrix of the M -q eigenvectors associated to the smallest eigenvalue of Ry and which spans the noise subspace. However, the form of h(φ) may not be known in practice. In the following, we will consider that the shape of h(φ) and ∆ are imperfectly known and that the only parameter to estimate is θ. We will then investigate the effect of an imperfect knowledge of h(φ) and ∆, it is to say of c h (φ), on the 1D-DSPE algorithm therefore denoted by:

θi = arg max θ 1 c H h (θ)V b 2 (4)

PERFORMANCE ANALYSIS

General case

The data covariance matrix Ry0 = E[y 0 y H 0 ] is estimated from N snapshots {y 0 (tn)}n=1,...,N by Ry0 ≈ 1 N N n=1 y 0 (tn)y H 0 (tn). It follows that the DOA estimator of θi is :

θi = arg max θ 1 c H h (θ) Vb0 2 , ( 5
)
where Vb0 is associated to the eigendecomposition of:

Ry0 = Vs0 Λs0 VH s0 + Vb0 Λb0 VH b0 , (6)
and c h (θ) is given by (2) but with a form of the angular spread distribution h(φ) which may be different from the true distribution of the actual sources h0i(φ). Note that the parameters related to the true values are indexed by 0 to distinguish them from parameters used by the estimator. The theoretical covariance matrix of the signal related to the actual steering matrix C h0 (θ) is thus defined as Ry0 = Vs0Λs0V H s0 + σ 2 b V b0 V H b0 , and we write Ry0 = R yh + ∆R yh , where R yh would correspond to the covariance matrix of the signal related to a steering vector c h (θ) and where

R yh = V sh Λ sh V H sh + σ 2 b V bh V H bh .
We introduce the error on the steering vector model ∆c = c h (θi) -c h0i (θi). We also note ∆V bh = V b0 -V bh where ∆V bh is an error on the noise eigenmatrix due to the model mismatch between c h and c h0 . Let us also introduce the error due to the finite number of snapshots N on the estimation of the data covariance matrix ∆Ry0 = Ry0 -Ry0 and the error on the noise corresponding eigenmatrix ∆V b0 = Vb0 -V b0 . The error due to the finite number of snapshots on the noise subspace projector is

∆Π b0 = Vb0 VH b0 -V b0 V H b0 , with Π b0 = V b0 V H b0 .
According to [START_REF] Weiss | Effects od modeling errors on the resolution threshold of the MUSIC algorithm[END_REF], the relation of ∆Π b0 and ∆R b0 is:

∆Π b0 = -Π b0 ∆Ry0Q -Q∆Ry0Π b0 , (7)
where

Q = Vs0(Λs0 -σ 2 b I) -1
Vs0, I is a q × q identity matrix. According to [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF], the DOA estimation θi satisfies that the first order derivative of the denominator of ( 5) equals zero which yields:

2Re{ċ H h ( θi) Vb0 VH b0 c h ( θi)} = 0, (8)
where ċh ( θi) = ∂c(θ) ∂θ | θi . Assuming that, θi is not far from θi, we introduce the second order Taylor series approximations of c h (θ) and ċh (θ) in ( 8), and exploiting [START_REF] Zoubir | Efficient subspace-based estimator for localization of multiple incoherently distributed sources[END_REF] to substitute Vb0 Vb0 yields:

A(θi)∆θ 2 i + B(θi)∆θi + C(θi) = 0, (9)
where the terms of order greater than 2 in ∆θi have been neglected, and the scalar A, B, C are defined with B = B1 + B2 and C = C1 + C2:

A(θi) =Re 1 2 ċh (θi) H V b0 V H b0 ch (θi) + ch (θi) H V b0 V H b0 ċh (θi) + 1 2 ... c h (θi) H V b0 V H b0 c h (θi) , B1(θi) =Re ċh (θi) H V b0 V H b0 ċh (θi) + ch (θi) H V b0 V H b0 c h (θi) , E[∆θ 2 i ] = B1(θi) 2 -2A(θi)C1(θi) -B1(θi) B1(θi) 2 -4A(θi)C1(θi) 2A(θi) 2 + σ 2 b N 1 4A(θi) 2 - B1(θi) 4A(θi) 2 B1(θi) 2 -4A(θi)C1(θi) ϕ(θi) + χ(θi) 2A(θi) B1(θi) 2 -4A(θi)C1(θi) , (11)
where:

ϕ(θi) =Re 4ċ(θi) H V b0 V H b0 ċ(θi) ċ(θi) H QRy0Qċ(θi) + 4c(θi) H V b0 V H b0 ċ(θi)c(θi)QRy0Qc(θi) + c(θi) H V b0 V H b0 c(θi)c(θi)QRy0Qc(θi) , χ(θi) =Re 2ċ(θi) H V b0 V H b0 ċ(θi)ċ(θi)QRy0Qc(θi) + c(θi) H V b0 V H b0 ċ(θi)c(θi)QRy0Qc(θi) . ( 12
) B2(θi) =Re 2ċ h (θi) H V b0 V H b0 ∆Ry0Qċ h (θi) +c h (θi) H V b0 V H b0 ∆Ry0Qc h (θi) , C1(θi) =Re ċh (θi) H V b0 V H b0 c h (θi) , C2(θi) =Re ċh (θi) H V b0 V H b0 ∆Ry0Qc h (θi) .
The expression of ∆θi can be obtained by solving the 2 nd order equation ( 9), and taking into account the convexity of the criterion to select the right root, so that:

∆θi = -B(θi) + B(θi) 2 -4A(θi)C(θi) 2A(θi) . (10)
In addition, according to the formulas in [START_REF] Marcos | Les méthodes à haute résolution, traitement d'antenne et analyse spectrale[END_REF], the MSE can be expressed by [START_REF] Swindlehurst | A performance analysis of subspace method in the presence of model errors : Case of MUSIC[END_REF], as shown at the bottom of the previous page, exploiting the fact that B2(θi) and C2(θi) are random terms related to ∆Ry0, while A(θi), B1(θi), C1(θi) are deterministic, and using E[∆Ry0] = 0.

It is interesting to see that the MSE is composed of two terms, one depends only on the model error, the other depends on the model error but with a factor σ 2 /N . It can be expected that the second term will be negligible when N increases.

First order approximation

In this subsection, we discuss the situation that the estimation error ∆θi is small enough, so that the second order terms in ∆θi can be negligible with respect to the first order terms. Keeping only the first order terms in [START_REF] Souden | A two-stage approach to estimate the angles of arrival and the angular spreads of locally scattered sources[END_REF], replacing V b0 by V bh + ∆V bh , neglecting the second order terms in ∆V b0 ∆θi and ∆V bh ∆θi, and using R yh V b0 = V b0 Λ bh , the estimation error can be simplified as:

∆θi = Re{ċ H h (θi)V bh V H bh ∆c} ċH h (θi)V bh V H bh ċh (θi) + Re{ċ H h (θi)V bh V H bh ∆Ry0Qc h (θi)} ċH h (θi)V bh V H bh ċh (θi) . ( 13
)
The MSE can be given by:

E[∆θ 2 i ] = Re{ċ H h (θi)V bh V H bh ∆c} ċH h (θi)V bh V H bh ċh (θi) 2 + σ 2 b 2N • Re{ċ H h (θi)V bh V H bh ċh (θi)c H h (θi)Ry0Qc h (θi)} (ċ H h (θi)V bh V H bh ċh (θi)) 2 . ( 14
)
As already noted for [START_REF] Swindlehurst | A performance analysis of subspace method in the presence of model errors : Case of MUSIC[END_REF], the second term decreases when N increases, which is not the case for the first terms of ( 13) and ( 14). Note that this simplification allows to express the first terms of ( 13) and ( 14) as a function of the model error ∆c.

Estimation error as an explicit function of the model error

In this subsection, we assume that the covariance matrix is perfectly known and we only study the deterministic terms of (10) so that B2 and C2 are null. We assume that in this case, the shape of h related to the actual signals sources is known, where all the sources have the same shape and the same angular spread dispersion ∆0. The model error is therefore caused by the error on ∆ or by the fact that the sources are assumed to be points when the classical MUSIC estimator is used. Assuming again that ∆ is not far from ∆0, and noting δ = ∆-∆0, we can introduce the second order Taylor series approximations in δ:

c h (θi) ≈ c h0 (θi) + δg h01 (θi) + 1 2 δ 2 g h02 (θi), (15)
where

g h01 (θi) = ∂c h0 (θ i ) ∂∆ |∆ 0 , g h02 (θi) = ∂ 2 c h0 (θ i ) ∂∆ 2
|∆ 0 . We can pay attention that g h01 and g h02 reveals the sensibility of our model to the variation of the angular dispersion of the actual signal. The estimation error can be thus given by:

∆θi = -(f4 + f5δ + f6δ 2 ) + x + yδ + zδ 2 + wδ 3 + kδ 4 2(f1 + f2δ + f3δ 2 ) , (16)
where fi(i = 1, . . . , 5) and x, . . . , k can be expressed in function of V b0 , derivatives of g and c. Details can be found in [START_REF] Xiong | Performance analysis of MUSIC for spatially distributed sources[END_REF].

Ignoring the second terms in ∆θi in ( 9) and extending Taylor series in δ to third order terms for an accurate approximation, the estimation error can be further expressed as a polynomial in δ:

∆θi = α(θi, ∆0)δ + β(θi, ∆0)δ 2 + γ(θi, ∆0)δ 3 , (17)
where α, β and γ are functions depending on the source signals and sensor parameters (see [START_REF] Xiong | Performance analysis of MUSIC for spatially distributed sources[END_REF] for details).

NUMERICAL SIMULATIONS

In this section, numerical examples are presented to illustrate the validity of the analytical results of the estimation performances established in section 3. In all simulations, a uniform linear array composed of M = 10 sensors spaced by d = λ/2 is considered, SN R = 10dB, and N = 1000 snapshots. Figure 1 exhibits the estimation error or bias on θ for two Uniform distributed sources at θ1 = 21 • and θ2 = 39 • , respectively, with an angular dispersion ∆0 = 10 • . The model error due to the dispersion parameter ∆ has been varied to study its effect on the DOA's estimation accuracy. As expected the bias is null when the exact model parameter ∆ = ∆0 = 10 • is used. We can observe that the estimation error is smaller when the parameter ∆ is smaller than ∆0, than when it is larger than ∆0. Focusing on the validity of the expression derived in the previous section, we can notice that for θ2, the DOA estimation error ∆θ2 obtained in [START_REF] Zoubir | Efficient DSPE algorithm for estimating the angular parameters of coherently distributed sources[END_REF] and the one obtained in [START_REF] Ferréol | Statistical analysis of the MUSIC algorithm in the presence of modeling errors, taking into account the resolution probability[END_REF] both match the simulation results. However, for θ1, which has an estimation error much bigger than θ2 when the model error increases, the result obtained in [START_REF] Zoubir | Efficient DSPE algorithm for estimating the angular parameters of coherently distributed sources[END_REF] with the second order approximation outperforms the one in [START_REF] Ferréol | Statistical analysis of the MUSIC algorithm in the presence of modeling errors, taking into account the resolution probability[END_REF] with the first order approximation.

In Figure 2, the RMSE of θ evaluated with 100 Monte Carlo simulations is plotted versus the snapshot number, with the model parameter set to ∆ = 15 • . The MSE decreases as well as N increases, and then converges to a non zero value whose expression is given in [START_REF] Swindlehurst | A performance analysis of subspace method in the presence of model errors : Case of MUSIC[END_REF]. This reveals that when there are two perturbations, the finite number of snapshots effect dominates in the case of a small number of snapshots whereas the model error effect dominates when the snapshots number is large. Again we can see that as in Figure 1, the result with the second order approximation outperforms the one with the first order approximation.

Figure 3 illustrates the approximation obtained in ( 16) and [START_REF] Xiong | Performance analysis of MUSIC for spatially distributed sources[END_REF], in order to express the estimation error as an explicit function of the model error, in the same previous scenario. We can observe that, in this case, [START_REF] Raich | Bearing estimation for a distributed source: Modeling, inherent accuracy limitations and algorithms[END_REF] fits better the simulation results, while Figure 4 shows the results of the absolute value of the DOA estimation error when the model error due to the dispersion parameter and the angle between the two sources both vary, with θm = 1 2 (θ1 + θ2) = 30 • . The green stars mark the region where the sources are not resolved, that is to say, the two sources are so close that the criterion 1D-DSPE gives the false appearance that there is only one source in the middle. These results allow to highlight two behaviors of the estimators: firstly, the closer the sources are, the more the model error impacts on the estimation accuracy; secondly, when the model error is small enough, high resolution is achieved. Taking into account that a distributed estimator with ∆ = 0 • means the classical estimator MUSIC with a point source model, and ∆ = 10 • is the 1D-DSPE without model error, the advantage of a distributed estimator with respect to a classical MU-SIC is highlighted. In addition to the fact that the theoretical results correspond to the simulation results, we can observe that again the result obtained in [START_REF] Zoubir | Efficient DSPE algorithm for estimating the angular parameters of coherently distributed sources[END_REF] works better than the result obtained in [START_REF] Ferréol | Statistical analysis of the MUSIC algorithm in the presence of modeling errors, taking into account the resolution probability[END_REF], even in the case where there are resolution problems.

CONCLUSION

In this paper, we have investigated the effects of both the angular dispersion of the source and the finite number of snapshots on the behavior of the DOA MUSIC-based estimator. New analytical expressions of the DOA estimation error and MSE as a function of these two perturbations have been given. Particularly, in the special case when the theoretical covariance matrix is available, expressions as an explicit function of the model error is proposed, which gives an easier way to analyze the influence of a model error, or to optimize the array configurations to reduce the DOA estimation error. Simulations which are carried out are in adequacy with the proposed theoretical results. The performance of MUSIC for coherently distributed sources can thus be analyzed.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: DOA estimation error vs. angular dispersion model error (2 sources with Uniform angular dispersion, ∆0 = 10 • , θ1 = 21 • , θ2 = 39 • )

Figure 4 :

 4 Figure 4: Absolute value of estimation error |∆θ1| vs. source angular separation |θ2 -θ1| and model angular dispersion ∆ (2 sources with Uniform angular dispersion, ∆0 = 10 • , 100 Monte-Carlo simulations)


Keywords: Smart Camera, Multi-core processor, Real-time Image processing, Development Framework, Shared Memory Management, Building Information Modeling (BIM)

A New Development Framework for Multi-Core Processor based Smart-Camera Implementations Serge Bobbia, Julien Dubois, Christophe Nicolle, Axel Moinet, Benoit Darties Laboratory Le2i, UMR CNRS 6306 University of Bourgogne, France Contact Email: julien.dubois@u-bourgogne.fr, http://le2i.cnrs.fr/-Membres-The exponential evolution of the smart camera processing performances is directly linked to the improvements on hardware processing elements. Nowadays, high processing performances can be reached considering hardware targets which enables a high level of task parallelism to be implemented. Highly regular tasks are good candidate for a reconfigurable logic implementation and less regular parts of the algorithm could be described on the processor. Meanwhile the prototyping time is related to the selected target and the associated development methodology. The implementation on reconfigurable logic is highly efficient in exploiting the intrinsic task parallelism nevertheless can be time consuming using traditional methodology (i.e. Hardware Language Description). Several approaches can be considered to decrease the prototyping time and to conserve high processing performances for instance implementation based on:

• heterogeneous architectures [START_REF] Mosqueron | Smart camera based on embedded HW/SW coprocessor[END_REF] that mixed reconfigurable logic (i.e. FPGA) and embedded processor,

• high-level abstraction description and the associated fast prototyping tools [START_REF] Thavot | Dataflow design of a co-processor architecture for image processing[END_REF][3][4],

• multi-core processor architectures such as Digital Signal Processors (DSP), Graphic Processor Units (GPU) or even Generic Purpose Processor (GPP).

In this paper, we propose to focus on implementation based on GPP due to the emergence of new generation of low-cost multi-core processors which enables high processing performances to be reached and therefore to match with some constraints of complex image-processing algorithms. The key idea of this development is to be able to propose fast prototyping using a low-cost smart camera based on this kind of target. Hence, we have developed a new framework dedicated to multi-core processor associated with an image sensor. The framework aims to offer a high degree of flexibility for managing the tasks and the memory allocation. Hence, the framework enables the priority and the allocation of each task to be controlled. Each task (or binary) is independent in terms of execution nevertheless it can be linked and controlled using a higher hierarchy level binary. The image acquisition task can be completely independent from the other processing tasks. One processor's core can even be dedicated to the acquisition task to guarantee a constant input data-flow to the image processing tasks. The data exchange is defined in POSIX, each binary can be therefore coded differently (for instance in C or C++, or in another languages) and offer a relative Operating System (OS) compatibility. The memory management enables a sequence of images to be automatically stored and a simultaneous access to be granted for several processings. The framework includes an interface dedicated to the management of the tasks: the user can add or suppress a binary during the runtime, logs or processing results can be visualised for each task.

The resulting framework has been used with the low-cost multi-core processor Raspberry2 platform coupled with its specific camera (1080p@30fps). The proposed architecture offers a Quad-Core ARM Cortex-A7 900 MHz associated to a Dual-Core VideoCoreIV GPU which enables to perform 1080p video decoding. The 1 Go memory space is double that of the previous Raspberry product. Using Operating Systems (OS), the resulting platform enables the user to access a low-cost embedded computer. The low-cost smart camera associated to the proposed framework enables fast prototyping of complex image processing. We focus on the extraction and the exploitation of visual information for smart building management. For instance, the people tracking would be a very useful processing in order to increase the level of security or services in the building. The Building Information Modeling (BIM) [START_REF] Volk | Building Information Modeling (BIM) for existing buildings Literature review and future needs[END_REF] regroups all the data which are involved in the building life-cycle management [START_REF] Vanlande | IFC and building lifecycle management[END_REF]. This priori information could be used, for instance, to simplify the tracking process in a multi-view configuration.