Louis Blanchard

Régis Duvigneau

Anh-Vu Vuong

Bernd Simeon

Shape Gradient Computation in Isogeometric Analysis for Linear Elasticity

Keywords: Shape derivative, Isogeometric analysis, Linear elasticity, Structural shape optimization * Opale Project-Team

The transfer of geometrical data from CAD (Computer Aided Design) to FEA (Finite-Element Analysis) is a bottleneck of automated design optimization procedures, yielding a loss of accuracy and cumbersome software couplings. Isogeometric analysis methods propose a new paradigm that allows to overcome these difficulties by using a unique geometrical representation that yields a direct relationship between geometry and analysis. In this study, we investigate its use for sensitivity analysis and more specifically shape gradient computations, in the framework of linear elasticity problems. The potential of isogeometric analysis methods for shape gradient computations is demonstrated for two-and three-dimensional design problems.

Calcul de gradient de forme en analyse isogéométrique pour l'élasticité linéaire

Résumé :

Le transfert de données géométriques depuis la CAO (Conception Assistée par Ordinateur) vers la simulation par éléments-finis est un goulot d'étranglement des procédures automatiques de conception optimale, conduisant à une perte de précision et à des couplages de codes pénibles. Les méthodes d'analyse isogéométrique proposent un nouveau paradigme, qui permet de surmonter ces difficultés en utilisant une unique représentation géométrique qui conduit à une relation directe entre la géométrie et la simulation. Dans cette étude, on explore son utilisation pour l'analyse de sensibilité et plus spécifiquement le calcul de gradient de forme, dans le cadre de problèmes en élasticité linéaire. Le potentiel des méthodes d'analyse isogéométrique pour le calcul du gradient de forme est démontrée pour des problèmes de conception bi-et tridimensionnels.

Mots-clés : Gradient de forme ; analyse isogéométrique ; Elasticité linaire ; optimisation de forme structurale

Introduction

This study concerns shape optimization problems for systems governed by PDEs (Partial Differential Equations) as encountered in engineering applications, e.g., in structural mechanics, fluid dynamics or electromagnetics. Multi-disciplinary optimization procedures are nowadays commonly used by engineers to solve complex problems in aeronautical or automotive industries, for instance. However, the processing of the geometrical data through the design loop is still a delicate issue: several different representations of the geometry coexist and several conversions are required, yielding additional overhead and extra difficulties from the theoretical and practical point of view.

Usually, the geometry of the engineering system of interest is defined using CAD (Computer Aided Design) software, on the basis of high-order representations. NURBS (Non-Uniform Rational B-Spline) basis functions are considered as standard in this context [START_REF] Sabin | Cad system components[END_REF]. Cubic functions are often used, exhibiting a C 2 regularity.

The predominant analysis tools used to solve PDEs, such as FEA (Finite-Element Analysis) software, rely on a grid to describe the computational domain. The construction of this grid accounts for the geometrical data from CAD. However, a different representation is employed, since the grid usually defines the geometry with piecewise linear functions and C 0 regularity only. The task of generating the grid from CAD data is one of the most time-consuming in industry. Moreover, it is remarkable that some errors are introduced before any physical analysis has started.

Once the grid is available, the PDEs can be solved, and the resulting solution fields allow to estimate the performance of the current design. The optimizer can finally propose an update of the geometry in the CAD system to improve the performance. Here again, some inconsistencies occur, since the performance estimation is based on the approximate piecewise linear geometry, whereas the optimizer updates the high-order NURBS representation. The same inconsistencies hold for gradient estimates, if a gradient-based optimization method is employed.

The coexistence of these two geometrical representations makes the practical implementation of the design loop unnecessarily complex, since additional methods and software must be introduced: the automated update of the grid for any geometry change, the exchange of data between CAD and mesh generation tools as soon as automatic refinement is done on the boundary, the exchange of data between the solvers (multidisciplinary problems) as soon as the grid at the interface is modified, re-computation of the grid interface between different domains for problems with rotating elements, etc.

In the context of gradient-based optimization, the use of two geometrical representations makes also the estimation of the performance derivatives more complex, since the derivatives of the grid motion with respect to the NURBS parameters should be computed in the whole computational domain. This task may require to build the adjoint system of the mesh deformation method [START_REF] Nielsen | Using an adjoint approach to eliminate mesh sensitivities in computational design[END_REF].

