Jérémie Allard

Franc ¸ois Faure

Hadrien Courtecuisse

Florent Falipou

Christian Duriez

Paul G Kry

Volume Contact Constraints at Arbitrary Resolution

Keywords: CR Categories: I.3.5 [Computer Graphics]: Physically based modeling-[I.3.7]: Computer Graphics-Animation physically based animation, contact forces, Coulomb friction, constraints

Figure 1: Left: A stack of tori, with the intersection volumes and normals on the GPU. Right: intersection volume constraints allow us to resolve frictional contact in challenging scenarios.

Introduction

Dealing with contact is a fundamental problem in physically based computer animation, and possibly one of the most challenging. Accurate simulation of frictional contact is important for a wide variety of natural phenomena and has numerous applications, from video games, to robotics. c ACM, 2010. This is the author's version of the work. It is posted here by permission of ACM for your personal use. Not for redistribution. The definitive version will be published in ACM Transactions on Graphics 29(3) (SIGGRAPH 2010).

The computation of contact forces involves several important difficulties. One such problem is collision detection and modeling. This requires not only detection of the intersection of two objects, which is a complex problem on its own, but also computation of the depth and direction. While the problem can be straightforward when geometric meshes are in close proximity, or in point contact, the solution is less clear in cases that involve deep interpenetration. This imposes restrictions on the time step size, or obliges the use of continuous collision detection methods that permit the processing of collision at the exact time of contact. For interactive simulation, however, it can be preferable to use methods that can handle interpenetrations gracefully, without additional computational cost.

Another difficulty in contact force computation is the solution of the contact equations themselves. Solving contact with Coulomb friction involves Lagrange multipliers and coupled inequalities. With the inclusion of dry friction, the problem is larger and particularly complex because we must determine which constraints are active and this can be a hard combinatorial problem. Furthermore, numerous simultaneous contacts can lead to constraints that may be redundant, or even inconsistent, which results in numerical problems. Processing contacts sequentially alleviates some difficulties but may create artifacts and raises convergence issues. Despite decades of research and significant progress, the computation of contact forces between complex geometries remains a hard problem. Interactive applications use simplified collision geometries such as spheres, cylinders, and cubes, which restricts the range of simulated objects or leads to visible artifacts. The simulation of contact between fine meshes is generally performed offline.

In this paper, we propose a completely new approach for contact force computation, suitable for all objects bounded by triangular meshes. The collision detection and modeling phase returns the size of the intersection volume and its derivative with respect to the vertices of the meshes. In contrast to traditional methods, this allows the non-interpenetration constraint to be modeled using a single scalar equation, i.e., the size of the intersection volume, rather than a (potentially large) number of distance constraints between geometric primitives. The result is a dramatically simpler system of equations. Figure 1 shows a visualization of this and a preview of our results. Effectively, this approach can be seen as using contact volumes as opposed to contact points to formulate the problem. Additionally, we show how to constrain the relative tangential motion, allowing frictional contact between arbitrarily complex objects. To allow simultaneous sticking and sliding behavior across the con-tacts between objects, we extend our geometric model by splitting the single contact volume into several contacts. We use a regular grid where the contact equations are generated from the cells that contain parts of the intersection volume. Tuning the resolution of the grid allows a trade off between accuracy and speed. To the best of our knowledge, this is the first method where the complexity of the contact model is completely independent from the complexity of the surface model without pre-processing.

Our contact equations are purely geometric and can be applied to arbitrary Lagrangian models including rigid and deformable solids. They can be solved using the traditional numerical methods. Moreover, no precomputation is needed, which makes the method suitable for animations with topology changes such as cutting, fracturing, or adaptive levels of detail. See Figure 2 for an example of this and other difficult scenarios that our method handles efficiently. Finally, a major advantage of our method is that it is no longer necessary to worry about collision proxies, signed distance fields, or other adapted collision models; collision processing and contact force computation can easily be done with the same highly detailed geometry used for rendering.

The remainder of this paper is organized as follows. We review related work in Section 2 and provide the necessary background in Section 3. Our new constrained volume method is presented in Section 4. This is extended to multiple volumes in Section 5. Results are presented and discussed in Section 6.

Related Work

In most previous work, the problem of collision detection and response are addressed separately. The amount of work done on collision detection alone is quite vast. Collision detection typically involves hierarchical bounding volume data structures, such as bounding boxes [START_REF] Gottschalk | OBB-Tree: a hierarchical structure for rapid interference detection[END_REF] or spheres [START_REF] Hubbard | Collision Detection for Interactive Graphics Applications[END_REF]]. In addition to rigid object collision, significant work has been done to address deformable models [van den Bergen 1997; [START_REF] James | BD-tree: output-sensitive collision detection for reduced deformable models[END_REF], and self-contact [START_REF] Volino | Collision and self-collision detection: Efficient and robust solutions for higly deformable surfaces[END_REF][START_REF] Provot | Collision and self-collision handling in cloth model dedicated to design garments[END_REF]]. See the survey by [START_REF] Teschner | Collision detection for deformable objects[END_REF] for an excellent overview of different methods. More recently, techniques which exploit the GPU have been proposed as this can result in important speed improvements through parallelization [START_REF] Vassilev | Fast cloth animation on walking avatars[END_REF][START_REF] Heidelberger | Real-time volumetric intersections of deforming objects[END_REF][START_REF] Heidelberger | Detection of collisions and self-collisions using image-space techniques[END_REF][START_REF] Baciu | Image-based collision detection for deformable cloth models[END_REF][START_REF] Wong | GPU-based intrinsic collision detection for deformable surfaces: Collision detection and deformable objects[END_REF][START_REF] Sud | Fast proximity computation among deformable models using discrete voronoi diagrams[END_REF]]. Most of the techniques return pairs of primitives, between which distance constraints can be formulated to resolve the contacts. [START_REF] Otaduy | Haptic display of interaction between textured models[END_REF] use the GPU to compute penetration depth in a given direction. Our work builds on that of Heidelberger et al. extended by Faure et al. [2008], where the GPU is used to return the volume of intersection and self-intersection between the objects, along with the derivatives of these volumes. This new type of collision data allows a new approach of collision response. Our work significantly improves upon that of Faure et al. in that we show how to formulate unilateral contact constraints instead of penalty forces, Coulomb friction instead of viscous friction, and we introduce an arbitrary-resolution formulation for resolving varying interaction forces across contact patches.

While it is possible to solve for contact response using penalty methods, it has long been recognized that there are advantages to formulating contact with constraints. For contact, the constraints are unilateral, which results in a linear complementarity problem (LCP). The addition of Coulomb friction involves a non-linear constraint relating the tangential force to the normal force. Most work in graphics uses a pyramid discretization of the friction cone to formulate the problem as an LCP [START_REF] Baraff | Coping with friction for non-penetrating rigid body simulation[END_REF][START_REF] Milenkovic | Optimizationbased animation[END_REF]], but it is also possible to use the exact cone (a non-linear complementarity problem) and to compute the solution using iterative methods [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF]]. While natural to formulate the frictional contact equations at the acceleration level, a velocity level formulation is preferred since it avoids the Painlevé paradox [START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF][START_REF] Stewart | An implicit timestepping scheme for rigid body dynamics with inelastic collisions and coulomb friction[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF]. However, the problem can still result in an NP-hard combinatorial problem to identify the active constraints [START_REF] Baraff | Coping with friction for non-penetrating rigid body simulation[END_REF]]. While some recent work strives to find solutions which are as correct as possible [START_REF] Harmon | Asynchronous contact mechanics[END_REF][START_REF] Kaufman | Staggered projections for frictional contact in multibody systems[END_REF]], we instead focus on an approximate solution that is fast and interactive. We exploit small interpenetrations to formulate constraints, which also makes our approach naturally robust to larger interpenetrations. Velocity level implicit integration is typical in many works in graphics because it can also lead to stable simulations with large time steps [START_REF] Baraff | Large steps in cloth simulation[END_REF]]. [START_REF] Otaduy | Implicit contact handling for deformable objects[END_REF] use velocity level implicit contact constraints in combination with iterative methods to simulate combined elastic, rigid, and thin shells. Velocity level LCPs have also been used for fluid-solid coupling [START_REF] Batty | A fast variational framework for accurate solid-fluid coupling[END_REF]]. Other relevant work on friction includes fast approximate models in video game systems [START_REF] Parker | Real-time deformation and fracture in a game environment[END_REF], and simulation of anisotoripic friction [START_REF] Pabst | Anisotropic friction for deformable surfaces and solids[END_REF]].

In graphics, the methods used to solve contact equations can be classified into two categories: iterative methods that repeatedly process the contacts sequentially [START_REF] Erleben | Velocity-based shock propagation for multibody dynamics animation[END_REF][START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF][START_REF] Otaduy | Implicit contact handling for deformable objects[END_REF], and direct methods that build and solve a single system of equations [START_REF] Baraff | Coping with friction for non-penetrating rigid body simulation[END_REF][START_REF] Baraff | Fast contact force computation for nonpenetrating rigid bodies[END_REF][START_REF] Pauly | Quasi-rigid objects in contact[END_REF]. Kaufmann et al. [2008] use a direct solver, but alternatively solve nonpenetration and friction constraints to avoid scaling and non-convexity problems associated with traditional direct methods, while improving on the slow convergence that can be typical with iterative methods. Alternatively, for frictional dynamics, Lemke's direct method can be modified to run faster by introducing frictional complementary constraints only when necessary [START_REF] Lloyd | Fast implementation of Lemke's algorithm for rigid body contact simulation[END_REF]]. But much of the recent work in graphics has shifted to iterative methods because the implementation of such methods is less complex and the computation time is easier to control. Our method can be applied to both approaches. In this work we use iterative methods, and our arbitrary-resolution frictional contact formulation helps address the problem of slow convergence by decoupling the resolution of the contact from that of the geometry. Models at multiple resolutions have been exploited in graphics for simulation of elastic materials and deformation of models [START_REF] Debunne | Dynamic real-time deformations using space and time adaptive sampling[END_REF][START_REF] Otaduy | Adaptive deformations with fast tight bounds[END_REF]Barbič and James 2007;[START_REF] Shi | A fast multigrid algorithm for mesh deformation[END_REF]Müller 2008]. The main contributions in previous work have been irregular nesting, haptic rates for frictionless contact, and speed improvements through the use of multigrid. Our work takes the unique approach of treating the contact constraints at an adjustable resolution, independent of the degrees of freedom. While there has been vast amounts of work on physically based deformation (see [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF]] for a good survey), the focus is typically deformation as opposed to contact.