The isogeometric analysis approach, proposed recently by Hughes et. al. [START_REF] Cottrell | Isogeometric analysis : towards integration of CAD and FEA[END_REF][START_REF] Cottrell | Studies of refinement and continuity in isogeometric analysis[END_REF][START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry, and mesh refinement[END_REF], promises to resolve or at least alleviate all these issues. It consists basically in solving PDEs with a variational approach on the basis of a NURBS representation originating directly from CAD. Thus, this approach promotes the integration of the two so far disjunct disciplines of FEA and CAD. Although this approach is especially beneficial for automated shape optimization procedures, as explained above, only a few studies can be found in the literature concerning its use for design problems [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF][START_REF] Seo | Isogeometric topology optimization using trimmed spline surfaces[END_REF][START_REF] Wall | Isogeometric structural shape optimization[END_REF]. Therefore, this study is devoted to the estimation of shape derivatives in this particular context. The concept of shape derivatives originates from Hadamard [START_REF] Hadamard | Mémoire sur le problème d'analyse relatif à l'équilibre des plaques élastiques[END_REF]. It has been studied from a theoretical point of view by several authors [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF], and also used to solve practical shape optimization problems (see for instance [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF]).

In this study, we consider the framework of the linear elasticity equations, as outlined in Section 2. The application of isogeometric analysis to linear elasticity problems is described in Section 3 and validated on a typical problem in Section 4. Then, we introduce the concept of shape derivatives and its use in the isogeometric analysis context in Section 5. The resulting shape optimization procedure is described in Section 6 and validated in Section 7, for a problem with an analytical solution. Finally, a three-dimensional application is shown in Section 8 to demonstrate the capability of the proposed approach to tackle realistic problems.

Linear elasticity problem

We consider an open domain Ω in d = 2 or d = 3 dimensional space, which represents a deformable solid subject to external forces. Its boundary is composed of three disjoint parts, ∂Ω = Γ N ∪Γ D ∪Γ with Γ∩Γ D = ∅, Γ∩Γ N = ∅ and Γ N ∩Γ D = ∅. Dirichlet (imposed displacement) and Neumann (imposed stress) boundary conditions are prescribed on Γ D and Γ N , respectively, whereas Γ is considered as optimization variable (moving boundary). A zero Neumann boundary condition is prescribed on Γ.

The governing equations are the linear elasticity equations under the small deformation assumption. In case of zero body loads, the displacement field u is therefore the solution of the following system

       -div σ(u) = 0 in Ω u = 0 on Γ D σ(u) • n = g on Γ N σ(u) • n = 0 on Γ (1)
where n is the outward unit normal vector and σ(u) the second-order stress tensor defined by Hooke's law σ(u) = 2 µ (u) + λ T r (u) Id.

Moreover, λ and µ are the Lamé parameters of the material and

(u) = 1 2 ∇u + ∇u T
is the linearized strain tensor. Like in the finite element method, the variational formulation of the linear elasticity problem defined by Eq. ( 1) is the starting point for isogeometric analysis. The weak form of (1) reads as follows: we seek for the physical displacement field u ∈

V = ϕ ∈ H 1 (Ω) d , ϕ = 0 on Γ D such that Ω 2 µ (u) : (v)+λ divu divv dΩ = Γ N g . v dΓ ∀v ∈ V. (3)
By projecting the solution u and the test functions v to a finite-dimensional subspace V h of V we next construct a numerical approximation.

Isogeometric analysis

The governing equations described in the previous section can be easily solved by means of a classical FEA approach. A suitable discretization of the domain is usually constructed by Inria introducing a mesh Ω h that approximates the geometrical domain Ω. As demonstrated by Hughes et al. [START_REF] Cottrell | Isogeometric analysis : towards integration of CAD and FEA[END_REF][START_REF] Cottrell | Studies of refinement and continuity in isogeometric analysis[END_REF][START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry, and mesh refinement[END_REF], this geometrical approximation of the domain has several drawbacks. In a few words, the use of an approximate geometry introduces some errors that could be avoided, the construction of the mesh is time consuming, and the coupling between different scientific disciplines is made unnecessarily complex.

Isogeometric analysis [START_REF] Cottrell | Isogeometric analysis : towards integration of CAD and FEA[END_REF][START_REF] Cottrell | Studies of refinement and continuity in isogeometric analysis[END_REF][START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry, and mesh refinement[END_REF] is aimed at a better integration of FEA and CAD methods by using a unique representation basis for the geometry and the discrete solution fields. Isogeometric analysis proposes to discretize the computational domain exactly by using a NURBS basis originating from CAD.

NURBS representation

NURBS basis functions are defined in a parametric domain Ω 0 . They can be represented in the physical domain Ω by introducing the transformation

F : Ω 0 → Ω, F (ξ) = x(ξ). (4)
Thus, any point of coordinates x = (x, y) T in the physical domain Ω is mapped to a point of parameters ξ = (ξ, η) T in the parametric domain Ω 0 , as illustrated by Fig. 1. For the sake of simplicity, only the two-dimensional case is described here, but the extension to the threedimensional case is straightforward. Two-or three-dimensional NURBS basis functions are defined as the bivariate or trivariate tensor product of one-dimensional basis functions. The latter have compact support and are constructed using the so-called knot vector Ξ = (ξ 0 , . . . , ξ l ) ∈ R l with l = n + p + 1, which consists of nondecreasing real numbers. Here, p is the degree of the functions and n the number of functions considered. Open knot vectors, i.e., knot vectors with first and last knots of multiplicity p + 1, are typically used [START_REF] Boor | A Practical Guide to Splines[END_REF][START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design[END_REF]. NURBS basis functions are rational extensions of B-Spline basis functions N i,p i = 0, • • • , n that are computed recursively as

N i,0 (ξ) = 1 if ξ i ≤ ξ < ξ i+1 0 otherwise (5) N i,p (ξ) = ξ -ξ i ξ i+p -ξ i N i,p-1 (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 N i+1,p-1 (ξ). (6)
Note that the quotient 0/0 is assumed to be zero.

One-dimensional NURBS basis functions of degree p are then given by

R i,p (ξ) = w i N i,p (ξ) j∈J w j N j,p (ξ) (7)
where w i ∈ R is the weight associated to the ith function, and J = {0, . . . , n}. Two-dimensional NURBS basis functions are constructed using a bivariate tensor product (we omit the degrees p ξ and p η for the sake of readability)

R k (ξ, η) = w k N k (ξ)N (η) i∈I j∈J w ij N i (ξ)N j (η) . ( 8
)
This representation requires two knot vectors Ξ ξ = (ξ 0 , . . . , ξ l1 ) and Ξ η = (η 0 , . . . , η l2 ). The extension to three-dimensional basis functions is straightforward. More details on the properties of NURBS can be found in [START_REF] Boor | A Practical Guide to Splines[END_REF][START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design[END_REF]. The transformation of the parametric domain Ω 0 to the physical domain Ω in Eq. ( 4) is defined by associating a control point to each basis function

x(ξ, η) = i∈I j∈J R ij (ξ, η)X ij , (9)
where X ij ∈ R 2 represents the coordinates of the control point of indices (i, j) in the physical domain. Using such a representation, the domain Ω is described as a single NURBS patch. The shape of the domain Ω is thus defined by the position of these control points, the weights, the knot vectors and the degrees of the basis. NURBS representations are nowadays considered as standard for geometric modeling in CAD tools [START_REF] Sabin | Cad system components[END_REF]. Thus, we assume this representation as exact. Isogeometric analysis proposes to employ such a parametric representation to solve the governing equations, without approximating the domain by a piecewise linear grid.

Galerkin projection

In the isogeometric paradigm, the same representation is employed for both the geometry and the physical fields. Therefore, according to Eq. ( 9), the discretized displacement field u h is constructed as linear combination of the NURBS that define the geometry, i.e.,

u h (x) = i∈I j∈J R ij (ξ)U ij = i∈I j∈J Rij (x)U ij (10)
where we define the basis functions Rij in the physical domain by

Rij (x) = Rij (x, y) = Rij • F (ξ, η) = R ij (ξ). (11)
Note that the unknowns U ij in this representation are two-dimensional vectors and comparable to control points. In general, they do not stand for displacements in specific nodes, as in the finite element method. Due to the usage of open knot vectors, however, the bivariate NURBS are interpolatory on the boundary, and thus the zero boundary conditions u = 0 on Γ D are easily enforced by setting those coefficients U ij to zero that belong to the corresponding knots on the boundary Γ D .

Inria

Eliminating those basis functions that are required to enforce the zero boundary conditions and renumbering the remaining basis functions and unknowns from = 1 to m, we can write the numerical approximation in a standard fashion as

u h (x) = m =1 N (x)U = N(x)U (12)
where U ∈ R m contains the unknown displacement coefficients and each column N (x) : R 2 → R 2 of the 2 × m matrix N(x) the NURBS basis functions R ij , split into the components for the x and y-directions. In this way, the finite-dimensional subspace for the Galerkin projection is given by

V h = span {N , = 1, . . . , m}. (13)
Consequently, the restriction of the variational formulation (3) to V h leads to a system of linear equations

K U = F, K ij = Ω 2 µ (N i ). (N j ) + λ div(N i ) div(N j ) dΩ, F = Γ N g . N dΓ. (14)
Compared to standard finite elements, the m × m stiffness matrix K is less sparse in general, which is due to the larger support of the NURBS. The entries of the stiffness matrix and also those of the force vector F ∈ R m are computed by applying classical quadrature rules. Actually, this integration is performed in the parametric domain, in accordance with the evaluation process of the NURBS in (5-7). Since the basis functions are defined on each knot interval, these knot intervals are often considered as elements in the isogeometric analysis terminology. Therefore, the above integration is carried out as usually in an element-wise way.

For the solution of the linear system with symmetric positive definite stiffness matrix, various standard methods can be applied. As long as the dimension remains moderate, a sparse direct solver, such as the multi-frontal method implemented in the UMFPACK library [START_REF] Davis | A combined unifrontal/multifrontal method for unsymmetric sparse matrices[END_REF], is a good choice. The reader will find more details on the implementation in [START_REF] Cottrell | Isogeometric analysis : towards integration of CAD and FEA[END_REF]. Techniques for local refinement are discussed, e.g., in [START_REF] Dorfel | Adaptive isogeometric analysis by local h-refinement with t-splines[END_REF][START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF].

Remarks on the convergence properties of isogeometric analysis

The variational framework of the finite element method applies also to isogeometric analysis. In particular, this means that the best approximation property

u -u h E ≤ u -v E ∀v ∈ V h
holds with respect to the energy norm

u 2 E := Ω 2 µ (u). (u) + λ div(u) div(u) dΩ,
which is equivalent to the semi-norm |u| H 1 (Ω) in H 1 (Ω). Moreover, the numerical approximation u h approaches the energy norm of the exact solution from below,

u h E ≤ u E .

RR n°8111

Note that the energy norm of the numerical solution is available for free since it can be computed from the quadratic form u h 2 E = U T KU. When refining a discretization, it must hold

u h1 E ≤ u h2 E
where h 1 stands for the coarse and h 2 for the refined mesh. This property is often used to cross-check the numerical solution from a practical point of view.

The theoretical convergence properties of isogeometric analysis are quite similar to those established for the finite element method. We outline here the results of [START_REF] Bazilevs | Isogeometric analysis: Approximation, stability and error estimates for h-refined meshes[END_REF], assuming an invertible geometry parametrization F that is at least of class C 0 . Let Q denote the analogue of an element in parametric space, i.e., Q is the cartesian product of d = 2 (or d = 3) non-empty knot spans. The function F maps this rectangle to K = F (Q), which defines an element in physical space, with possibly distorted shape. If the degree of the NURBS used for the geometry and the discretization is p, then it holds at the element level

|u -u h | H 1 (K) ≤ Ch p p+1 i=0 ∇F i-p-1 L ∞ ( Q) |u| H i ( K) . (15)
Here, C is a constant that depends on the geometry of Ω but not on the meshsize parameter h, which is directly related to the size of Q. Furthermore, Q and K, respectively, denote the support extension of Q and K, which is the union of the supports of those basis functions that are non-zero over Q.

By summing up over all elements, convergence of order p can then be established with respect to the H 1 -seminorm. We remark that this estimate requires bounds for the i-th order seminorm of the solution u, unlike for finite elements. Also, the parametrization enters the error bound and thus has a clear effect on the quality of the solution. If we measure the error in the H k -seminorm where k > 1 and if the smoothness of the parametrization is less than C k-1 , the global error estimate requires the concept of so-called bent Sobolev spaces that take the reduced smoothness of the parametrization along specific knot lines into account, see [START_REF] Bazilevs | Isogeometric analysis: Approximation, stability and error estimates for h-refined meshes[END_REF] for more details.

As an example we look at the infinite plate with a hole already introduced in [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry, and mesh refinement[END_REF]. The shape and the boundary conditions are illustrated in Fig. 2. Symmetry along the coordinate axes is used, and the exact traction is employed for the Neumann boundary conditions. The underlying geometry parameterization is C 1 . The convergence plot for the first principle stress σ x for different degrees is shown in Fig. 3 and turns out to be as predicted by Eq. ( 15): the convergence rate increases with the degree.

Shape gradient computation

Optimization problem

We consider as optimization problem the minimization of the compliance subject to a constant mass constraint. The cost function, denoted J(Ω), depends on the shape of the domain Ω and is defined as

J(Ω) = Γ N g . u dΓ. (16)
Therefore, the constrained shape optimization problem reads where U ad is the set of admissible shapes and V 0 is the volume of the initial domain, considered as a reference volume. We refer to [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF] to define precisely the conditions of existence and uniqueness of the solution for such a problem. To account for the constraint, we transform the problem defined in Eq. ( 17) into an unconstrained optimization problem using a Lagrangian formulation

min Ω∈U ad J(Ω) , U ad = Ω ⊂ R d , Ω dΩ = V 0 (17)
min Ω∈R n L(Ω, β) , L(Ω, β) = J(Ω) + β Ω dΩ -V 0 ( 18
)
where β > 0 is a positive Lagrange multiplier.

Shape derivative

In order to solve the problem in Eq. ( 18) using a gradient-based method, we introduce the concept of shape derivative. The theoretical foundations of this concept can be found in [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF]. Let v : R d → R d denote a sufficiently smooth admissible vector field that defines the shape deformation due to the optimization process. Let T t : R d → R d denote the transformation which maps the initial domain Ω into the moving domain Ω t = T t (v)(Ω) := Ω + tv(Ω). Note that this deformation is different from the mechanical deformation of the system. The mapping introduced above is defined in physical space and is a priori independent of the mapping used in isogeometric analysis between the parametric domain Ω 0 and the initial physical domain Ω. However, as we will see later, the same NURBS basis functions will be used in practice, according to the isogeometric paradigm.

We describe next the main ingredients that allow to establish the expression of the shape derivative. Further details can be found in, e.g., [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF]. Since the minimization of the compliance for a structural system governed by linear elasticity equations is self-adjoint, the computation of the shape derivative does not require the solution of an adjoint problem. By differentiating Eq. ( 16), under the assumptions that g, u ∈ H 2 (Ω), it can be shown that the shape derivative

dJ(Ω; V ) = - Γ 2 µ | (u)| 2 + λ |divu| 2 v.n dΓ. (19)
One can notice that the compliance can always be decreased by enlarging the domain: v.n > 0 implies that dJ(Ω; v) ≤ 0. Similarly, the shape derivative of the Lagrangian function reads

dL(Ω, β; V ) = Γ β -2 µ | (u)| 2 + λ |divu| 2 v.n dΓ. ( 20
)
Therefore, a possible choice to ensure a decreasing compliance is to consider a shape deformation v defined by

v = 2 µ | (u)| 2 + λ |divu| 2 -β n on Γ 0 on Γ D ∪ Γ N ( 21
)
where β is chosen to satisfy the constraint.

Eq. ( 21) defines the shape deformation of the boundaries ∂Ω = Γ∪Γ D ∪Γ N . However, it would be convenient to extend this definition on the whole domain Ω, in order to define the domain transformation T t . From a practical point of view, this can be used to update the geometry as a whole. Consequently, we choose to extend the deformation field v using a linear elasticity model, which enables to re-use the structural elasticity solver. Finally, the deformation field is defined by solving

   -div(σ(v)) = 0 in Ω v = 0 on Γ D ∪ Γ N σ(v) • n = (2µ| (u)| 2 + λ|divu| 2 -β)n on Γ (22)
One can notice that, in Eq. ( 22), we do not impose the deformation field v on the boundary Γ as a Dirichlet condition, but instead we consider a Neumann condition for v. Actually, this choice is justified by two reasons. First, from numerical point of view, the deformation v will be defined using NURBS basis functions. Usually, the deformation defined by Eq. ( 21) does not belong to the NURBS space. Second, it is well known that the algorithms based on shape derivatives may suffer from a loss of regularity [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF], that should be corrected by applying a smoothing operator. The use of such a Neumann boundary condition enables to increase the regularity of the deformation field v.

Inria

Concluding, the problem defined by Eq. ( 22) is solved by using exactly the same approach as for the linear elasticity problem given by Eq. ( 1). In particular, an isogeometric approximation is employed, yielding a deformation of the domain in terms of control point displacements,

v(x) = m =1 N i (x)V i . (23)
The components V i that constitute the motion of the control points are thus obtained by solving a linear system similar to Eq. ( 14),

K V = C 1 -β C 2 , C 1 i = Γ 2 µ | (u)| 2 + λ |divu| 2 N i . n dΓ, C 2 i = Γ N i . n dΓ. (24)
The stiffness matrix K is defined by equation ( 14), V is the vector of the deformation control points and C 1 , C 2 are the external force vectors, which include the terms for the shape derivative and constant volume constraint. Actually, the system above is solved twice to compute K -1 C 1 and K -1 C 2 . Then, for any Lagrange multiplier β, the deformation field V can be computed without extra cost by

V = K -1 C 1 -β K -1 C 2 . ( 25
)
6 Shape optimization algorithm

The shape optimization procedure used in the numerical experiments presented in the next sections can be summarized by the following algorithm:

1. Initialization of the domain defined by the control points X 0 and counter k = 0 2. Start of optimization loop 3. Solve system [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF] for the domain defined by the control points X k , yielding the structural displacement defined by the control points U k 4. Solve system (24) for the domain defined by the control points X k , yielding the shape deformation defined by the control points V k 5. Initialization of step length t k 6. Start of the line search 7. Compute the Lagrange multiplier β k such that the domain defined by the control points [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF] for the domain defined by the control points X k+1 , yielding the structural displacement defined by the control points U k+1

X k+1 = X k + t k V k is admissible 8. Solve system
9. Compute the compliance for the domain defined by the control points X k+1

10. If the compliance decreases then the step length t k is retained, else t k is reduced and goto

11. End of line search 12. If converged then stop, else k ← k + 1 and goto (3)

End of optimization loop

We underline the fact that, according to the isogeometric analysis paradigm, a unique basis is used to represent the geometry of the domain, the structural displacement field and the deformation field for the optimization. Nevertheless, the hierarchical property of NURBS can be used to refine this representation if required. Typically, an accurate resolution of the elasticity equations requires the use of a refined representation. On the other hand, the shape optimization problem can be solved by using a few degrees of freedom.

Inria 7 Shape optimization validation

A two-dimensional test-case is first considered to validate the shape gradient computation and the optimization procedure. It deals with a classical shape optimization problem in linear elasticity, for which the solution is known analytically. We consider a square flat plate with a hole located at its center, subject to uniform external normal forces. The objective is to determine the complete shape that minimizes the compliance for a constant plate area. The initial shape of the moving boundary is defined by four straight lines. The configuration of this test-case is illustrated by Fig. 5. Only a quarter of the geometry is considered for symmetry reasons. Additional geometrical constraints are thus introduced to maintain the extremities of the moving boundary on the symmetry axes.

We define the computational domain by a single bi-quadratic patch, which exhibits a singular point at the top left corner, as shown in Fig. 6. The 12 × 7 net of control points can be seen in this figure. Since quadratic basis functions are used, two superimposed control points permit to create the singularity at the top left corner. The optimization problem counts 24 variables that correspond to the coordinates of the 12 control points defining the moving boundary (in red in the figure). This test-case has already been studied using an isogeometric approach in [START_REF] Wall | Isogeometric structural shape optimization[END_REF], but without using the shape gradient concept. The shape gradient defined on the moving boundary is depicted in Fig. 7 for the initial domain and in Fig. 8 after 11 iterations. One can notice in particular the balance between the motion to reduce the compliance and the motion to maintain a constant area. As seen, the moving boundary shape evolves progressively towards a circle, which is the solution of the problem. The final domain is depicted in Fig. 9. As can be observed, the domain parameterization remains perfectly symmetric. One can notice a crossing of some control points in the vicinity of the hole. However, the parameterization remains injective which enables to solve the systems.

The fast convergence of the optimization procedure is illustrated in Fig. 10. Convergence is achieved in about 20 iterations, with results that are very close to those presented in [START_REF] Wall | Isogeometric structural shape optimization[END_REF]. This shape optimization exercise allows to assess the shape gradient computation which steers the optimization procedure in the isogeometric context.

Inria 8 Application to 3D problems

The proposed methodology is finally applied to a three-dimensional design problem. The objective is to optimize the shape of an open-spanner, which is subject to Dirichlet conditions at one extremity (zero displacement) and Neumann conditions at the other extremity (imposed external forces), as illustrated by Fig. 11. The initial domain is depicted in Fig. 12. A similar two-dimensional problem has been proposed in [START_REF] Wall | Isogeometric structural shape optimization[END_REF]. For this three-dimensional extension, two problems are considered successively: in a first optimization exercise, only the top and bottom boundaries are considered as optimization variables, whereas in a second exercise the shapes of the lateral boundaries are also optimized. The objective is again the reduction of the compliance, subject to a constant volume constraint. Additional constraints are also imposed to maintain the shape of the open-spanner symmetric with respect to horizontal and vertical planes. Since the boundary conditions are non-symmetric, the computational domain represents the whole open-spanner, as shown in Fig. 12. The domain is parameterized using linear basis functions in the vertical and crosswise directions and quadratic functions in the lengthwise direction. The net of control points counts 6 × 7 × 2 points in horizontal, vertical and crosswise directions, respectively, and is depicted in Fig. 13 for the initial domain. Thus, the two optimization exercises account for 48 variables (16 control points) and 192 variables (64 control points), respectively.

The evolution of the cost function during the optimization procedure is displayed in Fig. 14. Convergence is achieved in a few iterations and for a very low computational time. Obviously, a lower compliance value is reached for the second optimization exercise, which includes a larger number of variables. A comparison of initial and final domains is presented in Fig. 15 for the first exercise and in Fig. 17 for the second one. The displacement fields are depicted in Fig. 16 and Fig. 18. The optimization procedure yields an enlargement of the left part of the openspanner, whereas the right part becomes thinner to fulfill the constant volume constraint. For the first exercise, one can notice that the lateral boundaries are slightly modified, although they are not considered as optimization variables. This is due to the fact that the side boundaries are considered as Neumann conditions for the Eq. ( 22) that defines the shape evolution. The shape obtained is very close to that found in [START_REF] Wall | Isogeometric structural shape optimization[END_REF] for a two-dimensional problem. For the second exercise, the open-spanner is made thinner in the crosswise direction, which enables to enlarge it significantly. A closer look at the displacement fields shows that the maximum displacement value, at the right extremity, is significantly reduced during the optimization, especially in the second exercise.

Conclusion

This study has demonstrated that isogeometric analysis provides an attractive theoretical and practical framework for sensitivity analysis and shape derivative computation. The concept of shape derivatives has been set-up in this context for linear elasticity problems and used for practical shape optimization exercises.

The main ingredient is the use of the same NURBS representation for the computational domain, the structural displacement field, the shape gradient on the moving boundary and the domain deformation during the design process. The uniqueness of the geometrical representation makes the implementation straightforward and avoids the introduction of spurious errors in the design optimization loop.

Nevertheless, several improvements are still required to face industrial problems: from a practical point of view, the approach needs to be extended to multi-patch representations, to describe more complex geometries, and local refinement techniques should be introduced to improve computational efficiency. From a methodological point of view, more complex problems should be considered in the future, for instance exhibiting non-linear behavior.

Figure 1 :

 1 Figure 1: Mapping of the parametric domain Ω 0 to the physical domain Ω = F (Ω 0 ).

Figure 2 :

 2 Figure 2: Description of plate example

Figure 3 :

 3 Figure 3: Convergence plot for the plate with a hole for different ansatz degrees

Figure 4 :

 4 Figure 4: Mapping of the initial domain Ω to the moving domain Ω t = T t (v)(Ω).

Figure 5 :

 5 Figure 5: Description of the validation test-case.

Figure 6 :

 6 Figure 6: Initial domain.

Figure 7 :

 7 Figure 7: Shape gradient for initial domain.

Inria

Figure 8 :

 8 Figure 8: Shape gradient at the11th iteration.

Figure 9 :

 9 Figure 9: Final domain.

Figure 10 :

 10 Figure 10: Evolution of cost function.

Figure 11 :

 11 Figure 11: Description of the open-spanner test-case.

Figure 12 :

 12 Figure 12: Initial domain.

Figure 13 :

 13 Figure 13: Initial domain parameterization.

Figure 14 :

 14 Figure 14: Cost function evolution for the two optimization exercises with respect to CPU time.

Figure 15 :

 15 Figure 15: Comparison of the initial (white) and final (blue) domains for the first exercise.

Figure 16 :

 16 Figure 16: Comparison of the initial (bottom) and final (top) displacement fields for the first exercise.

Figure 17 :

 17 Figure 17: Comparison of the initial (white) and final (blue) domains for the second exercise.

Figure 18 :

 18 Figure 18: Comparison of the initial (bottom) and final (top) displacement fields for the second exercise.

3 . 2 6 4 7 5gradient computation 8 5 . 1 9 6

 3267519 representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Galerkin projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Remarks on the convergence properties of isogeometric analysis Shape Optimization problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 5.2 Shape derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Shape 2004 route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

RR n°8111

Acknowledgements

This study was supported by the European Union within the Projects 218536 "EXCITING" and 284981 "TERRIFIC" ( 7th Framework Program).

Inria


Stéphanie Lefèvre

Christian Laugier

Javier Ibañez-Guzmán

Risk Assessment at Road Intersections: Comparing Intention and Expectation

Intersections are the most complex and hazardous areas of the road network, and 89% of accidents at intersection are caused by driver error. We focus on these accidents and propose a novel approach to risk assessment: in this work dangerous situations are identified by detecting conflicts between intention and expectation, i.e. between what drivers intend to do and what is expected of them. Our approach is formulated as a Bayesian inference problem where intention and expectation are estimated jointly for the vehicles converging to the same intersection. This work exploits the sharing of information between vehicles using V2V wireless communication links. The proposed solution was validated by field experiments using passenger vehicles. Results show the importance of taking into account interactions between vehicles when modeling intersection situations.

I. INTRODUCTION

Intersection-related accidents account for 40-50% of road crashes in most countries, hence the interest shown by numerous safety improvement programs in the past ten years [START_REF] Fuerstenberg | Advanced intersection safetythe EC project INTERSAFE[END_REF], [START_REF] Maile | Objective testing of a cooperative intersection collision avoidance system for traffic signal and stop sign violation[END_REF], [START_REF] Sugimoto | Realization of head-on collision warning system at intersections-DSSS: driving safety support systems[END_REF]. Statistical studies of the causes of accidents at intersections have shown that 89% of them are due to driver error. The most common errors are perception failures, situation misunderstanding, and wrong decision [START_REF]Deliverable D2.2, Accident causation and preaccidental driving situations -In-depth accident causation analysis[END_REF].

The potential for Vehicle-to-Vehicle (V2V) technologies to address these situations is considerable. The sharing of information between vehicles over wireless links allows a vehicle to perceive its environment beyond the limits of the field of views of its on-board sensors (e.g. cameras or radars). This allows for an enlarged representation of the environment which can be used to assist the driver in perceiving and analyzing the situation. This paper focuses on the application of these technologies to the problem of situation and risk assessment at road intersections.

The problem can be formulated as follows. Each vehicle collects information over time about its own state and that of other vehicles, through proprioceptive sensors and V2V communication links. In order to assess the situation and the risk associated to it, algorithms are needed that can infer relevant information from this uncertain and incomplete information. One difficulty is that intersection scenes are highly dynamic and involve complex interactions between vehicles. Therefore physical models of the motion of vehicles are valid only for a short term and are insufficient for S. Lefèvre and C. Laugier are with Inria Grenoble Rhône-Alpes, 655 av. de l'Europe -Montbonnot, 38334 Saint Ismier Cedex, France, {stephanie.lefevre, christian.laugier}@inria.fr S. Lefèvre and J. Ibañez-Guzmán are with Renault S.A.S., 1 av. du Golf, 78288 Guyancourt, France, {stephanie.lefevre, javier.ibanez-guzman}@renault.com anticipatory risk evaluation. Instead, algorithms are needed that can reason at a higher level (e.g. maneuver intention).

In this paper we propose a probabilistic framework for reasoning about situations and risks at road intersections based on a semantic analysis of traffic situations. Vehicle motion at road intersections is modeled using a Dynamic Bayesian Network. Inference is performed on selected variables to estimate intention (i.e. what a driver intends to do) and expectation (i.e. what a driver is expected to do) jointly for all the drivers in the scene. Subsequently dangerous situations are detected by comparing the two estimates.

The remainder of this paper is organized as follows. Section 2 reviews related work. Section 3 describes our mathematical modeling of road intersection situations and our solution for risk assessment. Results from field trials using passenger vehicles equipped with off-the-shelf V2V communication links are presented and analyzed in Section 4.

II. RELATED WORK

A simple and intuitive approach is to define a set of rules which detect danger based on the context and on the current observations [START_REF] Worrall | Improving vehicle safety using context based detection of risk[END_REF]. Because context is explicitly taken into account, the characteristics of the environment can easily be incorporated. However an established limitation of rulebased systems is their inability to account for uncertainties or to reason on a high-level basis about a situation.

An alternative is to learn collision patterns from data so that dangerous configurations can be identified when they occur at a later time. Data mining techniques can be used to map the relationship between vehicles states (input) and collision risk (output) directly. A neural network was used in [START_REF] Chinea | Risk assessment algorithms based on recursive neural networks[END_REF], while the authors of [START_REF] Salim | Collision pattern modeling and real-time collision detection at road intersections[END_REF] applied the Expectation-Maximization algorithm to cluster data. Obtaining the data to learn from remains an issue, since real data is not available and simulations will not reflect real accident situations.

By far the most popular approach to collision risk estimation is the "trajectory prediction + collision detection" approach. In the first step, future trajectories are predicted for the objects in the scene using a motion model. The second step consists in checking these trajectories for intersection points. Numerous algorithms rely on a physical model of vehicles to do that [START_REF] Tan | DGPS-based vehicle-to-vehicle cooperative collision warning: engineering feasibility viewpoints[END_REF], [START_REF] Eidehall | Statistical threat assessment for general road scenes using Monte Carlo sampling[END_REF], [START_REF] Brännström | Model-based threat assessment for avoiding arbitrary vehicle collisions[END_REF], but they are not able to reason on a high-level basis about a situation and therefore are limited to short-term collision prediction. Other approaches estimate the maneuver intention of the drivers to better predict trajectories in the long term. The authors of [START_REF] Aoude | Threat assessment design for driver assistance system at intersections[END_REF], [START_REF] Aoude | Mobile agent trajectory prediction using bayesian nonparametric reachability trees[END_REF] used an evolution of the Rapidlyexploring Random Tree algorithm to generate possible future trajectories, after identifying maneuver intentions using a Support Vector Machine in [START_REF] Aoude | Threat assessment design for driver assistance system at intersections[END_REF] and Gaussian Processes in [START_REF] Aoude | Mobile agent trajectory prediction using bayesian nonparametric reachability trees[END_REF]. In [START_REF] Laugier | Probabilistic analysis of dynamic scenes and collision risks assessment to improve driving safety[END_REF] the maneuver intention is estimated using a Hierarchical Hidden Markov Model, and potential trajectories are modeled by Gaussian Processes. The first limitation of these methods is the computational cost of calculating all the possible trajectories and the pairwise probabilities that they intersect. The second one is that, in order to reduce the complexity, it is generally assumed that vehicles on different roads move independently from each other. As will be shown in Section III-A.1 and in Section IV, this assumption is not valid at intersections and affects an algorithm's ability to understand a situation and assess the risk.

III. PROPOSED APPROACH

The proposed approach focuses on intersection accidents caused by driver errors (89% of all intersection accidents). We propose to estimate jointly what each driver is expected to do in the current situation and each driver's actual intention. Risk is computed as the probability that expectation and intention do not match. The mutual influence between vehicle maneuvers is accounted for, as well as the uncertainty of the information. The following paragraphs describe the proposed model for the modeling of traffic situations at road intersections, as well as its applications to risk estimation.

A. Context-aware scene representation 1) On the importance of context: Vehicles operate on the road network and are therefore strongly constrained by its layout (geometry, topology). Besides, vehicles share the road and therefore their actions are strongly tied with other vehicles. Current approaches to situation and risk assessment for vehicles generally do not account for these facts, and assume that vehicles evolve independently from each other in a mostly unconstrained environment. Fig. 1 illustrates why this assumption, when applied to road intersection scenes, can lead to misinterpretations. The situation is the following: the dark vehicle slowed down when approaching the intersection, and is currently stopped at the entrance of the intersection. If independence between vehicles is assumed, this behavior will be (incorrectly) interpreted as an intention to make a turn. The correct interpretation is that the driver intends to yield to the white vehicle, which has priority. A more detailed study of the limitations of the independence assumption at road intersections can be found in [START_REF] Lefèvre | Modelling dynamic scenes at unsignalised road intersections[END_REF].

One challenge for representing traffic situations at road intersections is that the model should be comprehensive enough to represent complex dynamic situations such as the one described above, but inference on the relevant variables should still be tractable. With this goal in mind, the variables below are defined for a scene featuring N vehicles.

2) Physical variables (observed):

In this work the available observations are the position, heading, and speed of the vehicles. For each vehicle n ∈ N at time t, an observation variable is defined as:

O n t = (P n t S n t ), with • P n t = (X n t Y n t θ n t ) ∈ R 3 : the pose (i.e.

position and orientation w.r.t. a local frame),

• S n t ∈ R: the speed. For the experimental validation (see Section IV), this information originates from the vehicles' proprioceptive sensors and is shared via Vehicle-to-Vehicle communication (V2V). However, the method can be applied independently of the type of sensors that are used to observe the scene.

3) Behavioral variables (hidden):

The selection of the higher-level variables for representing the situation is a crucial step; a trade-off has to be found between representational power and complexity. We exploit the fact that the road network is a structured environment that constrains the motion of the vehicles, and automatically extract from the map the set of authorized maneuvers and the traffic rules (stop, give way, etc.). We use an "exemplar paths" representation, as illustrated by Fig. 2. An exemplar path is defined for each authorized maneuver as the typical path that is followed by a vehicle executing that particular maneuver in the intersection. In addition, we define a "stop point" for each exemplar path, which delimits the approaching and execution phases of a maneuver. The stop point is located at the entrance of the intersection when there is a stop or give way line, and inside the intersection for left turn across oncoming traffic maneuvers. The exemplar paths and associated stop points can be either automatically generated from the map, or learned by applying a path clustering technique to recorded data [START_REF] Morris | Learning trajectory patterns by clustering: experimental studies and comparative evaluation[END_REF], [START_REF] Hermes | Long-term vehicle motion prediction[END_REF]. For each vehicle n ∈ N at time t, a behavior variable is defined as:

B n t = (M n t D n t I n t ), with
• M n t ∈ {m i } i=1:N M : the maneuver intention of the driver, with {m i } i=1:N M the set of authorized maneu- vers at the intersection of interest. For each maneuver, an exemplar path is defined (see Fig. 2).

• D n t ∈ R: the distance traveled by the vehicle along the exemplar path of M n t , i.e. the curvilinear abscissa. D n t is negative when the vehicle has not yet reached the stop point, and positive after the vehicle passed the stop point.

• I n t ∈ {0, 1}: the driver's intention to stop at the intersection (= intention).

4) Expectation variable (hidden):

For each vehicle n ∈ N , the relevant traffic rules at time t are incorporated into the variable:

• E n t ∈ {0, 1}: whether or not the driver is expected to stop at the intersection (= expectation).

For more clarity in the equations, in the remaining of this paper factored states will be used to represent the states of all the vehicles in the scene, e.g. M t = (M 1 t ...M N t ) (and similarly for all the variables defined above). We argue that (M t D t I t E t ) is a relevant high-level representation of a road intersection traffic situation, since inference on these variables allow to estimate key features of the situation. The remaining of this section introduces the proposed model for linking behavioral, observation and expectation variables, and how it can be used for risk assessment.

B. Dynamic scene modeling

In the following paragraphs, a Dynamic Bayesian Network (DBN) is described that links together the variables defined above in order to model the evolution of a traffic situation at a road intersection. The use of the Bayesian formalism allows us to take into account uncertainties on the relationships between the variables.

1) Joint distribution:

The graphical representation of the DBN is shown for one vehicle in Fig. 3. For all the vehicles, the following local joint distribution can be derived:

× N n=0 P (E n t |Bt-1Ot-1) × P (B n t |B n t-1 E n t ) × P (O n t |O n t-1 B n t )
Further, we assume the following independencies:

• For the expectation variable:

P (E n t |Bt-1Ot-1) = P (E n t |Mt-1Dt-1St-1)
• For the behavioral variables:

P (B n t |B n t-1 E n t ) = P (M n t |M n t-1 ) × P (I n t |I n t-1 E n t ) × P (D n t |M n t-1 D n t-1 I n t-1 )
• For the physical variables:

P (O n t |O n t-1 B n t ) = P (P n t |M n t D n t ) × P (S n t |S n t-1 M n t D n t I n t )
By making I n t dependent of E n t , which itself depends on the other vehicles, we take into account the mutual influences between the maneuvers performed by the vehicles in the scene. The parametric form of the conditional probability terms are described below, along with the hypotheses they build on.

2) Pose:

The likelihood of a pose while executing maneuver m and being at distance d from the stop point is defined as a bivariate normal distribution with no correlation between the position and the orientation:

P (P n t |[M n t = m][D n t = d]) = 1 2πσ δ σ α × e -1 2 ( δ 2 σ 2 δ + α 2 σ 2 α )
where δ is the distance between the vehicle's position (x, y) and the point (x ′ , y ′ ) with curvilinear abscissa d on the exemplar path of maneuver m, α is the angle between the vehicle's orientation θ and the orientation of the exemplar path at point (x ′ , y ′ ). σ δ (resp. σ α ) is the standard deviation set for the distance (resp. the angle). They are set according to the precision of the pose sensor.

3) Speed: It is assumed that drivers adapt their speed to their intentions and to the geometry of the road. The distribution on S n t is normal and defined as:

P (S n t |[S n t-1 = s][M n t = m][D n t = d][I n t = i]) = N (µS, σS)
with the calculation of µ S detailed below.

A number of statistical analyses of the behavior of drivers approaching an intersection can be found in the literature, e.g. [START_REF] Berndt | Driver braking behavior during intersection approaches and implications for warning strategies for driver assistant systems[END_REF]. From these it is possible to derive generic speed Fig. 4: Example average (s A ) and maximum (s M ) speed profiles generated for a pair (M n t = m, I n t = 1) (left) and for a pair (M n t = m, I n t = 0) (right). When I n t = 0, the lowest speed in the profile is a function of the curvature of the path associated with the maneuver m.

profiles for vehicles negotiating an intersection. We define s A = f (d) the average speed profile at the intersection and s M = f (d) the maximum speed profile (i.e. the highest speed at which it is possible to negotiate the intersection). Additionally we take into account the geometry of the road: for each possible pair (M n t , I n t ), the generic speed profiles are adapted to match the constraints imposed by the curvature of the exemplar path. In Fig. 4 the general aspect of the speed profiles is shown. They serve as a basis for predicting the evolution of the speed of a vehicle given the driver's intention, following the equation:

µ S = s A (d) - s A (d) -s M (d) s A (d ′ ) -s M (d ′ ) × (s A (d ′ ) -s)
with d ′ the distance at the previous timestep calculated as