Numerical integration of differential equations with position level constraints requires stabilization to avoid numerical drift [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF]]. If small steps are taken then drift is not always a significant problem and can be left as a separate problem [START_REF] Kaufman | Fast frictional dynamics for rigid bodies[END_REF][START_REF] Kaufman | Staggered projections for frictional contact in multibody systems[END_REF]. Gundleman et al. [] provide a bottom up approach to deal with stacking. A related problem is the robust treatment of contact resolution [START_REF] Bridson | Robust treatment of collisions, contact and friction for cloth animation[END_REF][START_REF] Harmon | Robust treatment of simultaneous collisions[END_REF]. Our work uses larger step sizes, and both allows and exploits interpenetration, thus, we also include stabilization. [START_REF] Baumgarte | Stabilization of constraints and integrals of motion in dynamical systems[END_REF] introduced a popular stabilization method, which is widely used, but also widely criticized for the difficulty of tuning parameters. Post stabilization is a preferable alternative [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF]] and is straightforward to implement as the constraint gradient is available from the formulation of the constrained equations of motion. However, unilateral constraints should not be stabilized in the same manner as bilateral constraints. A separate LCP can be formed to solve the post step [START_REF] Cline | Post-stabilization for rigid body simulation with contact and constraints[END_REF], or in some cases the time stepping and unilateral post step can be combined into a single LCP [START_REF] Anitescu | Constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]. While our formulation can accommodate any of these approaches to stabilization, our solution most closely resembles that of [START_REF] Cline | Post-stabilization for rigid body simulation with contact and constraints[END_REF] as we solve an iterative LCP to perform the position correction step.

Background

In this section we introduce the necessary background material on which we build our method. This consists of volume and volumegradient computation from layered depth images (LDIs), and a review of constrained dynamics. We consider a physical system governed by the ordinary differential equation

Mq = f (q, q) + fex (1)
where vector q is the independent degrees of freedom (DOFs), M is the mass Matrix, f provides the internal forces, and fex are the external forces (such as those due to contact and friction).

Numerical Integration and Contact Constraints

Overall, the methods that we describe here have also been used or described in depth with small variations by many others [START_REF] Baraff | Large steps in cloth simulation[END_REF][START_REF] Erleben | Velocity-based shock propagation for multibody dynamics animation[END_REF][START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF][START_REF] Otaduy | Implicit contact handling for deformable objects[END_REF][START_REF] Kaufman | Staggered projections for frictional contact in multibody systems[END_REF][START_REF] Pabst | Anisotropic friction for deformable surfaces and solids[END_REF][START_REF] Parker | Real-time deformation and fracture in a game environment[END_REF].

Integration

We use an implicit integration method to solve the ODE in Equation 1 because our system may involve stiff elastic deformation (the velocity level formulation will also be of benefit for friction). To reduce computation time, we use a backward Euler scheme with one Newton iteration to solve for the state update [START_REF] Baraff | Large steps in cloth simulation[END_REF]]. This works well as the internal forces are often smooth, but we note that other numerical integration methods with multiple Newton iterations could also be used for better accuracy at additional cost.

Given a time step size h and the current state (q0, q0), we compute an implicit velocity update ∆ q by solving the linearized system

A∆ q = b, (2)
where A = M-hB-h 2 K, and b = hf (q0, q0)+h 2 K q0. Recall, B = ∂f ∂ q and K = ∂f ∂q are the damping and stiffness matrices. From the solution, the velocities at the next time step are computed as q0+h = q0 + ∆ q.

(3)

These new velocities are subsequently used to perform an implicit update on the position, i.e., q 0+h = q0 + h q0+h , or using the exponential map for rigid motion.

Contact When point contacts exist in the system, the contacts form separation distance constraints, g(q) ≥ 0. The transpose of the constraint gradient J = ∂g/∂q provides the contact force directions (i.e., contact normals), which we add to the system using Lagrange multipliers λ. Thus, Equation 2becomes

A∆ q = b + J T λ. (4)
The solution of this system must respect the Signorini condition.

That is, contacts do not produce attraction forces, λ ≥ 0, separation distance at contacts must remain non-negative, g(q) ≥ 0, and the constraint is active (λ nonzero) if and only if the separation distance is zero. We write this as 0 ≤ g(q 0+h ) ⊥ λ ≥ 0, and note that we must differentiate the constraint to produce complementarity conditions involving only the unknown velocity update and Lagrange multipliers,

0 ≤ J( q0 + ∆ q) ⊥ λ ≥ 0. (5)
The combination of Equations 4 and 5 produces a mixed LCP (MLCP), which can be solved with either direct or iterative methods (see discussion in Section 2).

Stabilization Discrete collision detection can not anticipate the apparition of new contacts, and numerical integration errors can produce drift. Thus, a post stabilization position update is computed at the end of each numerical integration step. However, the position update ∆q must also respect the complementarity conditions. Either there is a separation at a contact, g(q 0+h ) > 0, or the constraint provides a repulsive force, ζ > 0, but not both. This produces the MLCP

A∆q = J T ζ, (6) 0 ≤ g(q0) + J∆q ⊥ ζ ≥ 0, (7)
which is solved for ζ and ∆q. When g is smooth, one iteration of this linearized post step is sufficient to produce a good correction.

Friction Coulomb friction involves additional tangential forces at the contact points. These are typically written in a symmetric basis of 4 vectors at each contact point (negative and positive t1 and t2 directions for each contact coordinate frame {n, t1, t2}), which lets the friction force be written as a combination involving only positive weights. For contact i, let 0 ≤ βi ∈ R 4 be the coordinates of the friction force and T T i be the basis that applies the force on the DOFs in equal and opposite directions. Equation 4becomes

A∆ q = b + J T λ + T T i βi, (8)
but it is now subject to additional complementarity conditions. Specifically, there is either tangential sliding at the contact, or the contact force is not on the boundary of the friction cone, but not both. For contact i, the pyramidal friction cone conditions can be written

0 ≤ σi ⊥ µλi -eβi ≥ 0, (9) 0 ≤ βi ⊥ T T i ( q0 + ∆ q) -e T σi ≥ 0, (10)
where e = (1, 1, 1, 1), and the slack variable σi is nonzero when there is sliding at the contact. One only needs the contact normals to build a tangent space basis to set up these frictional contact equations. Note also that the friction constraint is not considered in the post stabilization step.

Volume-Based Penalty Force

Let us briefly review the image-based collision response method we extend in this work. The basic idea is to minimize the intersection volume between two polyhedra. The intersection volume is computed on the GPU based on bounding pixels in the rasterization of the object surfaces into Layered Depth Images (LDIs) [START_REF] Heidelberger | Real-time volumetric intersections of deforming objects[END_REF][START_REF] Heidelberger | Detection of collisions and self-collisions using image-space techniques[END_REF]]. This data structure uses a stack of images to represent the object sufaces as discrete height fields. The necessary number of images depends on the number of surface layers in the chosen viewing direction. Each pixel stores the surface depth in the viewing direction, as well as additional data such as normal orientation or object index. By sorting the depths stored in the different images at each pixel (i, j), we obtain the ordered list of surface intersections with a ray parallel to the viewing direction. Along the ray, each object volume is represented by one or several depth intervals between entry points and exit points. Collision detection is straightforward based on interval intersections. One LDI in an arbitrary viewing direction is sufficient to detect collisions between volumetric objects. Computing the gradient of the intersection volume requires three LDIs in mutually orthogonal directions [START_REF] Faure | Image-based collision detection and response between arbitrary volume objects[END_REF]. At left in Figure 1 is a depiction of LDI volume models using red, green, and blue pixels to denote x, y, and z viewing directions, respectively. The corresponding intersection volumes are shown in the second image of the same figure. Figure 3 represents a slice of an intersection volume. The rasterization can be done in any direction, but for simplicity we assume an orthogonal projection along one of the primary axes. Thus, the volume computed using a z projection is

V = a (i,j)∈S + z z + ij -a (i,j)∈S - z z - ij , ( 11
)
where a is the area of a pixel, z + ij and z - ij are the upper and lower pixel depths, and the sets S + z and S - z respectively contain the pixel locations (i, j) of the upper and lower contact surfaces.

The derivative of the volume with respect to the depth of a given triangle vertex k represents the variation of the volume size corresponding to a unit displacement of the vertex in the viewing direction,

∂V ∂p z k = a (i,j)∈S + z ∂z + ij ∂p z k -a (i,j)∈S - z ∂z - ij ∂p z k . ( 12
)
where p z k is the z coordinate of vertex number k. Conveniently, the scalar ∂zij(p)/∂p z k at a given pixel corresponds to the barycentric coefficient used to interpolate the depth value from the depth at vertex k (that is, the value is readily available in the GPU during rasterization). If the vertex p k is not part of the triangle containing the pixel then the partial derivative is zero. For most vertices, the gradient receives a contribution from at most one contact surface, either S + z or S - z , but some vertices receive contributions of both, such as p3 in Figure 3.

Given only a projection in z, it is not possible to use the barycentric weights to accurately compute the other partials, ∂zij(p)/∂p x k and ∂zij(p)/∂p y k . Thus, we use LDIs with projections along the other two axes to compute derivatives with respect to the other vertex coordinates. For each axis, the sums over pixels for the volume and gradient computation are done simultaneously. This is done three times for the three LDIs. That is, the volume is accumulated three times (and subsequently corrected by 1 3 ), while its derivative with respect to the coordinate of each vertex is accumulated into the gradient vector,

∂V ∂p = ∂V ∂p1 ... ∂V ∂pn , (13)
∂V ∂p k = ∂V ∂p x k ∂V ∂p y k ∂V ∂p z k , (14)
where x, y, z are the three successive LDI viewing directions. [START_REF] Faure | Image-based collision detection and response between arbitrary volume objects[END_REF] apply a repulsion force derived from a volumebased potential, E = 1 2 kV 2 , i.e., a soft constraint. Thus, the force on vertex pi is the transpose of

- ∂E ∂pi = -kV ∂V ∂pi (15)
In this paper, we exploit the volume and its gradient differently, but they remain the only necessary geometric values.

When the independent DOFs q are not directly the vertex positions p, as in the case of rigid bodies or detailed surfaces embedded in coarse deformable models, we use the chain rule to differentiate the intersection volume with respect to the DOFs: ∂V ∂q = ∂V ∂p ∂p ∂q .

Volume-Based Contact Constraints

In this section, we present a new way of handling contact and friction by defining constraints based on the intersection volume and volume derivative. We first reformulate the Signorini condition in order to obtain an equivalent complementarity problem between contact reaction pressures and interpenetration volume. This volume-based formulation can directly use the LDI collision output to define a very efficient constraint-based collision response. In Section 4.2 we address the subtleties of also adapting the Coulomb friction complementarity conditions to the volume-based formulation. In Section 4.3 we describe how these volume-based constraints are used in the main animation loop.

Signorini Condition for Volume Contacts

For simplicity, we consider a system consisting of two colliding objects, where the intersection volume V(q) is computed as described in Section 3.2. As stated in the background, a common way of dealing with contact between two objects is to build a complementarity problem between separation distances g and reaction force magnitudes λ. If two deformable objects are in contact over an area, then this involves using multiple contact constraints to track the separation distances across the contact patch.

Equivalence Recall that each unilateral point-contact constraint has conditions 0 ≤ gi(q 0+h ) ⊥ λi ≥ 0. The solution of the constrained system must respect the principle of virtual work. That is, constraint satisfaction must not add or remove energy from the system. Since we are solving the system at a velocity level, this can be written λi • Ji( q0 + ∆ q) = 0, where Ji = ∂g i ∂q is the separation distance gradient. However, we can associate a small area Ai with each contact i, and observe that

λi 1 Ai • AiJi( q0 + ∆ q) = 0,
where ρi ≡ λi 1 A i is a pressure and ∂V i ∂q ≡ AiJi is a volume gradient. Because the inequalities of the separation distance constraint still apply, we have equivalent complementarity conditions for the pressure and volume-based constraint formulation,

0 ≤ ρi ⊥ ∂Vi ∂q ( q0 + ∆ q) ≥ 0.
Approximation Now, instead of using many point-contact constraints, we use the total intersection volume to constrain the system and a single Lagrange multiplier ρ to compute the pressure acting on the contact patch. The intersection volume V(q) identifies the active constraints, its gradient can be seen as the sum of -AiJi of active constraints, and the forces on the DOFs due to ρ act in the direction of this volume gradient. We call this approximation the mono-volume contact constraint.

Note that the mono-volume approach computes a uniform pressure across the whole contact area, which is not always realistic. This limitation is addressed in Section 5, which proposes a multi-volume approach as the equivalence discussed above need not be reduced to only a single volume. Nevertheless, when solving our volumebased constraints, some important physical properties are guaranteed:

• the interpenetration volume cannot increase;

• the pressure can only act to push objects apart;

• the pressure acting on the contact surface is non-zero if and only if the interpenetration volume is not decreasing.

Stabilization It is also important to note that LDI is a discrete collision detection. In practice, we can only create our volumebased constraints if there is a non-zero intersection volume. If the intersection is zero, then there is simply no contact. While the intersection volume should theoretically always be zero, it is instead allowed to remain below some small amount of interpenetration volume V to permit continuous contact. Thus, we modify the complementarity conditions in our volume-based formulation. For the implicit velocity update, the MLCP of Equations 4 and 5 becomes

A∆ q = b + J T V ρ, (16) 0 ≤ JV ( q0 + ∆ q) ⊥ ρ ≥ 0, (17)
where JV = ∂V ∂q . For the post stabilization correction, the MLCP of Equations 6 and 7 becomes

A∆q = J T V ξ, (18)
-V ≤ -V(q0) + JV ∆q ⊥ ξ ≥ 0. (19)
The value of V is discussed in section 4.3. Notice that we are effectively using the equations given in Section 3 with g = -V, but with modifications to the position correction (as given in Equation 19) and the friction conditions (as described below in Section 4.2).

Friction

Friction forces counteract relative velocity in the tangent plane. In the standard contact model, a friction force is applied at the contact point. We need to modify this definition because our contact force is actually a pressure applied to a contact surface. The sum of the repulsion forces applied to the vertices of one object provides us with the approximate direction of the contact normal vector n. We then compute two orthogonal unit vectors t1 and t2 spanning the tangent plane using Gram-Schmidt orthogonalization. To define the equivalent of the point-contact relative velocity in our model, we need to compute a weighted difference of the velocities on the top and bottom contact surfaces C + and C -as defined by the normal direction (see Figure 4). Specifically, the weight of each vertex must take into account its associated contact surface area, which we approximate by the volume gradient with respect to the vertex position projected onto the volume-contact normal. Thus, we define

v rel ≡ k A k ṗk , with A k = n • ∂V ∂p k , ( 20
)
where ṗk is the velocity at vertex k, and A k is the area associated with vertex k projected to the tangent plane. The sign of A k depends on the orientation of the vertex's associated surface with respect to the contact normal n. This ensures that v rel will be zero when the two intersecting objects have the same velocity. Note that v rel has units of velocity times area, so the product of v rel with our contact pressure (force per unit area) gives power (velocity times force), like in the standard contact model.

Viscous contact forces are straightforward, and can applied proportionally to the velocity:

f k = -ν v rel A k .
To enforce sticking in Coulomb friction, we must constrain the projections of the relative velocity to the basis vectors of the tangent plane t1 and t2. At the velocity level, a no-slip constraint can be written as

A1t1 T • • • Ant1 T A1t2 T • • • Ant2 T    ∆ ṗ1 . . . ∆ ṗn    = -t1 T v rel -t2 T v rel ,
where the ∆ ṗk are the velocity corrections necessary to cancel the current sliding velocity. For Coulomb friction instead of using this no-slip condition, we can use Equations 8, 9, and 10, where the matrix T is now constructed by weighting the symmetric basis formed by the contact frame tangent vectors t1 and t2,

T =     A1t1 T • • • Ant1 T A1t2 T • • • Ant2 T -A1t1 T • • • -Ant1 T -A1t2 T • • • -Ant2 T     ∂p ∂q . ( 21
)
This will let us compute the appropriate relative velocity for our volume contact, as well as providing a basis for applying equal and opposite friction forces on the DOFs.

Our net friction forces are centered on the barycenter of each collision surface C + and C -. They are parallel with opposite directions but, in contrast with the repulsion forces, they are not aligned and this generates ghost torques. This is due to the non-physical situation where two objects intersect each other. We can compute the contact center, at the origin of the contact frame shown in Figure 4, as the barycenter of the contact surfaces:

c = k A k p k / k A k .
When volumetric velocity fields are available, forces applied to c can be propagated to the objects using the transpose of the velocity mapping. We can thus apply all the contact forces at this point rather than the surface vertices, obtain perfectly aligned opposite forces and avoid the ghost torque. We have not found this useful in practice, because the offset between the force directions (and thus the intensity of the ghost torque) is proportional to the intersection depth, which is small in our experiments.

Animation Loop

Let us now consider the case where multiple objects are in the system. We will use the vector V(q) to represent the interpenetration volumes between each pair of objects, and we set up the MLCP problems following the description in Sections 3.1 and 4.1. Our animation algorithm starts with the velocity MLCP of Equations 16 and 17, and additionally includes the frictional constraints. We then update the velocities and the positions. The velocities do not necessarily produce a position that will meet the V(q 0+h ) ≤ V constraint, thus we perform a post stabilization step. Due to the nonlinearity of the non-interpenetration constraint, we perform multiple iterations and monitor the intersection volume using Newton-Raphson's algorithm.

The post step position correction can ensure non-penetration between colliding objects; however, we must maintain a small amount of intersection rather than removing it completely. The volume derivatives used for the contact and friction formulations can lack precision if the interpenetration volume computed from the LDI involves only a few pixels. In practice, for non-breaking contacts we maintain an interpenetration layer that has a thickness of approximately half a pixel (except in the left of Figure 1, for a better view of the intersection volumes). But this measure of the intersection depth is not available in our contact model. Fortunately, we can estimate the surface of the contact by summing the positive A k of Equation 20. From this we compute the average intersection depth within a contact, and estimate the necessary volume variation to preserve the average depth within the desired range. Thus, we automatically regulate the value of V in the position MLCP of Equations 18 and 19 and we perform LDI collision detection to compute the new intersection volume and its Jacobian. This loop is repeated until the desired precision is obtained or a maximum number of iterations is reached. Each position correction deals only with contact volume correction because the friction constraints do not apply.

Algorithm 1 Animation Loop Input: current state q0, q0 Output: next state q 0+h , q0+h 1: ∆ q ← MLCPsolve velocity update with friction 2: q0+h ← q0 + ∆ q 3: q 0+h ← q0 + h q0+h 4: Compute LDI: intersection V(q 0+h ) and Jacobian 5: while V(q 0+h ) > V do 6:

∆q ← MLCPsolve position stabilization 7:

q 0+h ← q 0+h + ∆q 8:
Compute LDI: intersectionV(q 0+h ) and Jacobian 9: end while We implement the velocity update in two steps. The first step, prediction, performs an unconstrained implicit integration of the body forces. The second step, correction, modifies the velocities in order to meet the contact constraints. The position stabilization is implemented as a correction of positions. To solve the ML-CPs we use a Gauss-Seidel like method (see, e.g., [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF]), but a direct solver like Lemke's algorithm could also be used. We have good results with the integration scheme presented here (linearized backward Euler and post stabilization) but other time-stepping schemes could be used. The novelty lies in the use of volume constraints that are computed using LDI collision detection.

Multi-Volume Geometric Model

In the previous section, we have shown that the minimization of the intersection volume generates a single scalar equation for the repulsion between two objects. This is independent of their geometric resolution, which results in dramatic simplifications with respect to traditional distance-based methods. A limitation of this approach is that the contact equation holds for the whole contact area, resulting in a single behavior. For Coulomb friction, the whole contact is either sticking, or sliding, or the patch is breaking contact as a whole. We often expect more detailed behaviors. An example of this is shown in Figure 7, where a part of the snake is sliding while the other is sticking. The solution is to split the intersection volume in several parts with independent contact equations. We call this the multi-volume model, in contrast to the mono-volume model presented in the previous section.

We divide the intersection volume using a regular grid aligned with the LDI directions, as shown in the second image of Figure 1. This makes it easy to tune the spatial resolution of the contact. Figure 5 represents a cell of the grid. To set up separate contact equations, we need to compute the volume gradient of each cell. In the viewing direction, the boundary of the intersection volume in a cell may be The contributions of the object surface pixels to the gradient are given in Section 3.2. The remaining question is how to take into account the pixels of the cell walls. If the cell walls are considered fixed, then their pixels do not contribute to the volume gradient because their depth is independent of the object vertices. For example, the cell highlighted in Figure 5 would only apply a force to the red object, because it does not include the surface of the blue object. Computing the pressures independently in each cell would thus violate Newton's law of equal and opposite forces between the interacting objects. Physical soundness requires the pressure force applied to the cell walls to be dispatched on the object surfaces.

W + S + W - z W - ij z - ij z + ij 1-W + ij 1-® ®
For each pixel of the cell wall, we compute its barycentric coefficient inside the intersection volume in the viewing direction. We do this such that the depth of a pixel of a top cell wall W + is given by

W + ij = αijz + ij + (1 -αij)z - ij , (22)
while the depth of a pixel of a bottom cell wall W -is given by

W - ij = βijz + ij + (1 -βij)z - ij . (23)
An example of these barycentric coefficients is shown in Figure 5.

In each cell, the gradient of the intersection volume is thus

∂V ∂p z k = a (i,j)∈S + z ∂z + ij ∂p z k -a (i,j)∈S - z ∂z - ij ∂p z k +a (i,j)∈W + z αij ∂z + ij ∂p z k + (1 -αij) ∂z - ij ∂p z k -a (i,j)∈W - z βij ∂z + ij ∂p z k + (1 -βij) ∂z - ij ∂p z k (24)
Once the gradient is computed, the contact equations of the cell are straightforward and can be set up the same way as presented in Section 4.

Varying friction behaviors across the contact can greatly enhance the simulation of large deformable objects, as illustrated in the accompanying video. Moreover, separation can also be improved, as illustrated in Figure 6. The actual volume variation between convex objects is smaller than anticipated using the gradient, due to the fact that the area of the contact surface becomes smaller as the objects move apart, as can be seen in Figure 6. As such, some intersection will remain after the linear solve. A complete interpenetration position correction with the mono-volume approach requires several iterations to converge. Using the multi-volume approach, the most exterior contacts are progressively discarded during the MLCP solution, and the final contact occurs only where the intersection was initially the deepest.

We validated the multi-volume model on the example of a snake pulled laterally illustrated in Figure 7. In the mono-volume model, the snake spins around the contact center located approximately in the middle of its body, while it gently evolves to a nice stable bent shape using the multi-volume model. Additionally, when an object is cut as shown in Figure 2c, the multi-volume model allows the independent processing of each of the parts.

Results

We have successfully applied our method to a variety of scenes including complex rigid and deformable objects. The animation algorithm presented in Section 4.3 is used in all cases. For deformable objects, we use a co-rotational Finite Element Model [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF]]. We implemented our method within the open-source SOFA library [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF]]. When possible, we compare it with our implementations of distance-based contact models. We performed the simulations on a Intel Core i7 975 CPU with a Nvidia GeForce GTX 285 GPU. Table 1 presents the complexity and measured performance for several example simulations. The Collision Detection and Contact Modeling columns represent the time taken to compute the LDI and evaluate the intersection volumes and gradients (steps 4 and 8 in Algorithm 1). The Mechanical Motion column represents the unconstrained predictive motion due to internal mechanical forces, while the Constraint Solver column represents the velocity and position corrections.

The two rigid armadillos in Figure 8 are made of 345 944 triangles each and simulated at 41 fps (not including rendering), with 7 contacts on average, and at most 16 contacts. We also experimented with precomputed distance fields on this scene, since they are a popular and efficient tool for collision detection and response between rigid bodies. The distance fields detect and model collisions faster than our method, allowing refresh rates as high as 57 fps. However, they generate considerably more equations (typically 200 contacts in this example, hence 600 equations including friction), which results in longer solve times. This is especially true when a high resolution distance field is required to handle small geometric details. Our method does not suffer from this problem of mesh resolution sensitivity.

The advantages of our method are clear when distance fields are not available, for instance, when simulating deformable objects or topology changes (Figure 2c). A interesting feature of LDI-based methods is that self-intersections are no more costly to compute than inter-objects collisions. This is illustrated in Figure 2a, where 3 octopuses are simulated at interactive rates, while touching along multiple and complex contact and self-contact areas. The hand grasping the deformable object in Figure 2b is simulated at 54 fps, with 10 contacts, while the proximity-based method runs only at 23 fps due to the complexity of the collision detection and the larger number of contacts, up to 70 in this example.

The complex medical scene in Figure 1 is based on real data. Due to unavoidable data noise and approximations in the reconstruction process, the organs intersect each other at initialization time, and the distance-based method fails to repell them, making the simulation intractable. In contrast, our volume-based method is robust to deep intersections and has no problems separating the organs on the first time step. No spurious velocity is introduced for this correction since the position correction step does not modify the velocities. This simulation is difficult because the grasping relies on high pressure and friction forces. As expected, the tool eventually slips and releases the organ when pulled too far. This is an excellent example of the stability of static friction in our method; the contact forces stay within the friction cone during the first part of the movement, and then lie on the edge of the cone as the organ slips from the grasp of the tool. Note that in this case the geometry of the tool is thin with very fine teeth. It can be handled by our method but requires very small pixels, which is very costly as the LDI is currently rasterized with an uniform precision over the full scene.

While the method could be extended to locally adapt the resolution, which is an interesting avenue for future work, we instead replaced the collision mesh for this tool with a slightly simpler and thicker version. The scene includes 550 000 triangles, 9 500 degrees of freedom, and runs at 7 fps including rendering.

An obvious limitation of our approach is that interacting objects must have a volume, so very thin objects such as cloth are not supported. Likewise, since our LDI computations provide intersections at discrete time steps, it is also possible for small fast moving objects to tunnel through one another. For both of these cases, methods based on continuous collision detection are more suitable. Note that we need to maintain a small intersection to model the velocity constraints necessary for resting contacts with dry friction. For each contact volume, we are able to do this by controlling the average intersection depth. This average depth could be much smaller than the true maximum interpenetration depth, but we have not observed any problems related to this in practice.

Conclusion

We have presented a new formulation for contact using a hard constraint on the intersection volume. As such, the number of equations is small and independent of geometric complexity. We have extended the formulation to include Coulomb friction, and introduced a multi-volume grid model to better accommodate complex behavior across multiple contact patches. Our method is robust to deep interpenetrations and allows us to simulate frictional contact between complex deformable bodies at unprecedented rates.

Our results show examples that demonstrate that our method performs well for simulations involving sliding, rolling, and resting contact. In general, our multi-volume constraints can be used for a variety of difficult simulation scenarios involving both rigid and deformable objects in frictional contact, such as surgery simulation, stacking, grasping and cutting.

( a )Figure 2 :

 a2 Figure 2: Demonstration of several difficult scenarios for contact modeling that are efficiently handled by the proposed method.

Figure 3 :

 3 Figure 3: A 2D slice of an LDI showing the intersecting object volumes of two objects. Here, the LDI viewing direction is the z axis. Vertices are labeled p, while pixels are shown with horizontal lines in different columns. The intersection volume appears in purple, and is bounded by sets of surface pixels, S + on top, and S -on bottom (dotted lines).

Figure 4 :

 4 Figure 4: Two intersecting objects and the contact model. The gradient of the intersection volume appears as red lines at each vertex. The contact frame used for friction is built based on the estimated normal direction n.

Figure 5 :

 5 Figure 5: A slice of the intersection volume in a grid cell.

( a )

 a Initial State (b) Mono-Volume (c) Coarse Grid (d) Fine Grid

Figure 6 :

 6 Figure 6: Intersection after one LCP solution of volume constraint stabilization for grids of different resolution. Increasing the multivolume precision provides a more accurate linearization of the intersection volume gradients.

Figure 7 :

 7 Figure7: Multi-volume friction. Red, green, and blue arrows show the contact frames, while the yellow arrows denote forces. The table shows the measured difference (as fractions of the length of the snake) between the classical penetration-based approach and multi-volume constraints with varying cell sizes.

  Detailed surfaces (b) Distance fields (c) LDI rasterization (d) Contact volumes on the right arm of the green Armadillo

Figure 8 :

 8 Figure 8: Comparison with distance fields. Our method is able to handle highly detailed meshes, as visible in (d), without any precomputation.

Figure 9 :

 9 Figure 9: Stacking of 40 soft toys with complex geometry such as thin ears, sharp features, and a mixture of small and large triangles.