d ′ = d -s × ∆t.

4) Distance to stop point:

The distribution on D n t is normal and defined as:

P (D n t |M n t-1 D n t-1 I n t-1 ) = N (µ D , σ D )
where µ D and σ D are extracted from the speed profiles described in III-B.3.

5) Maneuver intention:

Continuity in the maneuver intention is assumed.

P (M n t |M n t-1 ) = 0.9 if M n t = M n t-1 0.1 N M -1 otherwise
The value 0.9 was set manually and matches the authors' interpretation of "continuity in the maneuver intention", but should ideally be learned from data.

6) Expectation to stop:

It is assumed that the necessity for a driver to stop at an intersection is a consequence of the context at time t -1 (priority rules, presence of other vehicles). The necessity for a vehicle to stop given the context is derived using probabilistic gap acceptance models found in [START_REF] Pollatschek | A decision model for gap acceptance and capacity at intersections[END_REF], [START_REF] Chan | Characterization of driving behaviors based on field observation of intersection left-turn across-path scenarios[END_REF]. If we take as an example a vehicle v n heading towards a give way intersection, the calculation is: i. Project forward the position of v n until the time t n when it reaches the stop point of the maneuver. A constant speed model is used, with a combination of the vehicle's current speed and the average speed profile s A of the maneuver. ii. For each vehicle v m whose maneuver is prioritary w.r.t.

the maneuver of v n , project forward the position of v m until the time t m when it reaches the stop point of the maneuver (following the same procedure as for v n ). iii. Select the smallest positive time gap available for v n to execute its maneuver: t min = min m (t mt n ), f or t mt n ≥ 0 iv. The necessity for v n to stop at the intersection is calculated as the probability p stop that the gap t min is not sufficient, using a probabilistic gap acceptance model (from [START_REF] Pollatschek | A decision model for gap acceptance and capacity at intersections[END_REF] for merging cases, from [START_REF] Chan | Characterization of driving behaviors based on field observation of intersection left-turn across-path scenarios[END_REF] for left turn across oncoming traffic cases):

P ([E n t = 0]|M t-1 D t-1 S t-1 ) = 1 -p stop P ([E n t = 1]|M t-1 D t-1 S t-1 ) = p stop
This context-aware reasoning about the necessity for a vehicle to stop at the intersection will allow us to detect vehicles running stop signs, or vehicles entering an intersection when they should have waited for another vehicle to pass. A similar calculation can be done for intersections ruled by traffic lights, but this is not the focus of this work.

7) Intention to stop:

The evolution model for the driver's intention to stop is based on the comparison between I n t-1 and E n t . If the driver's intention at time t -1 coincides with what is currently expected of him, it is assumed that chances are high that the driver will comply. Otherwise a uniform prior (0.5) is assumed.