Table 1 :

 1 Complexity and measured performance for example simulations. The Objects, Triangles, LDI and Contacts columns respectively present the number of simulated objects, triangles, and upper limits on the LDI dimensions and generated contacts. The rest of the table presents the average time spent in the four main computation steps, as well as the overall speed of the simulation (excluding rendering).

	Simulation	Objects Triangles	LDI Pixels	Layers	Contacts	Collision Detection Modeling Contact Mechanical Constraint Motion Solver	Iterations Per Second
	Medical (Figure 1)	9	550K	224×168	32	16	22.69 ms	1.35 ms	17.85 ms	32.45 ms	13.3
	Octopuses (Figure 2a)	3	46K	392×352	28	110	5.67 ms	1.20 ms	43.43 ms	12.34 ms	15.6
	Grasping (Figure 2b)	2	9K	112×112	14	10	2.73 ms	0.97 ms	7.36 ms	6.18 ms	54.0
							1	3.45 ms	1.01 ms	2.72 ms	1.53 ms	109
	Snake (Figure 7)	1	7K	560×336	10	6 18	3.41 ms 3.45 ms	1.00 ms 0.95 ms	2.87 ms 2.98 ms	1.86 ms 3.09 ms	104 91.9
							41	3.46 ms	0.84 ms	2.99 ms	7.62 ms	65.3
	Armadillo	Rigids	2	128×107×98 distance field 692K 280×240 22	255 16	10.46 ms 18.48 ms	0.16 ms 1.15 ms	0.38 ms 0.48 ms	6.68 ms 4.08 ms	56.7 41.1
	(Figure 8)	FEM	2	692K	208×192	26	28	18.47 ms	1.53 ms	13.54 ms	17.98 ms	19.3
			10	205K	456×456	24	99	16.61 ms	1.64 ms	71.24 ms	34.82 ms	7.29
	Soft Toys (Figure 9)	20	379K	1232×472	34	242	28.79 ms	2.37 ms	146.94 ms 207.02 ms	2.44
			40	725K	1232×528	46	544	44.45 ms	3.43 ms	309.48 ms 723.83 ms	0.88


Julien Forget
email: julien.forget@onera.fr

Frédéric Boniol
email: frederic.boniol@onera.fr

David Lesens

Claire Pagetti
email: claire.pagetti@onera.fr

A Real-Time Architecture Design Language for Multi-Rate Embedded Control Systems *

Keywords: C.3 [Special purpose and application-based systems]: Real-time and embedded systems, D.2.11 [Software Architectures]: Languages, D.3.2 [Language classifications]: Data-flow languages Real-time systems, Synchronous, Scheduling, Compilation

This paper presents a language dedicated to the description of the software architecture of complex embedded control systems. The language relies on the synchronous approach but extends it to support efficiently systems with multiple real-time constraints, such as deadline constraints or periodicity constraints. It provides a high-level of abstraction and benefits from the formal properties of synchronous languages. The language defines a small set of rate transition operators, which enable the description of user-defined deterministic multi-rate communication patterns between components of different rates. The compiler of the language automatically translates a program into a set of communicating real-time tasks implemented as concurrent C threads that can be executed on a standard real-time operating system.

INTRODUCTION

Our work is specifically targeted for embedded control systems. An embedded control system consists of a control loop including sensors, control algorithms and actuators that regulate the state of the system in its environment. Spacecraft and aircraft flight control systems are good examples of control systems. Their objective is to control the position, speed and attitude of the vehicle thanks to physical devices, such as control surfaces for an aircraft or thrusters for a spacecraft.

These systems are highly critical as a malfunction can lead to dramatic consequences. Thus, their implementation must be deterministic, not only functionally (producing the right outputs for given inputs) but also temporally (producing outputs at the right time). They are also often multiperiodic, as different devices of the system have different physical characteristics and must therefore be controlled at different rates. Finally, the system must respect deadline constraints corresponding to maximum end-to-end latencies between an observation (inputs) and the corresponding reaction (outputs).

Case study: The Flight Application Software

To motivate our work, we consider the programming of an adapted version of the Flight Application Software (FAS) of the Automated Transfer Vehicle (ATV) designed by EADS Astrium Space Transportation for resupplying the International Space Station (ISS). The FAS handles all the software functionalities of the system as long as no fault is detected. Its architecture is described informally in Fig. 1. Each operation is represented by a box and arrows between boxes represent data-dependencies. Arrows without sources represent system inputs and arrows without destination represent system outputs. Data-dependencies define a partial order between the different operations of the system, as the consumer of a data-flow must execute after the producer of the data-flow.

The FAS first acquires data: the orientation and speed from gyroscopic sensors (GyroAcq), the position from the GPS (GPS Acq) and from the star tracker (Str Acq) and telecommands from the ground station (TM/TC). The Guidance Navigation and Control function (divided into an upstream part, GNC_US, and a downstream part, GNC_DS) then computes the commands to apply while the Failure Detection Isolation and Recovery function (FDIR) verifies the state of the FAS and looks for possible failures. Finally, commands are computed and sent to the control devices: thruster orders for the Propulsion Drive Electronics (PDE), power distribution orders for the Power System (PWS), solar panel positioning orders for the Solar Generation System (SGS) and telemetry towards the ground station (TM/TC), etc. Each operation has its own rate, ranging from 0.1Hz to 10Hz. An intermediate deadline constraint is imposed on data produced by the GNC Upstream (300ms while the period is 1s). Operations of different rates can communicate, which is an essential feature of such systems and has very important impacts on the complexity of the design and implementation process.

Contribution

Control systems like the FAS are very complex, their development involves several teams, which separately define the different functions of the system. The functions are then assembled by the integrator, who specifies the real-time properties of the functions and the communication schemes between the functions.

We propose a Real-Time Architecture Design Language designed to support this integration process. The language transposes the synchronous semantics [START_REF] Benveniste | The synchronous approach to reactive and real-time systems[END_REF] to the architecture design level. Defining the language as an ADL provides a high-level of abstraction, while relying on a synchronous semantics benefits from well-established formal properties. A program of this language consists of a set of imported nodes, implemented outside the program using other existing languages (C or Lustre for instance), and data-flows between the imported nodes. The language allows to specify multiple deadline and periodicity constraints on flows and nodes. It also defines a small set of rate transition operators that allows the programmer to precisely define the communication patterns between nodes of different rates. A program is automatically translated into a set of concurrent tasks implemented in C, which can be executed by standard real-time operating systems.

Paper Outline

This paper provides a complete presentation of the language and focuses on the illustration of its expressiveness through a series of examples. The formal semantics of the language is not detailed here and the compilation is only briefly described, due to lack of space. We first compare the language with existing work in Sect. 2. Then, we define the formal foundations of the language in Sect. 3. The main features of the language are detailed in Sect. 4. The expressiveness of the language is illustrated in Sect. 5. Finally, Sect. 6 gives an overview of the compilation process.

RELATED WORK

Probably the most widely used technique for program-ming real-time systems is to directly program the system with traditional imperative languages. These languages were not originally designed for real-time systems and the realtime aspects are mostly supported by Application Programming Interfaces (API) for real-time, which are very close to the underlying Real-Time Operating System (RTOS). For instance, the POSIX extensions for real-time [START_REF][END_REF] are a popular real-time API for C or Ada. Most real-time APIs share the same general principles: each operation of the system is first programmed as a separate function and functions are then grouped into threads, scheduled concurrently by the OS. As threads are related by data-dependencies, when dealing with critical systems the programmer must control the order in which data is produced and consumed (otherwise the behaviour of the system might not be predictable). This requires to implement synchronization mechanisms between threads, which can be a tedious and error-prone process.

This gradually leads programmers to consider languages with a higher-level of abstraction, which rely on an automated code generation process that translates the input program into lower level code. This greatly reduces the duration of the development process and produces more reliable software, as the correctness of the lower-level program is ensured by the translation process. The synchronous approach [START_REF] Benveniste | The synchronous approach to reactive and real-time systems[END_REF] proposes such a high level of abstraction, based on mathematical foundations, which make it possible to handle the compilation and the verification of a program in a formal way. This approach was originally mostly targeted for mono-periodic systems, ie systems where inputs and outputs are all acquired or produced at the same rate, but was later extended with real-time primitives, which enables the support of multi-periodic systems (different rates for different inputs/outputs). [START_REF] Bertin | TAXYS=Esterel+Kronos. a tool for verifying real-time properties of embedded systems[END_REF] supports the definition of deadline constraints in Esterel, but generates fully sequential code, which leads to a poor processor utilization factor in the case of multi-periodic systems. [START_REF] Curic | Implementing Lustre Programs on Distributed Platforms with Real-Time Constraints[END_REF] defines primitives to specify multiple deadline and periodicity constraints on a Lustre program, but the compilation process is very specific to the target architecture (TTA). [START_REF] Alras | Model-based design of embedded control systems by means of a synchronous intermediate model[END_REF] extends this work and proposes a more generic compilation process, which generates a set of real-time tasks executable on a standard multi-task RTOS. Task communications are handled using the protocol proposed by [START_REF] Sofronis | A memory-optimal buffering protocol for preservation of synchronous semantics under preemptive scheduling[END_REF], which ensures the conservation of the synchronous semantics of communicating tasks thanks to databuffering mechanisms based on task priorities. Our language has a similar approach but focuses on the description of communication patterns between operations of different rates.

Architecture Description Languages (ADL) provide an additional layer of abstraction compared to synchronous languages, by focusing on the modelling of the interactions between high-level components of a system (typically, by assembling threads instead of functions). Threads specified in an ADL specification are considered as "black-boxes", implemented outside the ADL specification. AADL [START_REF] Feiler | The architecture analysis & design language (AADL): An introduction[END_REF], CLARA [START_REF] Faucou | An ADL centric approach for the formal design of real-time systems[END_REF], or GIOTTO [START_REF] Henzinger | Giotto: A time-triggered language for embedded programming[END_REF] are specifically dedicated to the modelling of embedded real-time systems and describe a system as a set of communicating high-level real-time components. As far as we know, automated code generation is not supported by CLARA yet. Concerning AADL, the programmer only has a limited choice for the description of the communication schemes between components and cannot define its own schemes. This choice is even more limited in GIOTTO, where tasks can only consume data produced before their release dates (instead of before the actual beginning of their execution), effectively leading to high end-toend latencies. Finally, the MARTE profile for UML includes a powerful Clock Constraint Specification Language (CCSL) [START_REF] André | Combining CCSL and Esterel to specify and verify time requirements[END_REF] for specifying the timing properties of a system. CCSL is far more expressive than our language, however it is targeted for system specification and system verification, not for automated code generation. Due to its expressivness, efficient code generation would most likely require to only consider a subset of the clock constraints provided by the language.

A possible approach to our problem would be to consider multi-periodic systems as a special case of Globally Asynchronous Locally Synchronous systems to reuse existing work [START_REF] Berry | Multiclock esterel[END_REF][START_REF] Le Guernic | Polychrony for system design[END_REF][START_REF] Carloni | Theory of Latency-Insensitive Design[END_REF][START_REF] Caspi | About the design of distributed control systems: The quasi-synchronous approach[END_REF][START_REF] Cohen | N -Synchronous Kahn Networks: a relaxed model of synchrony for real-time systems[END_REF]. However, multi-periodic systems are in essence completely synchronous, treating them as such allows to more clearly identify the real-time properties of a system (during the specification), provides better analysis possibilities (proving the correct synchronization of the operations, performing schedulability analysis) and allows the production of more optimized code (specific communication protocols).

Our language can be compared with other languages as follows:

• ADLs: The level of abstraction of the language is the same as that of existing ADLs but the language relies on a different semantics (the synchronous approach) and allows the programmer to define its own communication patterns between tasks. Furthermore, the complete compilation scheme is defined formally;

• Synchronous languages: Though the language has a synchronous semantics and shares similarities with Lustre [START_REF] Halbwachs | The synchronous data-flow programming language LUSTRE[END_REF], it addresses a different phase of the development (the integration) and enables efficient translation into a set of real-time tasks;

• Imperative languages: Though a program is in the end translated into (correct-by-construction) imperative code, the program is written in a data-flow style, which is easier to analyze and to verify.

SYNCHRONOUS REAL-TIME

This section presents the formal foundations of the language, which rely on the synchronous data-flow approach. The synchronous model is extended to support multi-rate systems more efficiently.

Synchronous Data-Flow

The synchronous approach simplifies the programming of real-time systems by abstracting from real-time and reasoning instead on a logical-time scale, defined as a succession of instants, where each instant corresponds to a reaction of the system. The synchronous hypothesis says that the programmer does not need to consider the dates at which computations take place during an instant, as long as the computations finish before the beginning of the next instant.

The synchronous approach has been implemented in several programming languages, our language is a synchronous data-flow language and is thus close to Lustre [START_REF] Halbwachs | The synchronous data-flow programming language LUSTRE[END_REF], Signal [START_REF] Benveniste | Synchronous programming with events and relations: the Signal language and its semantics[END_REF] or Lucid Synchrone [START_REF] Pouzet | Lucid Synchrone, version 3. Tutorial and reference manual[END_REF]. Exactly as in Lustre, the variables and expressions of a program are flows, that is to say infinite sequences of values. The clock of a flow defines the instants during which the values of the flow must be computed. A program consists of a set of equations, structured into nodes. The equations of a node define its output flows from its input flows. A very simple program is given below: imported node p l u s 1 ( i : i n t ) returns ( o : i n t ) wcet 5 ; node N( i : i n t ) returns ( o : i n t ) l e t o=p l u s 1 ( i ) ; t e l

The program first declares an imported node (plus1), implemented outside the program using other existing languages (C or Lustre for instance). The imported node has one integer input and one integer output. The duration of the node is also declared (this will be explained later). The program then defines a node N, with one input i and one output o and one equation. At each reaction of the program, the current value of the flow o is defined as the result of the application of plus1 to the current value of the flow i. Imported node calls follow the usual data-flow semantics: an imported node cannot start its execution before all its inputs are available and produces all its outputs simultaneously at the end of its execution.

Relating Instants and Real-Time

In a classic synchronous data-flow program, the program describes the behaviour of the system for each of its basic iterations. Such an approach is not very well suited to the case of multi-periodic systems. Indeed, if we consider the system of Fig. 2, which consists of a fast operation F with a rate of 10Hz and a slow operation S with a rate of 2Hz, it is unclear what a basic iteration should consist of: one iteration of F and a varying part of S or one iteration of S and five iterations of F ? time

F F F F F F F F S S
basic iteration ? (10Hz) basic iteration ? (2Hz)

Figure 2: Base rate for a multi-periodic system Therefore, we define a synchronous model that allows to relate instants to real-time. The main idea of this model is that a multi-periodic system can be considered as a set of locally mono-periodic synchronous processes assembled together to form a globally multi-periodic synchronous system. Locally, each process has its own logical time scale. When we assemble processes of different rates, we assemble processes of different logical time scales so we need a common reference to compare instants belonging to different logical time scales. This common reference is the real-time scale.

Our synchronous real-time model relies on the Tagged-Signal Model [START_REF] Lee | Comparing models of computation[END_REF]. Given a set of values V, a flow is a sequence of pairs (vi, ti) i∈N where vi is a value in V and ti is a tag in Q, such that for all i, ti < ti+1. The clock of a flow, is its projection on Q. A tag represents an amount of time elapsed since the beginning of the execution of the program. Following the relaxed synchronous hypothesis of [START_REF] Curic | Implementing Lustre Programs on Distributed Platforms with Real-Time Constraints[END_REF], at each activation a flow is required to be computed before its next activation. Thus, each flow has its own notion of instant and the duration of the instant ti is ti+1 -ti. Two flows are synchronous if the durations of their instants are the same. We can compare the durations of the instants of two flows to determine if a flow is faster than the other and more importantly, we can determine how much faster this flow is.

Strictly Periodic Clocks

A clock is a sequence of tags, we define a particular class of clocks called strictly periodic clocks as follows:

Definition 1. ( Strictly periodic clock). A clock h = (ti) i∈N , ti ∈ Q, is strictly periodic if and only if: ∃n ∈ Q + * , ∀i ∈ N, ti+1 -ti = n
n is the period of h, denoted π(h) and t0 is the phase of h, denoted ϕ(h).

A strictly periodic clock defines the real-time rate of a flow and is uniquely characterized by its phase and by its period:

Definition 2. The term (n, p) ∈ Q + * × Q denotes the strictly periodic clock α such that: π(α) = n, ϕ(α) = π(α) * p
We then define clock transformations specific to strictly periodic clocks, which produce new strictly periodic clocks: Definition 3. Let α be a strictly periodic clock, operations /., * . and →. are periodic clock transformations, that produce new strictly periodic clocks satisfying the following properties:

π(α/.k) = k * π(α), ϕ(α/.k) = ϕ(α), k ∈ N * π(α * . k) = π(α)/k, ϕ(α * . k) = ϕ(α), k ∈ N * π(α →. q) = π(α), ϕ(α →. q) = ϕ(α) + q * π(α), q ∈ Q
Divisions and multiplications respectively decrease or increase the frequency of a clock while phase offsets shift the phase of a clock.

Strictly periodic clocks are actually a subset of the usual clocks of synchronous languages, which are defined using Boolean activation conditions. However, by restricting to this particular class of clocks, we are able to specify the realtime properties of a system more explicitly and to compile a program efficiently into a set of real-time tasks.

LANGUAGE DEFINITION

This section provides a detailed presentation of the language, focusing on its real-time primitives.

Real-time constraints

Real-time constraints are specified on node inputs and outputs, thus they represent constraints related to the environment of the node instead of constraints related to implementation concerns. First, the programmer can specify the period and the phase of an input/output as follows:

node p e r i o d i c ( i : i n t r a t e ( 1 0 , 1 / 2 ) ) returns ( o : i n t r a t e ( 1 0 , 0 ) ) l e t . . . t e l i has period 10 and phase 1 2 , so its clock sequence is (5, 15, 25, 35, . . .) Time units correspond to whichever metric the underlying operating system uses for durations.

The programmer can also specify deadline constraints:

node d e a d l i n e ( i : i n t r a t e ( 1 0 , 0 ) b e f o r e 2 ) returns ( o : i n t r a t e ( 1 0 , 0 ) due 8 ) l e t . . . t e l

x: due d specifies that x has a deadline constraint of d. For each tag ti in the clock of x, the computation of the value of x for this tag must complete before date ti +d. x: before d specifies that x should be acquired by the program before the deadline d. In practice, this means that, before this deadline, the program must copy the value produced by the sensor to a local variable accessible to the other operations of the program. Finally, the programmer must specify the duration of each imported node of the program. For instance, the following declaration specifies that the worst case execution time of plus1 is 5 time units.

imported node p l u s 1 ( i : i n t ) returns ( o : i n t ) wcet 5 ;

Rate transition operators

We define a series of rate transition operators, based on periodic clock transformations. These operators enable the definition of user-specified communication patterns between nodes of different rates. First, the programmer can oversample or under-sample a flow periodically as follows:

node s a m p l i n g ( i : r a t e ( 1 0 , 0 ) ) returns node o f f s e t s ( i : r a t e ( 1 0 , 0 ) ) returns ( o : r a t e ( 1 0 , 1 ) ; p : r a t e ( 1 0 , 0 ) ; q : r a t e ( 1 0 , 1 / 2 ) ) l e t o=t a i l ( i ) ; p = 0 : : o ; q=i ˜>1/2; t e l tail (e) drops the first value of e. If e has clock α, then tail (e) has clock α →. 1, so tail (e) becomes active one period later than e. cst :: e produces cst one period earlier than the first value of e and then produces e. If e has clock α, then cst :: e has clock α →. -1. If e has clock α, e ∼> q delays each value of e by q * π(α) (with q ∈ Q + ). The clock of e ∼> q is α →. q. The behaviour of these operators is illustrated in Fig. 4