P ([I n t = 1]|I n t-1 E n t ) =      0.9 if I n t-1 = E n t = 1 0.1 if I n t-1 = E n t = 0 0.5 if I n t-1 = E n t
The values 0.9 and 0.1 were set manually and match the authors' interpretation of "chances are high that the driver will comply", but should ideally be learned from data. Results (see Section IV) show that this simple modeling of the interactions between vehicles leads to a better situation and risk estimation, compared with a constant uniform prior on I n t .

C. Bayesian risk assessment

Using the DBN described above, from the successive observations (pose and speed) it is possible to infer the intentions of the drivers as well as what they are expected to do. As an alternative to the conventional "trajectory prediction + collision check" approach to risk estimation, we propose to base the computation of the risk on the probability that expectation and intention do not match, i.e.:

P ([I n t = 0][E n t = 1]|P 0:t S 0:t ) (1)
In this work inference was performed using a particle filter. The results presented in Section IV were obtained with 400 particles. One advantage of this approach resides in its flexibility in terms of applications. An example of a safety-oriented application is the detection of hazardous vehicles: the system can compute a "hazard probability" for every vehicle in the scene using Eq. 1 and warn all the drivers in the intersection area when the probability is higher than a predefined threshold. Alternatively the model can be used to compute the risk of a specific maneuver for a vehicle, which is an important feature for autonomous driving. Another interesting aspect of the model is that it is predictive, therefore the future states of the vehicles can be estimated. The applications targeted by this paper do not exploit this feature, but it is relevant to numerous applications.

IV. EXPERIMENTAL EVALUATION

The objective is to evaluate the ability of the algorithm to assess the risk of a situation. The choice of an evaluation method is not trivial, since there exists no ground truth for the quantity "risk of a situation". A solution is to perform a comparative study between two algorithms, comparing the issuance time of warnings in dangerous situations and the frequency of false alarms in non-dangerous situations. Another relevant comparison is their ability to assess situations correctly, e.g. the maneuver intention of the drivers in the scene, for which the ground truth M GT is available. We selected two versions of our algorithm to be compared :

• Proposed model: corresponds to the model described in Section III. Interactions between vehicles are accounted for in P (I n t |I n t-1 E n t ). • Reference model: corresponds to a modified version the model described in Section III. Independence is assumed between the vehicles in the scene by replacing P (I n t |I n t-1 E n t ) by a uniform prior P (I n t ). Comparing these two models will allow us to evaluate the benefits of taking into account vehicle interactions when performing situation and risk assessment at road intersections.

A. Experimental setup

Two passenger vehicles were equipped with off-the-shelf V2V communication modems (802.11p) and shared their pose and speed information at a rate of 10 Hz. In each vehicle the pose information was obtained via a GPS + IMU unit with a precision of σ = 2m for the position. The CAN provided the speed information. In its current non-optimized state the algorithm runs at 10 Hz on a dedicated dual core 2.26 GHz processor PC.