Polymorphic clocks

As we have seen in some of the previous examples, the types of the inputs and outputs of a node can be left unspecified. In this case, they will be inferred by the compiler. Similarly, clocks can be left unspecified, in which case they will be inferred by the compiler too. Furthermore, rate transition and phase offset operators can be applied to flows the rate of which is unspecified. For instance, the following node definition is perfectly valid: node u n d e r s a m p l e ( i ) returns ( o ) l e t o=i / ˆ2 ; t e l This program simply specifies that the rate of o is half that of i (ie twice the period of o), regardless of what the rate of i may be. The actual rate of the flows will only be computed when instantiating the node. The node can even be instantiated with different rates in the same program, for instance: node p o l y ( i : i n t r a t e ( 1 0 , 0 ) ; j : i n t r a t e ( 5 , 0 ) ) returns ( o , p : i n t ) l e t o=u n d e r s a m p l e ( i ) ; p=u n d e r s a m p l e ( j ) ; t e l

The clock inferred for o is (20, 0), while the clock inferred for p is (10, 0). Such clock polymorphism increases the modularity with which systems can be programmed.

PROGRAMMING MULTI-RATE SYSTEMS

A key feature of the language is that it gives the programmer the freedom to choose the pattern used for the communications between operations of different rates, while most languages only allow a small set of predefined patterns. Rate transition operators apply very simple flow transformations but can be combined to produce various communication patterns. The formal definition of the language and the static analyses performed by the compiler ensure that only deterministic patterns are used.

We detail a series of communication patterns and show that common patterns can easily be defined as generic library operators that can be reused in different programs. We illustrate the different patterns on the simple example of Fig. 5, where a fast operation F (of period 10ms), exchanges data with a slow operation S (of period 30ms) and operation F exchanges data with the environment of the program.

Sampling

The simplest communication pattern that can be implemented with our language is based on data sampling. For the example of Fig. 5, we obtain the node sampling presented previously in Sect. 4.2. Notice that we need a delay either from S to F or from F to S, otherwise there is a causality loop between F and S, the compiler cannot find an order for the execution of the program that respects all the datadependencies. In the sampling version of the system, the delay avoids the reduction of the deadline of operation S.

Let us now consider an alternative implementation:

node s a m p l i n g 2 ( i : r a t e ( 1 0 , 0 ) ) returns ( o )

var v f , v s ; l e t ( o , v f )=F( i , v s * ˆ3 ) ; v s=S ( ( 0 fby v f ) / ˆ3 ) ; t e l
This version is also correct but an important difference with the previous version is that S must end before a relative deadline of 10ms, as it must end early enough for F to end before its own deadline (10ms). In the previous version, the delay avoids this restrictive deadline. The behaviours of these two examples are compared in Fig. 6 Using the operator tail, it is possible to choose which of the successive values of the fast operation should be consumed (ie not necessarily the first). For instance, we can consume the second value instead of the first: node s a m p l i n g t a i l ( i : r a t e ( 1 0 , 0 ) ) returns ( o ) var v f , v s ; l e t ( o , v f )=F( i , ( 0 : :

( ( 0 fby v s ) * ˆ3 ) ) ) ; v s=S ( ( t a i l ( v f ) ) / ˆ3 ) ; t e l
The behaviour of this program is illustrated in Fig. 7. The previous values are initialised to the value init. This is illustrated in Fig. 8. Then we define two more complex rate transition nodes: The node count_n used in this example is a counter modulo n. The node split handles communications from slow to fast operations. It takes an array of size n as input and splits it into n successive values. The outputs are produced n times faster than the input. The node join handles communications from fast to slow operations. It joins n successive values of its input into an array. The output is n times slower than the input. This behaviour is illustrated below:

i i0 i1 i2 i3 ... A [0, 0, i0] [0, i0, i1] [i0, i1, i2] [i1, i2, i3] ...
i i0 i1 i2 i3 i4 i5 ... i'=join(3,i,0) [0, 0, i0] [i1, i2, i3] ... split(3,i') 0 0 i0 i1 i2 i3 ...
We can then define a different version of the example of Fig. 5 that uses a queuing communication pattern: node q u e u i n g ( i : r a t e ( 1 0 , 0 ) ) returns ( o )

var v f , v s ; l e t ( o , v f )=F( i , s p l i t ( 3 , 0 fby v s ) ) ; v s=S ( j o i n ( 3 , v f , 0 ) ) ; t e l
This version requires flows vf and vs to be arrays of integers instead of simple integers and thus node S must return an array, not just an integer.

Mean value

We give a third example of communication pattern, which is an intermediate between sampling and queuing. Instead of queuing the successive values produced by the fast operations, we take their mean value: node mean n ( c o n s t n : i n t ; i , i n i t ) returns (m: i n t )

var n v a l u e s , n sum : i n t ˆn ; l e t n v a l u e s=s t o r e n ( n , i , i n i t ) ;

n sum [ 0 ] = n v a l u e s [ 0 ] ; n sum [ 1 . . n-1]= n v a l u e s [ 1 . . n-1]+n sum [ 0 . . n -2 ] ; m=(n sum [ n -1]/ n ) / ˆn ; t e l
Notice that there is no causality loop in the third equation as each n_sum[i] depends on n_sum[i-1], which are independent variables. The behaviour of this node is illustrated in Fig. 9.

i 1 2 3 4 5 ... n_values [0, 0, 1] [0, 1, 2] [1, 2, 3] [2, 3, 4] [3, 4, 5] ... n_sum [0, 0, 1] [0, 1 , 3] [1, 3, 6] [2, 5, 9] [3, 7, 12]
... m 0 3 ...

A Complete Example

Fig. 10 gives a program for the FAS of Fig. 1 using communication patterns based on data sampling. Each operation of the system is first declared as an imported node. The program then specifies the wcet of each sensor and each actuator (for each input/output of the node FAS). The node FAS is the main node of the program. It first specifies the inputs (gyro, gps, str, tc) and outputs (pde, sgs, gnc, pws, tm) of the system with their rates (for instance gyro has period 100) and deadline constraints (for instance gnc has deadline 300). The body of the node then specifies the real-time architecture of the system. For instance, the first equation says that the node Gyro_Acq computes the variable gyro_acq from the input gyro and from the variable tm. As tm is produced with a period of 10s while Gyro_Acq has a period of 100ms, we over-sample tm by a factor 100. We use a delay before performing the over-sampling, to avoid reducing the deadline of node TM_TC. Some flows are over-sampled without using a delay, for instance when the outputs of the acquisition nodes are consumed by faster nodes (gps_acq*^10 in the equation of the FDIR for instance). This implies that the acquisition nodes will have a deadline much shorter than their period, which means that they execute less frequently but they must respond fast. Communications from fast to slow nodes are handled using the under-sampling operator /^. For instance, the node TM_TC has a period of 10s and consumes the flow fdir_tm, which is produced with a period of 100ms, so we under-sample fdir_tm by a factor 100. Finally, in the equation of node PWS we apply a phase offset of one half to its input gnc_pws, which causes the node PWS to execute with a phase of half its period (it has period 1000 and phase 500).

COMPILATION

This section gives an overview of the compilation of a program.

Static Analyses

The compilation puts strong emphasis on the verification of the correctness of the program to compile. This consists of a series of static analyses, which are performed before code generation. The first analysis is the type-checking. The language is a strongly typed language, in the sense that the execution of a program cannot produce a run-time type error. The type-checking is fairly standard and follows [START_REF] Pierce | Types and programming languages[END_REF].

The causality check verifies that the program does not contain cyclic definitions: a variable cannot instantaneously depend on itself (i.e. not without a fby in the dependencies). For instance, the equation x=plus1(x); is incorrect, it is similar to a deadlock since we need to evaluate plus1(x) to evaluate x and we need to evaluate x to evaluate plus1(x). This analysis is similar to that of Lustre [START_REF] Halbwachs | Generating efficient code from data-flow programs[END_REF].

The clock calculus verifies that a program only combines flows that have the same clock. When two flows have the same clock, they are synchronous as they are always present at the same instants. Combining non-synchronous flows leads to non-deterministic programs as we access to undefined values. For instance we can only compute the sum of two synchronous flow, the sum of two non-synchronous flows is ill-defined (ambiguous). The clock calculus ensures that a synchronous program will never produce ill-defined values. The clock calculus of the language is more specific than other analyses. It extends the work of [START_REF] Colaço | Clocks as first class abstract types[END_REF] and is defined imported node Gyro Acq ( gyro , t c : i n t ) returns ( o : i n t ) wcet 3 ; imported node GPS Acq ( gps , t c : i n t ) returns ( o : i n t ) wcet 3 ; imported node S t r A c q ( s t r , t c : i n t ) returns ( o : i n t ) wcet 3 ; imported node FDIR( gyr , gps , s t r , gnc : i n t ) returns ( t o p d e , t o g n c , t o t m : i n t ) wcet 1 5 ; imported node GNC US( f d i r , gyr , gps , s t r : i n t ) returns ( o : i n t ) wcet 2 1 0 ; imported node GNC DS( us : i n t ) returns ( pde , s g s , pws : i n t ) wcet 3 0 0 ; imported node TM TC( f r o m g r , f d i r : i n t ) returns ( cmd : i n t ) wcet 1 0 0 0 ; imported node PDE( f d i r , gnc : i n t ) returns ( p d e o r d e r : i n t ) wcet 3 ; imported node SGS( gnc : i n t ) returns ( s g s o r d e r : i n t ) wcet 3 ; imported node PWS( gnc : i n t ) returns ( p w s o r d e r : i n t ) wcet 3 ; s e n s o r g y r o wcet 1 ; s e n s o r gps wcet 1 ; s e n s o r s t r wcet 1 ; s e n s o r t c wcet 1 ; a c t u a t o r pde wcet 1 ; a c t u a t o r s g s wcet 1 ; a c t u a t o r gnc wcet 1 ; a c t u a t o r pws wcet 1 ; a c t u a t o r tm wcet 1 ; node FAS( g y r o : r a t e ( 1 0 0 , 0 ) ; gps : r a t e ( 1 0 0 0 , 0 ) ; s t r : r a t e ( 1 0 0 0 0 , 0 ) ; t c : r a t e ( 1 0 0 0 0 , 0 ) ) returns ( pde , s g s ; gnc : due 3 0 0 ; pws , tm ) var g y r o a c q , g p s a c q , s t r a c q , f d i r p d e , f d i r g n c , f d i r t m , gnc pde , g n c s g , gnc pws ; l e t g y r o a c q = Gyro Acq ( gyro , ( 0 fby tm ) * ˆ1 0 0 ) ; g p s a c q = GPS Acq ( gps , ( 0 fby tm ) * ˆ1 0 ) ; s t r a c q = S t r A c q ( s t r , 0 fby tm ) ; ( f d i r p d e , f d i r g n c , f d i r t m ) =FDIR( g y r o a c q , g p s a c q * ˆ1 0 , s t r a c q * ˆ1 0 0 , ( 0 fby gnc ) * ˆ1 0 ) ; gnc=GNC US( f d i r g n c / ˆ1 0 , g y r o a c q / ˆ1 0 , g p s a c q , s t r a c q * ˆ1 0 ) ; ( gnc pde , g n c s g s , gnc pws )=GNC DS( gnc ) ; pde = PDE( f d i r p d e , ( 0 fby g n c p d e ) * ˆ1 0 ) ; s g s = SGS( g n c s g s ) ; pws=PWS( gnc pws ˜>1/2); tm = TM TC( t c , f d i r t m / ˆ1 0 0 ) ; t e l Figure 10: The FAS, with data-sampling in details in [START_REF] Forget | A multi-periodic synchronous data-flow language[END_REF]. If all the static analyses succeed, then the program has a well-defined, deterministic semantics and can be translated into lower-level code.