Experiments were conducted at a give-way intersection. Five scenarios were defined, involving a Subject Vehicle (SV) driving on the main road and an Other Vehicle (OV) performing various maneuvers. They are illustrated in Fig. 5.

• Scenario 1: OV approaches on the secondary road and makes a right turn without stopping.

• Scenario 2: OV approaches on the secondary road, and proceeds to make a left turn when SV is about to reach the intersection. • Scenario 3: OV approaches on the secondary road, stops at the give way line, waits until SV has passed and makes a left turn. • Scenario 4: OV approaches on the main road from the opposite direction, stops in the intersection, and proceeds to make a left turn when SV is about to reach the intersection. • Scenario 5: OV approaches on the main road from the opposite direction, stops in the intersection, waits until SV has passed and makes a left turn. In the rest of this section we consider an application whose goal is to issue a warning whenever it detects that a situation is dangerous. A warning is triggered iff:

∃n ∈ N : P ([I n t = 0][E n t = 1]
|P 0:t S 0:t ) > λ Scenarios 2 and 4 are dangerous scenarios where emergency braking is the only way to avoid an accident, therefore we expect the application to issue a warning as early as possible. Scenarios 1, 3 and 5 are non-dangerous scenarios where a warning would be considered a false alarm.

In total ninety trials were carried out, with five different drivers for OV. In order to generate some variations in the scenarios, the drivers of OV were not given clear instructions about the execution of the various maneuvers and were only told to create what they felt were dangerous and nondangerous situations.

B. Qualitative results

In order to appreciate and understand the performance differences between the two models, it is useful to analyze the results on sample instances. To this end we display in Fig. 6 the estimation results for the maneuver intention of the drivers and for the probability that OV is dangerous in example instances of Scenarios 1, 2 and 3.

1) Scenario 1: In this scenario OV can execute his maneuver without considering SV since the vehicles always stay on different lanes. t 1t 2 : the model which assumes independence between the vehicles has no way to tell whether the intention of OV's driver is to turn right or left when he enters the intersection. Since executing a left turn maneuver at this time is clearly a violation of the priority rules, the probability of danger peaks at 0.5 for a short period before resuming to a lower value. The approach which models the influences between vehicle maneuvers favors the right turn as the most likely intention of the driver of OV, and as a result the probability of danger always remains lower than 0.2. This difference has a direct consequence on the implementation of the application. The threshold λ has to be set with the double objective to issue a warning as early as possible in dangerous situations and to keep the rate of false alarms low. In the case of the Reference model, λ has to be set to a value higher than 0.5 in order to avoid systematic false alarms in This means that the sensitivity of the risk estimation [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF] is significantly higher with the Proposed approach, since in dangerous situations the value of the risk varies in the range [0.3 1.0], against [0.65 1.0] with the Reference approach.

2) Scenario 2: In this scenario OV violates the priority, and without an emergency braking the situation would have resulted in a crash. Using the values of λ defined in the previous paragraph, we can see that there is no difference between the two models on the time of issuance of the warning (t 3 on the graph).

3) Scenario 3: In this scenario both vehicles behave as expected given the priority rules and there is no danger. Neither of the models produces a false alarm, however there is a noticeable difference in their ability to assess the situation correctly. Indeed, in the interval t 4 -t 5 the Proposed approach is able to interpret the behavior of OV (waiting at the entrance of the intersection while SV proceeds) as an indication that the driver's intention is to turn left. The Reference approach does not have this ability, since it assumes independence between the vehicles.

C. Statistical results

We compute statistical results over the ninety test instances in order to verify if the differences observed on specific instances (see previous paragraph) can be generalized. The performances of the two models are measured by looking at their ability to:

• Correctly estimate the maneuver intention of the drivers: this is measured by P (M t = M GT ). • Correctly assess the situation in general: this will be reflected by N ef f , the effective sample size of the particle filter (divided by the number of particles so that a relative value is obtained). • Correctly assess the risk. In the case of a dangerous situation, this means issuing a warning at time t warn as early as possible. In the case of a non-dangerous situation the risk should not be overestimated: this feature can be evaluated by monitoring ∆P ([ results obtained by the Proposed model. They are displayed for dangerous and non-dangerous situations in Fig. 7 and commented below.

I t = 0][E t = 1]).
1) Maneuver intention estimation: For dangerous situations, the two approaches perform similarly on average. When there is no danger, the two approaches are able to estimate the maneuver intention equally well 55% of the time. However the Proposed approach leads to a better estimation on average and performs better 40% of the time.

2) Neff: Accounting for the interactions between the vehicles has a positive and significant impact on N ef f both when the situation is dangerous and when it is not. In both cases the performance is increased 60% of the time and mostly equivalent the rest of the time. These results validate our modeling of vehicle interactions since the value of N ef f reflects the quality of the tracking, which itself is dependent on how well the model matches the reality.

3) Risk estimation: For every instance of a dangerous situation during the field trials, both approaches were able to issue a warning early enough that the collision was avoided by braking. In 55% of the cases the warning was issued at the exact the same time by the two approaches. When the situation was not dangerous, neither models produced any false alarm during the experiments. However, the plot of ∆P ([I t = 0][E t = 1]) shows that the risk is systematically higher when estimated using the Reference model.

V. CONCLUSIONS AND FUTURE WORK

This paper presented a novel framework for reasoning about situations and risk at road intersections. The risk of a situation is assessed based on the comparison between what drivers intend to do and what they are expected to do. This intuitive formulation of risk provides a lot of flexibility for safety applications and is relevant to both ADAS and autonomous driving.

The approach was validated by field experiments in the context of communicating vehicles, using passenger cars and off-the-shelf V2V communication links. The results demonstrated the ability of the algorithm to issue a warning in dangerous situations, and the benefits of taking into account interactions between the vehicles when reasoning about situations and risk at road intersections.

In its current form the method may be applied to any intersection layout and any number of vehicles, but this was not demonstrated in this paper and will be evaluated in future work. The performance of the algorithm should also be compared with that of state-of-the-art works. Further, the estimation of drivers' intentions could be improved by taking into account information about drivers' actions such as steering angle and pedal pressure.

Fig. 1 :

 1 Fig. 1: Illustration of the importance of taking context into account.

Fig. 2 :

 2 Fig. 2: Representation of a road intersection: the exemplar paths (black arrows) and associated stop points (red dots) originating from one road are displayed.

Fig. 3 :

 3 Fig. 3: Graphical representation of the DBN for one vehicle. Bold arcs correspond to multi-vehicle dependencies.

Fig. 5 :

 5 Fig. 5: Test scenarios.

Fig. 6 :

 6 Fig. 6: Estimated maneuver intention of OV (left) and danger introduced by OV (right) in a sample instance of Scenario 1 (top), 2 (middle) and 3 (bottom).

Fig. 7 :

 7 Fig. 7: Comparative results for dangerous (top) and nondangerous (bottom) situations.

ACKNOWLEDGEMENT

The authors would like to thank D. Vasquez and E. Mazer for useful discussions about the model.