Muti-Task Code Generation

The compiler translates the program into C code. However, the translation differs from the classic "single-loop" sequential code generation, as we deal with multi-periodic systems. Let us consider the example below: imported node A( i : i n t ) returns ( o : i n t ) wcet 1 ; imported node B( i : i n t ) returns ( o : i n t ) wcet 5 ; node m u l t i ( i : r a t e ( 3 , 0 ) ) returns ( o ) var v ; l e t v=A( i ) ; o=B( v / ˆ3 ) ; t e l Out of three successive repetitions of A, only the first one must execute before a repetition of B. As shown in Fig. 11, a sequential execution of the program cannot respect all the deadline constraints of the program, the second repetition of operation A will miss its deadline (at date 6). We cannot find a sequence of execution for operations A and B that will respect the periodicity constraints of the two operations. However, if we implement operations A and B as two separate tasks and execute them concurrently with a preemptive RTOS, we can find a schedule for the execution of the tasks that respects the periodicity constraints of the two operations. Thanks to preemption, it is possible to start the execution of the second repetition of A earlier and thus to find a schedule that respects the deadline constraints of both tasks, as shown in Fig. 12. Marks on the time axis of B represent task preemptions and restorations: at date 3, B is suspended to execute the second repetition of A. After A completes, B is restored (at date 4) and resumes its execution where it stopped.

Overview of the Code Generation

The compiler translates the program into a set of communicating real-time tasks. Each imported node call of the program is translated into a task. The real-time attributes (Ti, Ci, ri, di) of a task τi are obtained as follows: the period Ti is the period of the clock of the corresponding imported node, the initial release date ri is the phase of the clock, the execution time Ci is the wcet specified for the imported node and the relative deadline di is either that specified for an input/output or by default equal to the period. For instance, the node multi described above is translated into two tasks: τA of attributes (3, 1, 0, 3) and τB of attributes [START_REF] Chetto | Dynamic scheduling of real-time tasks under precedence constraints[END_REF][START_REF] Berry | Multiclock esterel[END_REF]0,[START_REF] Chetto | Dynamic scheduling of real-time tasks under precedence constraints[END_REF].

Tasks are related by data-dependencies, each task consumes data produced by other tasks (the inputs of the imported node) and produces data for other tasks (the outputs of the imported node). This implies precedence constraints between tasks: a task cannot start its execution before all the tasks that produce its inputs complete their execution and the task produces all its outputs simultaneously at its completion. A simple way to handle such inter-task communications is to allocate a buffer for each communication in a global memory, the producer of the data writes to this buffer when it executes and the consumer reads from this buffer when it executes. However, to preserve the semantics of the initial program, communications must respect two important properties. First, the producer must write before the consumer reads. Second, as we deal with periodic tasks, we must ensure that a new instance of the producer does not overwrite the value produced by its previous instance, if this previous value is still needed by some task.

To satisfy the first property, the compiler extends the technique proposed by [START_REF] Chetto | Dynamic scheduling of real-time tasks under precedence constraints[END_REF] to the case of multi-periodic tasks. This consists in encoding task precedences by adjusting task real-time attributes. The deadline of each task is adjusted so as to be lower than that of its successors, and the release date of each task is adjusted so as to be higher that of its predecessors. The encoded task set can then simply be scheduled by the Earliest-Deadline-First policy [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard-real-time environment[END_REF]. As priorities are affected based on deadlines, each task will have a priority higher than that of its successors and lower than that of its predecessors, thus the first requirement is satisfied.

To satisfy the second property, the compiler implements a communication protocol based on data-buffering mechanisms, which ensure that the inputs of a task remain available until its deadline. The protocol is similar to that of [START_REF] Sofronis | A memory-optimal buffering protocol for preservation of synchronous semantics under preemptive scheduling[END_REF] but is specifically targeted and optimized for multi-periodic systems. An important feature of this protocol is that it requires no specific synchronization primitives (semaphores).

CONCLUSION

We proposed a language for programming multi-rate embedded control systems. It transposes the synchronous semantics to an Architecture Design Language. It enables the description of a set of high-level communicating components with multiple deadline and periodicity constraints. The language focuses on the precise description of deterministic multi-rate communication patterns. The compiler translates a program into a set of communicating real-time tasks that respects the semantics of the original program. A prototype of the compiler has been implemented in OCaml (approximately 4000 lines of code). It generates a set of concurrent C threads and relies on the POSIX extensions for real-time to handle the real-time aspects.

Figure 1 :

 1 Figure 1: The Flight Application Software

  ( o : r a t e ( 1 0 , 0 ) ) var v f , v s ; l e t ( o , v f )=F( i , ( 0 fby v s ) * ˆ3 ) ; v s=S ( v f / ˆ3 ) ; t e l e/ ˆk only keeps the first value out of each k successive values of e. If e has clock α then e/ˆk has clock α/.k. on the opposite, e * ˆk duplicates each value of e, k times. If e has clock α, then e * ˆk has clock α * . k. The fby operator is borrowed from Lucid Synchrone. cst fby e first produces cst and then the values of e delayed by the period of e.The clock of cst fby e is the same as the clock of e. The behaviour of these operators is illustrated in Fig.6

Figure 3 :

 3 Figure 3: Periodic sampling operators Three different operators based on clock phase offsets are available, as illustrated in the following program:

Figure 4 :

 4 Figure 4: Phase offset operators

Figure 5 :

 5 Figure 5: Simple multi-periodic communication loop

Figure 6 :

 6 Figure 6: Data sampling: the place of the delay

Figure 7 :

 7 Figure 7: Data sampling: starting the sample on the second value5.2 QueuingCombining rate transition operators and Lustre arrays, we can implement different communication patterns based on data queuing. First, we define a node that stores n successive values of a flow:node s t o r e n ( c o n s t n : i n t ; i , i n i t ) returns (A: i n t ˆn ) l e t A[ n-1]= i ; A [ 0 . . ( n -2)]=( i n i t ˆ( n -1)) fby (A [ 1 . . ( n -1 ) ] ) ) ; t e lint^n denotes an integer array of size n. A[n] denotes the n th value of array A. init^(n-1) denotes an array of size n-1 where each value of the array is set to init. A[1..(n-1)] denotes a slice of array A, from cell 1 to cell n -1. The behaviour of this node is the following. The n successive values of flow i are stored in the array A of size n. At each iteration, the node queues the current value of i at the end of A, moving the previous values of i by one cell towards the beginning of the array. At each instant, the array A contains the current value of i along with its n -1 previous values. The previous values are initialised to the value init. This is illustrated in Fig.8.

Figure 8 :

 8 Figure 8: Storing successive values of a flow (n = 3, init = 0)

Figure 9 :

 9 Figure 9: Computing the mean value of 3 successive "fast" values, (with init = 0) We can then define a different version of the example of Fig. 5 that takes the mean value of 3 successive fast values (flow vf) and over-samples fast values (flow vs): node s a m p l i n g ( i : r a t e ( 1 0 , 0 ) ) returns ( o ) var v f , v s ; l e t ( o , v f )=F( i , ( 0 fby v s ) * ˆ3 ) ; v s=S ( mean n ( 3 , v f , 0 ) ; t e l

9 BFigure 11 :

 911 Figure 11: Sequential execution of node multi

Figure 12 :

 12 Figure 12: Preemptive scheduling of node multi

  .

	date	0	5	10	15	20	25	30	...
	i	i0		i1		i2		i3	...
	o=tail(i)			i1		i2		i3	...
	p=0::o	0		i1		i2		i3	...
	i∼>2		i0		i1		i2		...

  .

	date	0	10	20	30	40	50	60	...
	vf	vf0	vf1	vf2	vf3	vf4	vf5	vf6	...
	vs	vs0			vs1			vs2	...
	vf/^3	vf0			vf3			vf6	...
	(0 fby vs)*^3	0	0	0	vs0	vs0	vs0	vs1	...
	(0 fby vf)/^3	0			vf2			vf5	...
	vs*^3	vs0	vs0	vs0	vs1	vs1	vs1	vs2	...

* This work was funded by EADS Astrium Space Transportation