Andrea Del Prete

Nicolas Mansard
email: nmansard@laas.fr

Addressing Constraint Robustness to Torque Errors in Task-Space Inverse Dynamics

Task-Space Inverse Dynamics (TSID) is a wellknown optimization-based technique for the control of highlyredundant mechanical systems, such as humanoid robots. One of its main flaws is that it does not take into account any of the uncertainties affecting these systems: poor torque tracking, sensor noises, delays and model uncertainties. As a consequence, the resulting control trajectories may be feasible for the ideal system, but not for the real one. We propose to improve the robustness of TSID by modeling uncertainties in the joint torques as additive white random noise (similarly to LQG). This results in a stochastic optimization problem, in which we can maximize the probability to satisfy the inequality constraints (i.e. to be feasible). Since computing this probability is computationally expensive, we propose three ways to approximate it that are much faster to compute and that we can then use for online control (resolution time below 1 ms). Simulation results show that taking robustness into account greatly increases the chances to have feasible control trajectories (even when the uncertainties affecting the system are not the one modeled in the controller). Robot Model Classic TSID Robust TSID p ind p(τ ) Ineq. Violation p(τ ) Ineq. Violation Perfect 60.29 No 98.92 No Corrupted 1 53.

I. INTRODUCTION

Task-space inverse dynamics (TSID) has become an increasingly popular way to control humanoid and quadruped robots [START_REF] Koolen | Summary of Team IHMC s Virtual Robotics Challenge Entry[END_REF][START_REF] Faraji | Versatile and Robust 3D Walking with a Simulated Humanoid Robot (Atlas): a Model Predictive Control Approach[END_REF][START_REF] Herzog | Balancing experiments on a torque-controlled humanoid with hierarchical inverse dynamics[END_REF][START_REF] Nakanishi | Operational Space Control: A Theoretical and Empirical Comparison[END_REF][START_REF] Nori | iCub Whole-Body Control through Force Regulation on Rigid Non-Coplanar Contacts[END_REF]. This technique offers a number of attractive features: it is theoretically sound [START_REF] Del Prete | Prioritized Motion-Force Control of Constrained Fully-Actuated Robots[END_REF], it works well in simulation [START_REF] Sentis | Synthesis of wholebody behaviors through hierarchical control of behavioral primitives[END_REF][START_REF] Oscar E Ramos | Towards reactive vision-guided walking on rough terrain: an inversedynamics based approach[END_REF], and it is fast enough to be applied for online control. However, as usual, the gap between simulation and real world is large and can be explained through countless unmodeled uncertainties affecting these systems, such as poor torque control, model uncertainties, sensor noises and delays. This results in control trajectories that are feasible for the ideal system, but not for the real one.

To improve the robustness of the control trajectories, we propose to account for additive random noise on the joint torques. First, Section II introduces the problem arising from solving an optimization problem without accounting for uncertainties. Then Section II-A focuses on the case in which the problem variable is affected by additive random noise. Unfortunately, the resulting stochastic optimization problem 1 is too computationally expensive to be used for online control. In Sections II-B and II-C, we discuss three different ways to approximate a general stochastic optimization problem with linear inequality constraints. Section III shows how these ideas relate to the specific TSID control problem. We then 1 In a stochastic optimization problem some of the variables are random and follow a known probability distribution

Minimum-norm solution

Active-set solution

Inequality constraints

Set of solutions of Ax = a Bx + b 0 Fig. 1: 2D example of an inequality-constrained least-square problem, such as [START_REF] Koolen | Summary of Team IHMC s Virtual Robotics Challenge Entry[END_REF], solved by an active-set method.

discuss how to use the three presented approximations to get different formulations of robust TSID (Sections IV and V). Finally, in Section VII we validate the proposed methods on a simulated HRP-2 humanoid robot. We show that taking robustness into account greatly increases the chances to have feasible trajectories for a system subject to uncertainties (even when these uncertainties are not the one modeled in the controller). Moreover, we verify that we can solve the proposed optimization problems in less than 1 ms on a standard CPU, so that these formulations are suitable for online control.

II. FROM DETERMINISTIC TO STOCHASTIC OPTIMIZATION

Let us consider a general constrained least-square optimization problem of the form:

minimize x ||Ax -a|| 2 subject to Bx + b ≥ 0 (1)
In the next section we will see how this relates to the problem of controlling a robot with TSID. In robotics typically (1) has infinite solutions (because A is a fat matrix) so the determined solution depends on the technique we use to solve it. The most common approach to solve problem (1) in robotics is through active-set methods [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF], mainly because they are easy to warm-start 2 (contrary to interior-point methods). Take for instance the 2D example depicted in Fig. 1. The activeset algorithm first computes a solution of the unconstrained problem (i.e. red dot). Since this solution violates an inequality constraint, it adds this constraint to the so-called active set, which is the set of constraints that are satisfied as equalities at the optimum. The new solution is represented by the blue dot. Clearly this solution has little robustness because infinitely-small changes in x, B or b could lead to violations of the active inequality constraint. Because of their working principle, active-set algorithms tend to find the worst solutions in terms of robustness. Intuition suggests that we could instead choose a solution that has a higher chance to satisfy the inequalities by moving towards the internal part of the feasible solution space.

A. Introducing Random Variables

Our idea is to model uncertainties as Gaussian noise e ∼ N (0, Σ) (with a diagonal covariance matrix

3 Σ = diag( σ 2 1 . . . σ 2 n ))
affecting the problem's variable x:

minimize x ||A(x + e) -a|| 2 subject to B(x + e) + b ≥ 0 (2)
Since e is a random variable, both cost and constraints are now random variables, so (2) does not make sense. Rather than minimizing the cost function we can minimize its expected value, but since e has zero mean, this actually does not change the problem:

I E||A(x + e) -a|| 2 = ||Ax -a|| 2
The inequalities are less trivial and consequently less frequently considered. A classic approach is to replace them with their probability to be satisfied [START_REF] Henrion | Introduction to Chance Constrained Programming[END_REF]:

p(x) = P(B(x + e) + b ≥ 0) (3)
Rather than directly using p(.) in our optimization problem, we define another function, R(.), which is easier to optimize than p(.) 4 . Once we have defined R(.), we can use it in three different ways. 1) Priority to cost: We can maximize robustness in the null space of the cost:

minimize x∈S R(x) S = {x|Ax = Ax * }, (4)
where x * is a solution of the non-robust problem (1). This approach does not affect the cost, it only exploits its redundancy to maximize the robustness of the inequality constraints. This means that the solution may have little robustness if that is necessary to get a small cost.

2) Priority to robustness: We can minimize the cost subject to the constraint of ensuring a minimum robustness R min :

minimize x ||Ax -a|| 2 subject to R(x) ≤ R min (5)
Then of course we could still maximize robustness in the null space of this minimization. To select R min so that the constraint is feasible we might first minimize R(x) to find its minimum value.

3) Trade off cost and robustness:

minimize x ||Ax -a|| 2 + wR(x) subject to Bx + b ≥ 0, (6)
where w ∈ R weighs the importance of robustness with respect to cost. Here, we kept the deterministic inequalities to prevent the solution from violating them in favor of minimizing the cost (which may happen if w is not large enough).

Regardless of which of these three formulations we use, we need to evaluate the cumulative density function (cdf) of the multivariate random variable e B = Be, that is P(e B > -b -Bx). In general there is no analytical expression to compute this cdf, so we must resort to computationally-expensive numerical techniques [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF]. This makes the resolution too slow for applications in control. Aside from introducing robustness in TSID, the main contribution of this paper is to propose three approximations of (3) that are much faster to compute and that provide satisfying precision and robustness in practice.

B. Approximation 1 -Individual Constraints p ind

The first way to simplify (3) is by considering the probabilities of the single inequalities rather than the probability of all of them:

p ind (x) = m i=1 P(B i (x + e) + b i ≥ 0), (7)
where B i is the i-th row of B. When e is Gaussian, this is equivalent to neglecting the off-diagonal terms of the covariance matrix of e B . This approximation fastens to resolution of (3) because we only need m univariate cdfs -rather than one multivariate cdf. To get an intuition of why p ind (x) is a good approximation of p(x) let us look at a simple 2D example. Fig. 2a depicts the probability p(x) to satisfy a set of 5 linear stochastic inequalities, with e having a standard deviation σ 1 = σ 2 = 1.44. Fig. 2b shows the approximated probability p ind (x) obtained with [START_REF] Sentis | Synthesis of wholebody behaviors through hierarchical control of behavioral primitives[END_REF], while Fig. 2d shows the probability of a single inequality. The overall shapes of the approximated and the real probability are quite similar and it is hard to spot the differences. To highlight the errors Fig. 2c shows the difference between p(x) and p ind (x). The errors are concentrated at the intersections of the inequalities: when the angle between the inequalities is < 90 • the error is negative, when the angle is > 90 • the error is positive, whereas when the angle is exactly 90 • the error is void. Fig. 2: 2D example comparing p(x), namely the probability of satisfying a set of affine stochastic inequalities, with its approximation p ind (x) [START_REF] Sentis | Synthesis of wholebody behaviors through hierarchical control of behavioral primitives[END_REF].

C. Approximation 2/3 -Largest-Enclosed Hyper-Rectangle p rect /p box In the first approximation we exploited the fact that it is easy to compute the probability of a single inequality. Another case in which we can easily compute the probability is when all the inequalities are simple bounds (i.e. they define a hyperrectangle aligned with the main axes). In this case, the joint probability is the product of n probabilities of univariate random variables (no more coupling). Our idea is then to approximate the real polyhedra with a hyper-rectangle U (s) (where s is a parametrization of the hyper-rectangle) that is enclosed in it:

p(x) ≈ P((x + e) ∈ U (s)) (8)
Of course for any U (s) enclosed in the feasible set, the probability to be in U (s) is lower than the probability to be in the feasible set p(x). It follows that, among all the enclosed hyper-rectangles, the one resulting in the best approximation of p(x) is the one that maximizes P((x + e) ∈ U (s)):

p rect (x) = maximize s P((x + e) ∈ U (s)) subject to Bz + b ≥ 0 ∀z ∈ U (s), (9)
where the (infinite) constraints ensure that U (s) is enclosed in the feasible set. Fig. 3 shows the hyper-rectangle maximizing the probability for different parametrizations s. In Fig. 3a U (s) is a general hyper-rectangle and hence, we need 2n values to represent it. In Fig. 3b U (s) is a hyper-cube so we can represent it with only 1 scalar value (in this case we refer to the probability approximation as p box ). Clearly these different parametrizations trade off the quality of the approximation for its simplicity. Fig. 3c and3d show the value of p rect and p box over the solution space, whereas Fig. 3e and3f show the approximation errors (p(x)p rect (x)) and (p(x)p box (x)). While these approximations may seem much coarser than the first one, in Section VII we will show empirically that they perform well in practice. Moreover, in Section V we will prove that the single-variable parametrization results in a linear optimization problem, which is easier to solve then the nonlinear problem resulting from our first approximation.

III. ROBUST TASK-SPACE INVERSE DYNAMICS

Various formulations of the TSID optimization problem exist and are often equivalent [START_REF] Del Prete | Prioritized Motion-Force Control of Constrained Fully-Actuated Robots[END_REF]. We use here the formulation of [START_REF] Saab | Dynamic Whole-Body Motion Generation Under Rigid Contacts and Other Unilateral Constraints[END_REF], written as an optimization problem of y = v f τ :

minimize y ||Ay -a|| 2 subject to By + b ≥ 0 J c 0 0 M -J c -S   v f τ   = -Jc v -h , (10)
where v are the base and joint accelerations, f are the contact forces, τ are the joint torques, J c is the constraint Jacobian, M is the mass matrix, h contains the bias forces and S is the selection matrix. The inequality constraints (defined by B and b) contain the torque limits, the (linearized) friction cones, the ZMP bounds and the joint acceleration limits. The bounds of the joint positions and velocities are typically converted into joint acceleration bounds [START_REF] Rubrecht | Constraints Compliant Control : constraints compatibility and the displaced configuration approach[END_REF]. The cost function represents the error of the task, which is typically an affine function of v (i.e. a task-space acceleration):

Ay -a = J 0 0 A y -(ẍ des -Jv) a = ẍ -ẍdes
The task may be to track a predefined trajectory of the center of mass of the robot, of a link, or to regulate its angular momentum. This problem is rather similar to the one we considered in the previous section, apart from the fact that it has equality constraints. The uncertainty e ∈ R n in the joint torques prevents us from selecting a value of y that satisfies the equality constraints. For this reason, we need to reformulate [START_REF] Henrion | Introduction to Chance Constrained Programming[END_REF] with respect to τ alone by expressing v and f as functions of τ 5 :

  v f τ   =   M -1 N c S Λ c J c M -1 S I   C τ +   -M -1 (N c h + J c Λ c Jc v) Λ c (J c M -1 h -Jc v) 0   c ,
where

Λ c = (J c M -1 J c ) -1 and N c = I -M -1 J c Λ c J c .
Then the problem takes on the following form:

minimize τ ||Dτ -d|| 2 subject to Gτ + g ≥ 0, (11)
where

D = AC, d = a -Ac, G = BC, g = Bc + b.
Introducing stochasticity we measure the robustness of a solution of (11) by some approximations of:

p(τ ) = P(G(τ + e) + g ≥ 0)
Among the three discussed options to introduce robustness in the optimization (i.e. ( 4), ( 5) and ( 6)) we decided to look for a trade-off between performances and robustness, which is:

minimize τ ||Dτ -d|| 2 + wR(τ ) subject to Gτ + g ≥ 0, (12)
where w ∈ R weighs the importance of robustness with respect to performances.

IV. APPROXIMATION 1 -INDIVIDUAL CONSTRAINTS p ind

Our first idea to approximate p(x) is to consider the constraints individually: 5 In this paper we assume that Jc is full row rank, but these results can be extended to the case of Jc being rank deficient [START_REF] Del Prete | Partial Force Control of Constrained Floating-Base Robots[END_REF].

p(τ ) ≈ p ind (τ ) = m i=1 P(G i (τ + e) + g i ≥ 0)
While most distributions have an analytical expression to compute the cdf in the univariate case, the Gaussian distribution does not. However, expressions exist to approximate it with high accuracy and low computational cost [START_REF] Vazquez-Leal | High accurate simple approximation of normal distribution integral[END_REF] (e.g. polynomials). To compute p ind we need to evaluate P(G i (τ + e) + g i ≥ 0). Since e is Gaussian, we have e Gi = G i e ∼ N (0, σ Gi ), where

σ Gi = σ 2 i G i G i . Hence: P(G i (τ + e) + g i ≥ 0) = P(e Gi ≥ -G i τ -g i ) = = P(e Gi ≤ G i τ + g i ) = F Gi (G i τ + g i ),
where F Gi (.) is the cdf of e Gi . We then define the robustness function as:

R ind (τ ) = -log p ind (τ ) = - m i=1 log F Gi (G i τ + g i )
This function is convex and twice differentiable, so we can easily minimize it using any variant of Newton's method [START_REF] Boyd | [END_REF] (see the Appendix for the expressions of the gradient and the Hessian of R ind ). The final robust TSID problem is then a convex optimization:

minimize τ ||Dτ -d|| 2 -w m i=1 log F Gi (G i τ + g i ) subject to Gτ + g ≥ 0 (13) V. APPROXIMATION 2 -LARGEST-ENCLOSED HYPER-RECTANGLE p rect /p box
Our second idea is to approximate the polytope defined by the inequalities with a hyper-rectangle. We can compute this approximation by solving this optimization problem:

p rect (τ ) = maximize s P(e ∈ U (s)) subject to G(τ + z) + g ≥ 0, ∀z ∈ U (s) (14)
First we will define U (s) for a general hyper-rectangle, and then we will show how things simplify for the hyper-cube. When U is a general hyper-rectangle (aligned with the main axes) and consequently s ∈ R 2n , we define U (s) as:

U (s) = {z ∈ R n : s i+n ≤ z i ≤ s i i = 1 . . . n}
According to this definition, the first n values of s are the upper bounds, while the second n values are the lower bounds of the hyper-rectangle. The fact that U (s) is aligned with the main axes and that Σ is diagonal results in this simplification:

P(e ∈ U (s)) = m i=1 P(s i+n ≤ e i ≤ s i ) = = m i=1 (F i (s i ) -F i (s i+n )),
where F i is the cdf of e i . Despite this simple way of computing the probability of the inequalities, using p rect (x) in place of p(x) may seem rather complex. First, we need to optimize a function (p rect ) that is itself the solution of an optimization problem and second, ( 14) cannot be solved in this form since it has an infinite number of constraints. Despite appearances, this is not more complex to solve than our first approximation; on the contrary, for a particular choice of the parametrization of U , this boils down to a simple linear program. Let us now show how we can achieve this.

A. Reduction of the Infinite Number of Constraints

We can represent the infinite constraints of ( 14) as a finite number of constraints:

l i (s) ≥ 0 i = 1 . . . m,
where l i is the solution of an optimization problem:

l i (s) = minimize z G i (τ + z) + g i subject to s i+n ≤ z i ≤ s i i = 1 . . . n
In layman's terms, we are saying that if (and only if) an inequality is satisfied for the minimum value of its LHS (over all of the possible uncertainties), then it is satisfied for all of the possible uncertainties. Thanks to the simple shape of U that we selected, this is a Linear Program with solution:

l i (s) = G i τ + G np i s + g i ,
where

G np i = G n i G p i , G p i ∈ R n contains
only the positive elements of G i and zero elsewhere, while G n i ∈ R n contains only its negative elements. The rationale behind this (apparently magic) simplification is that we do not check that an inequality is satisfied for all the values inside the hyperrectangle: we only verify that it is satisfied for its worst corner. The worst corner is the one that will eventually collide with the hyper-plane defined by the inequality if you enlarge the hyper-rectangle. By the negative elements of G i for the product with the upper bounds of U and the positive elements of G i for the product with the lower bounds of U , we are selecting the worst corner. This allows us to reformulate [START_REF] Del Prete | Partial Force Control of Constrained Floating-Base Robots[END_REF] as a convex optimization problem with linear constraints whose solution will be our robustness measure:

R rect (τ ) = minimize s - m i=1 log(F i (s i ) -F i (s i+n )) subject to Gτ + G np s + g ≥ 0 s i > s i+n i = 1 . . . n, (15)
where we constrained each upper bound of U (s) to be greater than the lower bound to be able to evaluate the logarithm 6 .

B. Robust TSID

Now that we got rid of the infinite number of constraints, we need to understand how we can minimize p rect (τ ) with respect to τ , p rect (τ ) being the solution of an optimization problem. The answer is surprisingly simple: we perform both optimizations at the same time, which gives us: [START_REF] Boyd | [END_REF] This means that we look at the same time for the solution of the original problem τ and for the "best" enclosed hyperrectangle. This is again a convex optimization problem with linear constraints. To solve it with a Newton's method, we need the gradient and the Hessian of the cost function c(τ, s), which we can compute in a few microseconds thanks to their analytical expressions (see Appendix).

minimize τ,s 1 2 ||Dτ -d|| 2 -w m i=1 log(F i (s i ) -F i (s i+n )) subject to Gτ + G np s + g ≥ 0 s i > s i+n i = 1 . . . n (

C. Single-Parameter Hyper-Rectangle p box

If we parametrize U (s) with a single variable s ∈ R we can then get a remarkable simplification in the resulting optimization problem. Let us define U (s) as:

U (s) = {z ∈ R n : -k i s ≤ z i ≤ k i s i = 1 . . . n},
where k ∈ R n encodes the fixed ratio between the n sides of the hyper-rectangle. Contrary to s, k is fixed and given by the user, so it is not a variable of the optimization. In our tests we have always set k = σ 1 . . . σ n . The robustness measure then becomes the solution of this linear program:

R box (τ ) = minimize s -s subject to Gτ -|G|ks + g ≥ 0,
where |G| is a matrix containing all the absolute values of the elements of G. This simplification is possible thanks to the fact that maximizing the probability over the hyper-rectangle is equivalent to maximizing its volume, which is in turn equivalent to maximizing s. The robust TSID optimization problem is then a simple Quadratic Program:

minimize τ,s ||Dτ -d|| 2 -ws subject to Gτ -|G|ks + g ≥ 0 (17)
To solve this problem we can use a standard QP solver, as is usually the case for the classic TSID.

VI. RELATED WORKS

Considering the robustness of the control scheme is a well-identified problem, but it remains largely unanswered. The classic approach is to consider additive noise disturbing the objective function. When the noise is Gaussian and the objective is quadratic, we can solve the problem with the linear-quadratic-Gaussian regulator (LQG), which has been exploited for animating complex robots [START_REF] Li | Iterative linearization methods for approximately optimal control and estimation of non-linear stochastic system[END_REF][START_REF] Rawlik | On stochastic optimal control and reinforcement learning by approximate inference[END_REF][START_REF] Del Prete | Prioritized Optimal Control[END_REF]. However, considering uncertainties only in the objective cannot ensure any margin [START_REF] Doyle | Guaranteed margins for LQG regulators[END_REF] and is therefore of limited interest when the system is constrained by strict limitations. Modeling noise in the inequality constraints leads to resolution times that are too long for robot control. The attempts to address this issue, mostly in trajectory optimization, have led to results of limited practical interest. In [START_REF] Mombaur | Open-loop stable running[END_REF][START_REF] Westervelt | Feedback control of dynamic bipedal robot locomotion[END_REF], the robustness is expressed by considering the stability of the resulting trajectories inside a locomotion cycle. In [START_REF] Mordatch | Ensemble-CIO : Full-Body Dynamic Motion Planning that Transfers to Physical Humanoids[END_REF], several variations of the uncertain parameters are simultaneously considered, resulting in a trajectory that is valid for any of the variations. Despite their clear interest, none of these approaches could nowadays be used for online control.

TSID is strongly influenced by parameters that are difficult to identify (e.g. body inertias, joint frictions). Adaptive control [START_REF] Kelly | On adaptive impedance control of robot manipulators[END_REF] and time-delay estimation [START_REF] Jin | Robust compliant motion control of robot with nonlinear friction using time-delay estimation[END_REF] try to correct online the estimation of the major dynamic effects, while other schemes do not make use of the dynamic model of the robot, and therefore are robust to misestimations [START_REF] Dietrich | Multi-objective compliance control of redundant manipulators: Hierarchy, control, and stability[END_REF]. However, these schemes cannot consider strict constraints: in other words, they cannot provide any margin or guarantee on the quality of the control.

In most practical cases, the objective assigned to the robot leaves redundancy in the choice of the decision variables. We suggested to use this redundancy to select the most robust solution. Originally, TSID [START_REF] Khatib | A unified approach for motion and force control of robot manipulators: The operational space formulation[END_REF] selects the motion variables answering to the principle of least action, i.e. the closest to the free fall [START_REF] Bruyninckx | Gauss' principle and the dynamics of redundant and constrained manipulators[END_REF], which neither stabilize the system nor avoid violation of the constraints. A common heuristic is to impose a secondary objective to keep the robot posture close to a reference one [START_REF] Sentis | Synthesis and control of whole-body behaviors in humanoid systems[END_REF]. Even if this heuristic is rarely commented, it tends to keep the movement away from the boundaries (typically, joint limits) and therefore it increases robustness. Few works yet exist on handling the force redundancy, which only sometimes exists (e.g. in humanoid doublesupport stances, infinitely-many torques produces the same joint accelerations [START_REF] Sentis | Compliant control of multicontact and center-of-mass behaviors in humanoid robots[END_REF]). Once more, a heuristic exists that tends to produce robust contact forces [START_REF] Righetti | Inverse Dynamics With Optimal Distribution of Ground Reaction Forces for Legged Robots[END_REF].

VII. SIMULATIONS

In this section, we present a series of simulation results that try to answer to the following questions:

• What improvement can we get in terms of probability to satisfy the inequalities by using robust TSID? • Which of the proposed formulations performs better?

• Can we solve these optimization problems in under 1 ms?

• Is it worth to use robust TSID even in case of unmodeled uncertainties, such as errors in the inertial parameters? All tests concern the control of the 30-degree-of-freedom humanoid robot HRP-2, which is always standing on its two feet in plane rigid contact with the environment. Table I lists all the simulation parameters. If not stated otherwise, the inequality constraints are: 22 constraints for the contact forces at the two feet (friction cones, centers of pressure, minimum normal force), 60 constraints for the torque limits (upper and lower bounds), and 60 constraints for the joint acceleration limits (these include also the joint position and velocity bounds). For the simulations we integrated the equations of motion of the system with a first-order Euler scheme with fixed time step dt. We assumed that the joint torques were affected by Gaussian noise with standard deviation σ proportional to the relative maximum torque τ max .

A. Test 1 -Comparing Probability Approximations

This test aims to compare the different approximations of the probability to satisfy a set of linear inequalities subject to additive noise on the decision variables. We performed this comparison on the TSID inequality constraints, since we are actually interested in how the proposed approximations perform on this particular problem. We generated a state trajectory (i.e. q, v) by controlling the motion of the CoM of the robot with classic TSID [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF]. For each state, we computed the probability p(τ ) of the joint torques to satisfy the inequality constraints. We purposely asked for a demanding motion of the CoM (20 cm in 1.6 seconds), which caused several constraints to be saturated, so that p(τ ) covered the whole range 0-100 (see Fig. 4). We then compared p(τ ) with the three approximations p ind , p rect and p box . While p ind is always quite close to p, p rect and p box are often far below p. The average error |pp ... | is 2.6 for p ind , 60 for p rect and 68 for p box . However, as we will see in the next tests, even maximizing an inaccurate approximation of p can lead to great improvements in robustness.

B. Test 2 -Comparing Robustness

The goal of this second test is to compare different TSID formulations in terms of robustness of the inequality constraints. To do so, we used the same state trajectory generated for Test 1. For each state we solved the associated control problems using five TSID formulations and we measured p(τ ).

The main task was to follow a minimum-jerk reference 3D trajectory with the CoM while maintaining as close as possible to the initial joint posture 7 . For the robust formulations, the secondary task was to optimize robustness, i.e. to minimize R(τ ). We summed the two cost functions together using a weight w for the robustness task. For the classic TSID formulation we also added the force-regularization task proposed in [START_REF] Righetti | Optimal distribution of contact forces with inverse dynamics control[END_REF] (which aims to minimize tangential forces and moments).

Finally the task was:

  J com 0 0 w post S 0 0 0 w f A f orce 0     v f τ   =   ẍdes com -Jcom v w post qdes j 0   ,
where:

A f orce = diag( w f orce w f orce ) ẍdes com = ẍref com + K com p (x ref com -x com ) + K com d ( ẋref com -ẋcom ) qdes j = K post p (q ref j -q j ) -K post d qj
These are the five formulations that we tested:

• Classic TSID with w f = 0 (11)

• Classic TSID with w f = 10 -4 (11)

• Robust TSID with p ind approximation [START_REF] Rubrecht | Constraints Compliant Control : constraints compatibility and the displaced configuration approach[END_REF] • Robust TSID with p rect approximation [START_REF] Boyd | [END_REF] • Robust TSID with p box approximation (17) 7 The postural task is necessary to maintain stability [START_REF] Peters | A unifying framework for robot control with redundant DOFs[END_REF]. Table II and Fig. 5 report the results. In terms of probability to satisfy the inequalities the three robust formulations greatly outperform the two classic formulations. The force regularization (i.e. w f > 0) slightly improves the overall probability. Among the robust formulations, the optimization of p box leads to a probability slightly lower than the others, which we expected because of its simplicity. Robust and classic formulations differ the most in the probability of the force inequalities. All the formulations lead to small errors for the CoM task (< 10 -3 ). Interestingly, the p rect formulation took on average more iterations and computation time (about ×20) than the p ind formulation, while achieving similar results 8 .

C. Test 3 -Comparing Computation Time and Overall Behavior

The previous test focused on comparing the robustness of the solutions found by different TSID formulations for the same robot state. Now we focus instead on the robustness of the whole generated trajectory and the computation time. The main task was to reach a point far in front with the right hand (see Fig. 6). To avoid falling we constrained the capture point of the robot to lie inside the support polygon [START_REF] Oscar E Ramos | Whole-body Motion Integrating the Capture Point in the Operational Space Inverse Dynamics Control[END_REF]. We ran a simulation for each solver, in which we used its solution (i.e. τ ) to simulate the system and get its new state.

Table III reports the average probabilities p(τ ) for each test. Fig. 7 shows that robust TSID with p ind approximation leads always to a higher probability than the others. Robust TSID with p box approximation leads to a higher probability than classic TSID. The task-tracking error is kept small by all the Fig. 7: Test 3. Probability to satisfy all the inequality constraints p(τ ) using classic and robust TSID.

formulations, so the gains in robustness do not lead to a loss in performances.

D. Warm-start

To speed-up the resolution we exploited the warm-start capabilities of qpOases [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], the active-set QP solver that we used. For the nonlinear problem we implemented a Sequential Quadratic Programming (SQP) algorithm [START_REF] Wright | Numerical optimization[END_REF]. We initialized the SQP search with the last solution, which most of the times led to convergence in one Newton's iteration. We used a line-search algorithm that enforces strong Wolfe conditions [START_REF] Wright | Numerical optimization[END_REF]. The algorithm stopped as soon as the squared newton decrement [START_REF] Boyd | [END_REF] was less than the desired accuracy (∆x newton H∆x newton < 2 accuracy ) or the computation time exceeded the limit t max . The computation time only included the time taken by qpOases, which means that it neglects the line search and the computation of Hessian and gradient of the cost function. This choice was motivated by two facts. First, the time to solve the QP typically dominates the time taken by the other operations. Second, these operations were implemented in Python, so their computation time is much longer than it would be with a C++ implementation (which will be mandatory for its application on a real robot). Table III shows the results: thanks to the warm start we got an average speed-up of ∼ 30×. Apart from a few outliers (maybe due to the python interface of qpOases), the computation time was always below t max = 0.8 ms.

E. Test 4 -Robustness to Noise Probability Distribution

This test verifies whether we can benefit from robust TSID when the noise probability distribution is not known exactly. We repeated the same motion of Test 3 using different values for the noise standard deviations in the controller. We then computed the probability to satisfy the constraints p(τ ) using the real standard deviations. As shown in Table IV the gain in robustness with respect to the classic TSID (which led to p(τ ) = 26.04) is large regardless of the value of σ.

F. Test 5 -Robustness to Inertial Parameters

We generated ten models of HRP-2, differing from each other only for the link masses (each mass was multiplied by a random number between 0.97 and 1.03). We performed the same reaching task described in Section VII-C with classic TSID and p ind robust TSID using the corrupted models. The results reported in Table V show that robust TSID always led to a higher probability to satisfy the constraints and, most importantly, it violated the deterministic constraints only 3 times, compared to 7 times for classic TSID.

VIII. CONCLUSION

This paper extended the TSID control framework to take the robustness of the inequality constraints into account. The inclusion of additive random white noise on the joint torques gives rise to a stochastic optimization problem that is too hard to solve in few ms. To overcome this issue we proposed three ways to approximate this problem and we showed through simulations that they lead to large increases in the probability to satisfy the inequality constraints. We also empirically showed that the same holds true even if the uncertainties affecting the system are not the same considered by the controller.

While this work focused on the theoretical results and their validation, its motivations lie in the desire to control real robots. We plan to test the presented control algorithms on HRP-2 (on which we recently implemented torque control [START_REF] Del Prete | Implementing Torque Control with High-Ratio Gear Boxes and without Joint-Torque Sensors[END_REF]) and empirically measure the robustness improvements. We are also interested in testing robust TSID with strict priorities, which will require the implementation of a hierarchical nonlinear convex solver.

  Joint inequalities probability p(x) = P(B(x + e) + b > 0).

  Individual inequalities probability p ind (x).

  Difference between Fig.2a and 2b: p(x) -p ind (x); mean error 1.8%.

  Probability of a single inequality constraint P(B 0 (x + e) + b 0 ≥ 0).

  Hyper-rectangle probability error: p(x) -prect(x); mean error 4.7%.

Fig. 3 :

 3 Fig.3: 2D example comparing p(x) with its approximations p rect (x) and p box (x)[START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF].

Fig. 4 :

 4 Fig. 4: Test 1. Probability of the inequalities of the TSID control problem computed by the three approximations proposed in this paper.

Fig. 5 :

 5 Fig. 5: Test 2. Probability of feasibility p(τ ) over time of the control determined by different TSID formulations.

Fig. 6 :

 6 Fig. 6: Test 3. Snapshots of a reaching task.

TABLE I :

 I Simulation parameters.

	Symbol	Meaning	Value
	dt	Simulation time step	0.005 s
	vmax j v max j	Max joint acceleration Max joint velocity	100 rad s -2 0.5 rad s -1
	µ	Force friction coefficient	0.3
	µm	Moment friction coefficient	0.1
	f min σ	Minimum normal force Torque standard deviation	1 N 0.003 τ max
	w	Robustness weight	10 -5 (Test 1,2) / 1 (others)
	wpost	Posture weight	10 -3
	w f orce

TABLE II :

 II Results of Test 2. For each formulation we report the average values over 700 tests.

	Formulation	Classic	Classic	Robust	Robust	Robust
		TSID	TSID	TSID	TSID	TSID
		w f = 0	w f > 0	p ind	p box	prect
	Probability p(τ )	25.1	27.4	75.7	66.5	72.3
	Force prob.	28.2	31.0	86.0	72.2	81.1
	Joint-accelerat. prob.	85.2	85.2	85.2	85.3	85.9
	Joint-torque prob.	100	100	100	100	100
	Active inequalities	3.37	2.83	0.6	0.95	0.01
	Iterations	1.06	1.05	2.06	1.11	3.96
	Comput. time [ms]	0.23	0.19	0.31	0.2	5.62

TABLE III :

 III Results of Test 3. We report the average values for each tested TSID formulation.

	Formulation			Classic	Robust	Robust
				TSID	TSID	TSID
				w f > 0	p ind	p box
	Probability p(τ )			26.04	99.23	82.96
	Force probability			65.29	99.99	90.97
	Joint-acceleration probability		67.88	99.25	83.68
	Joint-torque probability		100	100	100
	Capture-point probability		59.06	99.98	96.11
	Active inequalities			1.83	0	0.76
	Mean comput. time [ms]		0.18	0.33	0.24
	Max comput. time [ms]		1.8	1.06	1.94
	Mean comput. time (no warm start) [ms]	10.06	9.36	11.55
	Max comput. time (no warm start) [ms]	22.7	23.89	40.57
	(a) 0 s	(b) 1 s	(c) 2 s		(d) 7 s

TABLE IV :

 IV Results of Test 4. Probability to satisfy the constraints p(τ ) when the noise distribution used in the controller differs from the real one σ = 0.003τ max .

	Standard deviation σ	Robust TSID	Robust TSID
	used in controller 0.001 τ max	p ind 84.23	p box 86.42
	0.002 τ max	98.05	79.3
	0.003 τ max	99.23	82.96
	0.005 τ max	95.93	84.73
	0.006 τ max	87.27	85.89

TABLE V :

 V Results of Test 5. Probability to satisfy the constraints p(τ ) when the robot model used in the controller differs from the real one.

Warm-starting the resolution of an optimization problem consists in exploiting the solution of a similar problem (which was already computed) to speed-up the computation.

All results can be generalized to any probability distribution, with arbitrary mean and nondiagonal covariance matrices, as long as we can compute its probability density function (pdf) and cumulative density function (cdf).

For instance, rather than maximizing a probability we maximize its logarithm because it is a concave function.

In practice we set s i ≥ s i+n + , where is a small positive constant, e.g. 10 -6 .

We could reduce the computation time of the prect formulation by exploiting the sparsity of its Hessian.

ACKNOWLEDGMENTS

This work was supported by Euroc FP7 Grant Agreement 608849 and Entracte ANR grant agreement 13-CORD-002-01. The authors would like to thank Dr. René Henrion for his introductory course on chance constraint optimization.


Micha Ël Lauer
email: michael.lauer@laas.fr

Fr Éd Éric Boniol

Claire Pagetti
email: claire.pagetti@onera.fr

J Ér Ôme Ermont

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

Jérôme Ermont
email: jerome.ermont@enseeiht.fr

End-to-end latency and temporal consistency analysis in networked real-time systems

Keywords: real-time systems, embedded systems, worst-case analysis, latency, temporal consistency, mixed integer linear programming, MILP

Critical embedded systems are often designed as a set of real-time tasks, running on shared computing modules, and communicating through networks. Because of their critical nature, such systems have to meet strict timing properties. To help the designers to prove the correctness of their system, the real-time systems community has developed numerous approaches for analysing the worst case scenarios either on the processors (e.g., worst case response time of a task) or on the networks (e.g., worst case traversal time of a message). These approaches provide results only for local components behaviours. However, there is a growing need for having a global view of the system, in order to determine end-to-end properties. Such a property applies to functional chains which describe the behaviour of sequences of tasks. We propose an approach to analyse worst case behaviour along functional chains in critical embedded systems. It is based on mixed integer linear programming (MILP) and is general in the sense that it can be applied to a variety of end-to-end properties. This paper focuses on two essential properties: end-to-end latency and temporal consistency.

Introduction

Problem description

Most embedded systems rely on distributed architectures whatever the applicative domains. For instance, in the automotive domain [START_REF] Navet | A review of embedded automotive protocols[END_REF], embedded architectures are composed of several electronic control units (ECUs) connected via a several communication networks such as controller area network (CAN), FlexRay, or time triggered protocol (TTP). The software running on each ECU is often composed of real-time tasks scheduled following a fixed priority policy. In the civil aeronautics domain, embedded architectures follow the integrated modular avionics (IMA) standard (ARINC 653, 1997). An IMA architecture relies on a platform composed of modules (also knows as CPMs -core processing modules) communicating through an avionics full dupleX (AFDX switched ethernet) network. The scheduling is based on the reservation of static pre-defined time windows for the execution of each task. In the aerospace or military aeronautics domain, embedded architectures involve shared processors with fixed priority scheduling (under the RTEMS operating system for instance) interconnected with SpaceWire or Mil-1553B buses. These embedded architectures host highly critical functions such as X-by-wire functions (where 'X' denotes any safety application, such as steering, braking or guidance). Typically, these functions have to meet hard real-time requirements in order to interact correctly with their environment (i.e., others systems, sensors, actuators, human pilot, etc.). Failure to meet these requirements may lead to hazardous or even catastrophic situations. Hence, real-time analysis is a major issue for the design of safe and reliable embedded systems.

In general, the real-time analysis is decomposed in three steps:

1 Verification of the temporal behaviour of each task. This is done by proving that in any case and in particular in a worst case scenario, an execution always lasts within some time intervals.

2 Evaluation of the network worst case traversal times for each message crossing the network.

3 The combination of the last two analyses to verify end-to-end properties.

The first and second steps are already abundantly addressed in the literature [START_REF] Sha | Real time scheduling theory: a historical perspective[END_REF]. The first step is usually done by combining a worst case execution time (WCET) analysis of a task, with a worst case response time (WCRT) analysis in order to take into account the interference of the other tasks executing on the same module. The second step can be treated with various formal methods, depending on the nature of the network. For instance, the trajectory approach [START_REF] Martin | Worst case end-to-end response times of flows scheduled with FP/FIFO[END_REF] has been successfully applied to AFDX networks in [START_REF] Bauer | Applying and optimizing trajectory approach for performance evaluation of AFDX avionics network[END_REF]. The network calculus (Le [START_REF] Boudec | Network calculus[END_REF] has been extended to multi-priority switched ethernet network (Sofack and Boyer) as well. Similarly, [START_REF] Herpel | Stochastic and deterministic performance evaluation of automotive can communication[END_REF], [START_REF] Chokshi | Performance analysis of flexray-based systems using real-time calculus, revisited[END_REF] and [START_REF] Ferrandiz | Worst-case end-to-end delays evaluation for spacewire networks[END_REF] propose methods to evaluate lower and upper bounds of communication delays in other classical real-time networks such as CAN, Flexray, and SpaceWire. Generally speaking, these methods involve the communication path parameters (route in the network, throughput of the network nodes, maximum size of the messages allocated to the path, . . . ), and they determine for each path its minimal and maximal traversal times. We do not describe these analysis techniques in the following. Readers interested in worst case traversal times analysis are invited to consult the provided references.

In this paper, we focus on the third step by considering two specific properties: end-to-end latency and end-to-end temporal consistency along a functional chain. A functional chain a0 → τ 1 a1 → . . . τ R aR → (task τ 1 receives data a 0 and produces a 1 for task τ 2 . . . until a final output a R ) is a sequence of tasks deployed on computing modules, and communicating through a network or a set of buses.

Latency analysis

The end-to-end latency of a functional chain f = a0 → τ 1 a1 → . . . τ R aR → is defined as the time elapsed from an input event a 0 of f to the related response of the system a R (the one related to the received a 0 ). Usually a requirement regarding an end-to-end latency is expressed as a δ-latency requirement: for any occurrence of a 0 , the reaction a R must always be produced before a time δ at the latest. Noting t 0 the arrival time of an occurrence of a 0 and t R the related production of a R , then we must always have t R ≤ t 0 + δ. This is equivalent to impose that the maximal end-to-end latency of the functional chain f is less than δ.

For instance, a fly-by-wire control system must apply the pilot's orders to the flight control surfaces within a maximum delay depending on the speed of the aircraft. For civil aircraft, this maximum delay is about 50 ms. For military aircraft, it is about 5 ms. Such a requirement is essential for embedded systems and has to be carefully verified. However, because of the distributed nature of the fly-by-wire control function, and because of the increasing complexity of the embedded architecture (i.e., the increasing number of functions running in parallel on the architecture), analyse end-to-end latency becomes a difficult challenge for realistic systems.

Temporal consistency analysis

The second type of property is the notion of temporal consistency between outputs. Informally, if an input results in multiple outputs, then we need to guarantee that the outputs are temporally consistent, i.e., all outputs must occur within a certain time window. A little more formally, let us consider two functional chains

f 1 = a0 → τ 1 a1 → . . . τ R a R → and f 2 = a0 → τ ′ 1 a ′ 1 → . . . τ ′ R a ′ R →.
If those chains belong to a critical part of the system, they may have to satisfy a temporal δ-consistency requirement: whenever an instance of a 0 arrives at t (note that a 0 belongs to the two chains), then t aR the date of a R (related to the received a 0 ) and t a ′ R the date of a ′ R (related to the same received occurrence of a 0 ) must satisfied | t aR -t a ′ R |≤ δ. For instance, a fly-by-wire control system must apply the pilot's orders to the two opposite ailerons within a maximum time window the length of which is 3 ms (for civil aircraft). Put differently, the dates at which the order is applied on each side of the aircraft must be temporally separated by at most 3 ms.

In this work, we do not consider consistency of data with regard to their values, but only with regard to time. Thus, from now on, the adjective temporal may be implicit and we will either mention temporal consistency or simply consistency. The aim of the rest of the paper is to propose an efficient formal method, based on linear programming, for verifying δ-latency and δ-consistency requirements on networked real-time systems.

Related work

Numerous works can be found in the literature on latency and worst case response time analysis. The holistic approach [START_REF] Tindell | Holistic schedulability analysis for distributed hard real-time systems[END_REF][START_REF] Spuri | Holistic Analysis for Deadline Scheduled Real-Time Distributed Systems[END_REF] has been introduced to analyse worst case end-to-end response time of whole systems. This approach can be pessimistic as it considers worst case scenarios on every component, possibly leading to impossible scenarios. Indeed, a worst case scenario for a functional chain on a component does not generally result in a worst case scenario for this functional chain on any component visited after this component. Illustration of this pessimism is given in Section 4.

The real-time calculus by [START_REF] Thiele | Real-time calculus for scheduling hard real-time systems[END_REF] (a variation of the network calculus defined by Le [START_REF] Boudec | Network calculus[END_REF] has been proposed as an efficient method to determine worst case latencies in real-time systems. However, similarly to the holistic approach, worst case end-to-end delay is computed by adding the worst case delay of each crossed element, leading to very pessimistic results. The trajectory approach [START_REF] Martin | Worst case end-to-end response times of flows scheduled with FP/FIFO[END_REF] has been developed to deal with such over-approximations but can only be applied to evaluate network traversal time. Thus it cannot be used as such to compute real-time properties along functional chains.

The authors of [START_REF] Carcenac | A formal framework for verifying distributed embedded systems based on abstraction methods[END_REF] have developed a first methodology to compute upper-bounds of end-to-end properties in a networked embedded system. To do so, they model the functional chains and the networked architecture as a set of timed automata. In order to cope with the combinatorial explosion, they propose several abstractions. However, this work suffers from two shortcomings with respect to our objective:

1 their model does not take into account the real-time behaviour of each computing module 2 the abstractions remain insufficient and not efficient enough for realistic systems.

The study of temporal consistency between data is relatively sparse in the literature.

Temporal consistency is studied in real-time databases in [START_REF] Jha | Mutual consistency in real-time databases[END_REF] but the authors do not study the notion of functional chains. This notion can be found in [START_REF] Song | Maintaining temporal consistency: pessimistic vs. optimistic concurrency control[END_REF] and [START_REF] Pontisso | Analysis of distributed multi-periodic systems to achieve consistent data matching[END_REF]. However, both papers do not deal with IMA assumptions and [START_REF] Song | Maintaining temporal consistency: pessimistic vs. optimistic concurrency control[END_REF] evaluates temporal consistency through a simulation angle and thus cannot be used to provide safe guarantees to a critical system. Works linking real-time systems and linear programming can also be found. For instance, the authors of [START_REF] Sagaspe | Constraint-based design and allocation of shared avionics resources[END_REF] propose an automatic method based on integer linear programming (ILP) for allocating functional specifications on a platform while ensuring some safety (i.e., fault-tolerance) requirements. In Al [START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF], the authors also propose a technique based on mixed ILP to automatically allocate and schedule avionics functions on an IMA platform. Their objective is the maximisation of the spare resource on each module in order to meet future resource demand growth. Also, the authors of Cucu-Grosjean and Buffet (2009) use a constraint satisfaction problem (CSP) to find a feasible real-time schedule of tasks in a multi-processor context.

Contributions

In [START_REF] Boniol | End-to-end latency analysis in networked real-time systems[END_REF] we proposed a latency analysis method for a specific class of real-time systems, in which tasks are scheduled with a preemptive fixed priority policy, and are supposed to satisfy a strong assumption: the worst and best execution times of a task are equal, i.e., the execution time does not vary. It is easy to show that under such an assumption, the system is deterministically and statically scheduled. Outputs were produced deterministically (with respect to the initial time of the task) at the end of each job. Put differently, this first work only considered deterministic systems. In the next sections, we extend [START_REF] Boniol | End-to-end latency analysis in networked real-time systems[END_REF] in two directions.

• We continue to consider statically scheduled systems in which tasks run within time intervals. However, the system model is now more general and fits a wider variety of embedded real-time systems. First, we consider non-deterministic systems: indeed, task execution time can vary from a lower bound (zero in this article, however any positive value could be used as well) to the WCET. Secondly, we consider that a task can be interrupted and resumed later in its next time interval. Formally, each job of each task has got several non-contiguous time intervals to complete its execution, contrary to [START_REF] Boniol | End-to-end latency analysis in networked real-time systems[END_REF] in which each job completely runs in a single deterministic time interval. Such an extension allows to capture more general behaviours, and some kind of non-determinism.

• The approach has been extended to analyse an other real-time property: temporal consistency. We show that such a property can be computed in the same framework, following a mixed integer linear programming (MILP) modelling.

We illustrate these contributions on an industrial case study built from the Airbus A380 configuration, and we show using several benchmarks that the method remains scalable even for temporal consistency properties and for large functional chains.

Organisation

Section 2 details the main assumptions of the real-time systems we consider in the paper and gives a model for the attributes required for the analysis. This model is also instantiated on a motivational case study taken from the aerospace domain. In Section 3 we present our approach by describing the modelling and then the computation of worst case properties. Several aspects of our approach are discussed in Section 4: the tightness of its results compared to a local approach, its scalability and its possible extension to best case analysis. Section 5 concludes the paper and opens some research directions.

The different concepts introduced in the paper are illustrated on the case study.

Model and main assumptions

In this section, we introduce a model for the systems under consideration and the assumptions we use. We also instantiate the model on a motivational case study taken from the aerospace domain.

Model

Critical embedded systems are often composed of tasks statically scheduled on shared computing resources and communicating through a shared network. This is the case for modern aircraft such as the Airbus A380 or the Boeing B787. These embedded systems follow the IMA standard ARINC 653 (1997). The scheduling on each computing module is time-triggered, meaning that each task periodically executes at fixed and predetermined time intervals. However, in order to avoid the use of complex synchronisation protocols, modules are globally asynchronous. Such systems can be considered as globally asynchronous and locally time-triggered (GALTT). Executions of the tasks on a module describe a periodic pattern. Each instance of a task in this pattern is called a job and a time interval (or a set of time intervals) is strictly reserved for its execution.

In the following we consider GALTT systems composed of N tasks T = {τ 1 , . . . , τ N } running on a set of m modules M = {M 1 , . . . , M m } communicating via a shared network. We note T (M i ) the set of tasks hosted by module M i . The assumptions made for the system under analysis are:

• Module: each module M i ∈ M is characterised by a period H i , i.e., the hyper-period of the tasks running on the module. The hyper-period is the least common multiple of the hosted tasks periods, H i = lcm(τ j ) τj ∈T (Mi) .

• Task: each task τ j ∈ T (M i ) is characterised by a set of jobs τ k j for k = 0..n j , meaning τ k j is the k th job of the task τ j in the period of M i . Each job is characterised by a set of time intervals [b k j (l), e k j (l)] for l = 0..n k j , where b k j (l) and e k j (l) are the beginning and ending dates respectively of each interval.

More precisely, the k th job of the task τ j starts at b k j (0). If it is not finished at e k j (0), then the job is suspended until b k j (1). If it is not finished at e k j (1), then the job is suspended until b k j (2), and so on until e k j (n k j ). It is assumed that a job always finishes during one of its intervals. During any interval [b k j (l), e k j (l)], only the job τ k j is executing on the module. These dates are relative to the beginning of the period of the module M i . Note that the sum of the intervals of a job is equal to the task WCET.

• TCommunication scheme: tasks communicate in an asynchronous way. Each job τ k j consumes input data arrived before the beginning of its first reserved time interval: between b k j (0) and b k-1 j (0). Inputs received after the beginning of the job will be consumed by the next job. This is consistent with typical queuing communication scheme. Each job produces output data at any time during its execution, meaning during any interval [b k j (l), e k j (l)] for l = 0..n k j . • TCommunication delay: tasks running on the same module communicate through the memory of the module, i.e., without delay. Tasks running on different modules communicate through the shared network. The communication delays in the network are supposed to be bounded.

• TGlobal asynchronism: finally, we suppose that the modules M are globally asynchronous, i.e., they can be shifted by an arbitrary amount of time.

Nevertheless, these offsets are supposed constant.

In the paper, we do not take into account any drift between the clocks of the module.

Although clock drift is a major issue in synchronised systems where a shared time reference needs to be established, in asynchronous systems, by definition, such time reference is not required. Still, the discrepancies in clock frequencies which cause the clock drift could have an impact in our analysis. Some modules could run a little faster (or slower). This may modify the actual tasks periods and executions times. However, worst case latency or consistency are usually measured in hundreds of ms. A high-quality quartz typically used in avionics systems is assumed to lose at most 10 -8 seconds per second. Hence, clock drift could not significantly impact our results. To the best of our knowledge, it is an implicit assumption in every performance evaluation papers dealing with asynchronous systems.

Motivation: an avionics case study

Let us consider an avionics case study depicted in Figure 1. This case study is a part of a flight management system (FMS).

System description

The system manages some displays in the cockpit. It provides some information on a waypoint requested by the crew (distance and estimated time of arrival). The pilot and the co-pilot use a personal keyboard and two displays to interact with the FMS. Information displayed on both screens must be similar although they are not processed by the same components. The FMS uses a redundant implementation of its functions which are segregated on each side of the plane (named side 1 and side 2). This system is composed of 18 tasks mapped onto five modules and communicating through a network (not shown in the figure). M 11 and M 12 are two redundant modules. They implement the same redundant tasks. They manage the interface, keyboard and display, for the pilot (M 11 ) and for the co-pilot (M 12 ). Each module M 1i hosts three tasks: KC i (keyboard control) which reads data sent by the pilot (if i = 1) or the co-pilot (if i = 2) through their respective navigation keyboard; MFD i (MultiFunction Display) which manages the data to be displayed on side i, and a background autotest task. Two redundant modules M 21 and M 22 host the tasks of the FMS 1 and FMS 2 : FlightM i is the main task, computing periodically the flight plan; CockpitReqM i manages the interactions between the crew and the FMS i ; it receives/sends messages from KC i / to MFD i ; WayPointM i elaborates the data of the flight plan waypoint, and it computes distance and estimated time of arrival of the waypoint requested by the pilot to be displayed; and FMAutoTest i is the background autotest of M 2i . Finally, a single module M 3 hosts the navigation database (NDB), which contains information related to the flight plan. This database is common to the two redundant sides. It is composed of four tasks: a background autotest; NDBReqM manages the request received from other modules; NDBServ processes the request; and finally NDBRep builds the response message. These tasks are scheduled on their corresponding modules as followed:

• KC i = 3 jobs in 120 ms: ([19, 31[), ([62, 74[), ([99, 111[) • MFD i = 3 jobs in 120 ms: ([0, 19]), ([43, 62[), ([80, 99[) • CockpitReqM i = 2 jobs in 120 ms: ([31, 41[), ([76, 80[, [95, 101[) ([141, 156[).

• WayPointM i = 2
For example, the task CockpitReqM i executes two jobs during the period of the module. Its first job executes in a single interval while the second has two intervals. According to the previously introduced notation, the first interval of the second job starts at b 1 CockpitReqM i (0) = 76 and it finishes at e 1 CockpitReqM i (0) = 80. Then, the second interval of the second job starts at b 1 CockpitReqM i (1) = 95 and it finishes at e 1 CockpitReqM i (1) = 101. This system is part of the whole avionics system composed of almost 100 modules communicating via more than 100,000 data through almost 1,000 communication paths.

System behaviour

Let us consider the following scenario. The pilot (on side 1) requests some information on a given waypoint. This request can be entered at any time via his navigation keyboard. KC 1 (on M 11 ) controls this keyboard. When KC 1 receives a request, it broadcasts it to CockpitReqM 1 and CockpitReqM 2 , which then query the navigation database for information on this waypoint (task NDBReqM). Each query is processed by tasks NDBServ and NDBRep (in that order). The answer (information on the waypoint) is sent to WayPointM 1 on module M 12 and to WayPointM 2 on module M 22 . Upon the reception of this message, each WayPointM i computes two complementary dynamic data: the distance to the waypoint, and the estimated time of arrival (ETA). These data are periodically sent to MFD i on module M 1i which periodically elaborates the pages to be displayed on the screen. This communication flows, from KC 1 to MFD 1 and MFD 2 via the database, form two functional chains summarised in Figure 2.

System requirements

These chains have to meet the following end-to-end temporal requirements:

• φ L : pages relative to a pilot's request must be displayed on each side within a time window of 700 ms, meaning that the duration between t req , the date of the pilot's request, and t dispi , the time at which the requested information is displayed, must satisfy t dispi -t req ≤ 700 ms for i = 1, 2

• φ C : the two pages on side 1 and 2 relative to a same pilot's request must be displayed in a bounded time window, more precisely:

| t disp1 -t disp2 |≤ 300 ms.
φ L is a worst case end-to-end latency requirement. φ C is a worst case temporal consistency requirement.

Worst case properties

In this section, we present our method to analyse temporal properties in networked real-time systems. We first model the behaviour of each element involved in the functional chain (modules, tasks, communication delays, . . . ) by a set of variables and constraints. The behaviour of the whole system is obtained by the conjunction of all these constraints. Using an objective function measuring the latency of the functional chain, we define a MILP which optimal solution gives the worst case latency of the chain. In the following, all variables used for offsets and dates are reals. All variables used to designate a specific hyper-period or a job interval are integers. We also use Boolean variables for technical reasons. Although not mandatory, we only use integers for parameters in order to improve readability.

To capture the temporal behaviour of a functional chain,

F = a0 → τ 1 a1 → . . . τ R a R
→, we need to 'follow' the propagation along the chain of an occurrence of an input a 0 until a related output a R is produced.

We then describe the case of latency in Section 3.2 and the specificities of consistency in Section 3.3.

Constraints and variables definition

We define the constraints and variables allowing to capture the expected behaviour of the modules, tasks and communications.

Module

Let M i be a module. The only variable which characterises M i is its possible offset with respect to the other modules. Modules are asynchronous, thus the time origin of their execution frame may be shifted by an offset O i . This offset may be arbitrary. However, as we are interested in the regular behaviour, and not the specific case of the initialisation phase, it is not necessary to consider offsets greater than the maximal hyper-period of the system. The first constraints for M i are then:

O i ∈ R, 0 ≤ O i ≤ max k=0...m H k (1)

Task

In the following, we focus on a single task of the chain τ j . The link with its previous and next tasks in the chain is later realised with the constraints on the communication delays (see 3.1.3). Let τ j ∈ F be a task of the functional chain executing on module M i , such as u → τ j v →, meaning the task consumes data u in order to produce data v. We first describe how to capture the job and hyper-period of the task in which a specific occurrence of u is consumed. Then, we show how to constrain the production of the occurrence of v related to this u in order to follow the propagation of u through v.

Consumption of data

Let t u be the date where an occurrence of u is available on M i . As presented in Section 2.1, depending on the location of the task preceding τ j , it may be delivered by the network, or put in a shared memory of the module. To capture the temporal behaviour along the functional chain we need to link t u to the job of τ j that consumes this occurrence of u. From the task model, we have for each start of a job of τ j : ∃n ∈ N, ∃k ∈ 0..n j such that

O i + nH i + b k j (0)
, where O i is the offset of the module M i , H i its hyper-period, n the number of the current hyper-period and b k j (0) the beginning of the first interval of the k th job of the task τ j . Hence, the occurrence of u arriving at t u is consumed by the k th job of the task τ j during the n th hyper-period of M i if and only if:

t u ≤ O i + nH i + b k j (0) { O i + nH i + b k-1 j (0) < t u , ∀k = 1..n j O i + (n -1)H i + b nj j (0) < t u , for k = 0 (2)
as shown in Figure 3 there is a special case for k = 0, i.e., the first job of the hyper-period.

t u O i + nH i + b k j (0) O i + nH i + b k 1 j (0) O i + nH i + b 0 j (0) O i + (n 1)H i + b nj j (0) ⌧ k 1 j (0) ⌧ k j (0) ⌧ k j (n k j ) ⌧ k 1 j (n k 1 j ) job ⌧ k 1 j job ⌧ k j job ⌧ nj j job ⌧ 0 j u ... ... n th hyper-period t u ⌧ 0 j (0) (n 1) th hyper-period ⌧ nj j (0) ⌧ 0 j (n 0 j ) ⌧ nj j (n nj j )
In order to take into account these constraints into the MILP, we define a set of Boolean variables C j , for all k = 0..n j (one Boolean per job, C stands for consumption). These are indicator variables: C k j is true iff job τ k j consumes the occurrence of u arrived at time t u . The Boolean variables are constrained by:

nj ∑ k=0 C k j = 1 (3)
meaning that one and only one job consumes u, and in order to respect the constraints (2), the Boolean variables are also constrained by:

t u ≤ O i + nH i + nj ∑ k=0 C k j • b k j (0) O i + nH i + nj ∑ k=1 C k • b k-1 j Infty • C 0 j < t u O i + (n -1)H i + Infty • (C 0 j -1) + b nj j < t u (4)
where Infty is a large integer representing infinity. We use this large value in order to deal with the two cases of constraints (2). Indeed, if C 0 j = 1 (meaning the consumer job is τ 0 j ) then the second constraint is always true and the only active constraint is the third one. Otherwise, if C 0 j = 0 (meaning the consumer job is τ k j , k = 1..n j ) the third constraint is always true and the second one is the only one active. As a result, for a given O i and for a given t u (arrival of u), the previous constraints determine a single n and a single i.e., the hyper-period and the job consuming u.

Production of data

We now capture the behaviour for the production of data v by task τ j following the consumption of an occurrence of u arrived at t u . Let t v be the date of this production of v. Data can only be produced during one of the job intervals. If t v happens during the n th hyper-period of M i and the l th interval of the k th job, then necessarily we have:

O i + nH i + b j (l) ≤ t v ≤ O i + nH i + e k j (l) (5)
In order to take into account these constraints into the MILP, we use a set of Boolean variables P k j (l), for all k = 0..n j and for all l = 0..n k j (one Boolean per job and per interval, P stands for Production). These are indicator variables:

P k j (l) is true iff job τ k j produces the occurrence of v at time t v during interval [b k j ′ (l), e k j ′ (l)].
Since the job producing v is the same than the one which consumed u arrived at t u , the Boolean variables are constrained by

C k j - n k j ∑ l=0 P k j (l) = 0, ∀k = 0..n j (6)
We show that the constraints (6) model the link between the job consuming u arrived at t u and the interval of the job producing v. Let k = 0..n j be a job, if C k j = 0 then necessarily, ∀l = 0..n k j , P k j (l) = 0, meaning that no interval of job k produces the occurrence of v we want to follow. Conversely, if C k j = 1 then necessarily, ∃l = 0..n k j , P k j (l) = 1, and ∀l ′ = 0..n k j , l ′ ̸ = l, P k j (l ′ ) = 0 meaning that one and only one interval of job k produces v. From constraint (4), we know that only one job k can consume u. Thus, the constraints model the expected behaviour.

We now use the Boolean variables P k j (l) to constrain t u to happen during the interval selected by the previous constraints:

O i + nH i + nj ∑ k=0 n k j ∑ l=0 P k j (l) • b k j (l) ≤ t v t v ≤ O i + nH i + nj ∑ k=0 n k j ∑ l=0 P k j (l) • e k j (l) (7)

Communication delay

The communication constraints link the tasks in the chain. We describe this link between a task τ j producing v and τ j ′ consuming v: τ j v → τ j ′ . Two cases arise for the communication delay. Either the communication between τ j and τ j ′ happens through the network or, if the tasks execute on a common module, through a shared memory. Let t out v be the time at which τ j produces data v was noted t v in the previous constraints) and t in v be the time at which v is available for consumption in the module of τ j ′ .

Communication through a timed channel

As presented previously, communications through the network are abstracted with a timed channel characterised by a communication delay in [δ min , δ max ]. Hence, t out v and t in v are constrained by:

t out v + δ min ≤ t in v t in v ≤ t out v + δ max (8)

Communication through shared memory

If the communication is realised through a shared memory, then it is considered instantaneous. Hence, we have:

t out v = t in v (9)
a shared memory is similar to a timed channel characterised by the interval [0, 0].

Objective function: worst case latency

Let F = a0 → τ 1 a1 → . . . τ R aR
→ be a functional chain. The possible behaviour of this chain is obtained by the conjunction of all the constraints of all the tasks and the communication involved in the chain. The set of constraints thus obtained forms a MILP (because some variables are reals and some are integers). Let us note t 0 the arrival date of a 0 , and t R the production date of a R . Then the latency of the chain is

L = t R -t 0
The worst case latency is obtained on a particular behaviour maximising L. This behaviour can be found by using a MILP solver with the objective function: maximise: t R -t 0

Objective function: worst case consistency

In this section, we detail the specificities of the proposed MILP formulation for the worst case consistency analysis. We consider a set of functional chains F 1 , . . . , F K .

For all c = 1..K, we note R c the number of task in F c . Every functional chain has the following form:

F c = ac,0 → τ c,1 ac,1 → . . . ac,R c-1 → τ c,Rc ac,R c → (10)
We are interested in the temporal consistency of divergent chains, thus all chains share a common first element:

∀(c, c ′ ) ∈ [1..K] 2 , a c,0 = a c ′ ,0 .
Temporal consistency measures the distance between the outputs of functional chains relating to a common input. Let t 1,R1 , . . . , t K,RK be the dates of such outputs. All these dates depend on a common initial date t 0 which corresponds to an arrival at the input of the chains. t 0 is also a variable of the problem. Moreover, the dates are constrained with the same constraints as for the latency. The worst case consistency is then the maximum distance between the largest t c,Rc , ∀c = 1..K and the smallest one. More precisely, the objective is:

maximise: max c=1..K {t c,Rc } -min c=1..K {t c,Rc } (11)
However, this objective is not linear and cannot be used in a MILP. In the following, we propose a linearisation of this objective function. First, we define two technical variables (maxT, minT) ∈ R 2 . They represent respectively the maximum and minimum dates of the outputs. The objective becomes:

maximise: maxT -minT (12)
Then, we add some constraints to ensure that at least one of the dates is equal to maxT and one of the dates is equal to minT. Since the objective is a maximisation, the solver will choose the largest value for maxT and the smallest for minT. We start with maxT.

For each t c,Rc , c = 1..K, we add the two following constraints:

maxT ≤ t c,Rc + (1 -B maxT c ) • Infty maxT ≥ t c,Rc -(1 -B maxT c ) • Infty (13)
where B maxT c is a Boolean variable and Infty is a large integer representing infinity. If B maxT c = 1 then necessarily t c,Rc = maxT. Otherwise the variable t c,Rc is free. We also add the following constraint to ensure that at least one variable B maxT c is equal to 1:

K ∑ c=1 B maxT c ≥ 1 (14)
We proceed in a similar manner to constrain the value of minT. For each t c,Rc , c = 1..K, we add the following constraints:

minT ≤ t c,Rc + (1 -B minT c ) • Infty minT ≥ t c,Rc -(1 -B minT c ) • Infty K ∑ c=1 B minT c ≥ 1 (15)
The linearisation of the objective function results in the addition of two real variables, 2K Boolean variables and 4K + 2 constraints.

The purpose of this section is to discuss several aspects of our approach. First, we benchmark our global approach against a local approach for the latency and consistency properties. The gain we achieve is intuitively explained in a simple example. The improved tightness of the analysis comes at a price: complexity. Our global approach involves solving a MILP. It is obviously more complex than a local approach that typically involves analytical formulas. We run several experiments to test the scalability of our approach. Finally, we consider that a key feature of our work is its flexibility. We show how it can be extended to analyse other properties. As an example, we propose an extension to best case latency and best case consistency evaluations.

Local vs. global analysis

For each property, we first give results for a local approach. These results are similar to formulas used in Al [START_REF] Sheikh | Strictly periodic scheduling in IMA-based architectures[END_REF] for latency and in [START_REF] Pontisso | Analysis of distributed multi-periodic systems to achieve consistent data matching[END_REF] for consistency. However, we specialise these local approaches to our model.

Worst case latency

As described in the related works, a local approach consists in determining worst case response time (WCRT) of each component visited by the functional chain, then the end-to-end latency is the sum of the WCRT.

Local Worst Case Latency

The WCRT of a timed channel c i is the upper bound of the communication delay:

δ max i
. If we note WCRT(τ j ) the worst case response time of a task τ j , then with a local approach an upper bound of worst case latency of a functional chain F = a0 → τ 1 a1 → . . . τ R aR → which uses a set c 1 , . . . , c l of timed channels to communicate is:

WCL(F) = R ∑ j=1 WCRT(τ j ) + l ∑ j=1 δ max j (16)
In our model, the WCRT of a job τ k j of a task τ j happens when a data waits the longest possible time before being acquired by the job and the processing in the job takes the longest possible time. This scenario occurs when the data arrives just after the beginning of the first interval of the previous job, and the production happens at the end of the last interval, as shown in Figure 4 (u represents the data and v the resulting data of its processing by τ k j ). Taking into account the hyper-period of the module M i hosting the task, the WCRT of this job is then:

( e k j (n k j ) -b p(k) j (0) ) mod H i , p(k) = { k -1 if k > 0 n j if k = 0 (17)
where p(k) the index of the job preceding job k.

Figure 4 Local worst case scenario for the job τ k j (see online version for colours)

⇣ e k j (n k j ) b p(k) j (0) ⌘ mod H i b k j (n k j ) e k j (n k j ) b p(k) j (0) ⌧ k j (n k j ) ⌧ k j (0) ⌧ p(k) j (n p(k) j ) u v ... ... job ⌧ p(k) j job ⌧ k j ⌧ p(k) j (0)
The worst case response time of a task τ j is the maximum response time of all its jobs: To be fair with the local approach, we use the following idea: because of the scheduling of the navigation database tasks, the result of the first (resp. second) job of NDBReqM is necessarily processed in the first (resp. second) job of NDBServ and which result is then processed in the first (resp. second) job of NDBRep. Thus, this sequence of tasks can be considered as only one task with two jobs per hyper-period. The beginning of the first (resp. second) job is the beginning of the first job (resp. second) of NDBReqM and the end of the first job (resp. second) job is the end of the first job (resp. second) of NDBRep. In the following, we name this task NDB and it is characterised by the intervals: [0, 54], [102,156]. Table 1 provides the WCRT considered for the tasks.

WCRT(τ j ) = max k=0..nj {( e k j (n k j ) -b p(k) j (0) ) mod H i } (18)

Local vs. global worst case latency

We compare the global approach against the local one by varying the upper bound of the communication delays (for simplicity we assume that all the δ max i are equals). The results are plotted in Figure 5. According to equation ( 16) the worst case latency determined with the local approach is linear (straight dashed line). The results of the global approach form a step function which is right continuous: for example [15,22), W CL = 443. . . We can see on the results that the global approach provides more accurate results than the local approach. The curve of the global approach varies by steps. One of the reason is that the involved functional chain crosses twice the modules M 1 and M 2 . Consider a part of the worst case scenario depicted in Figure 6. The task CockpitReqM (noted CRM) on module M 2 sends a query to the task NDB on module M 3 , the answer is then processed by the task WaypointM (noted WP) on module M 2 resulting in an output data wpInfo. We can see that variations in δ max value (up to a certain point) do not change which job of WP processes NDB answers. Thus, the end-to-end latency is unaffected by variations of the worst case communication delay. In the case of the local approach, it is the job with the worst case response time which is always assumed to process the data. Which also explains the improvement obtained by the global approach. System designers could take advantage of this more accurate evaluation technique: within certain range they could increase the network load with no impact on high-level requirement, which could be crucial for future evolution of the system.

∀δ max i ∈ [0, 15), W CL = 403, ∀δ max i ∈
To confirm the trend in this results, we extend the case study to build two experiments:

• A module M ′ 3 , similar to M 3 , is added. The tasks of M ′ 3 correspond to the navigation database tasks and as for M 3 we consider them as a bundle of tasks noted NDB'. The functional chain considered is then: • A module M ′ 2 , similar to M 2 , is added. The tasks of M ′ 2 correspond to the flight management tasks of M 2 . In particular, we add the tasks CockpitReqM' and WayPointM' to the functional chain. The goal is to add a module which is crossed twice. The functional chain considered is then: The results of these experiments are plotted in Figures 7 and8 respectively.

The improvement of the global approach is more important with the experiment 2. Indeed, the module added to the case study is crossed twice by the functional chain and so the benefit of the approach (as presented in Figure 6) is more significant. These results indicate that the more modules are crossed multiple time by the functional chain, the more the proposed approach may be relevant.

Worst case temporal consistency

An upper bound on the worst case consistency is given by the maximum difference between the worst case and best case latency of the functional chains. Thus, to express this bound we first need to define a lower bound on the best case latency of a chain.

Local best case latency

For the local approach, the best case response time (BCRT) of a job of a task happens when data arrives just before the beginning of the job and the data is processed with minimum delay. In our model, no minimum delay is defined, thus processing can be instantaneous. The best case latency is then reduced to the sum of the lower bound of the communication delays. For a functional chain F = a0 → τ 1 a1 → . . . τ n an → which uses a set c 1 , . . . , c l of timed channels to communicate, we have:

BCL(F) = l ∑ j=1 δ min j (19)
Note that a minimal processing delay could be included in our model without technical difficulty.

Local worst case consistency

With a local approach, we can obtain the following upper bound for the consistency of the functional chains F 1 , . . . , F K :

WCC(F 1 , . . . , F K ) = max c=1..K { WCL(F c ) } -min c=1..K {BCL(F c )} (20)

Local vs. global worst case consistency

As for the worst case latency, we compare the global approach against the local one by varying the upper bound δ max of the communication delays. The lower bound of the communication delays is assumed to be constant δ min = 1ms. The results are plotted in Figure 9.

Scalability

We have developed a prototype which takes the description of a system and a set of functional chains as inputs. It automatically generates the MILPs to determine the worst case latency of a chain or the worst case consistency of a set of chains. Then, the MILPs are solved with the open source solver lp solve. Experiments run on a Intel Core 2 Duo 2.53 GHz processor.

To find the limiting factors of the approach, we further extend the case study. In Figure 10(a) we add up to four modules with navigation database tasks and we plot the required computation time to the corresponding MILPs. In Figure 10(b) we add up to four modules with flight management tasks and we plot the required computation time to solve the corresponding MILPs. For the consistency, we both increase the length of the functional chains by adding NDB modules and increase the number of functional chains by duplicating the modules. The corresponding computation times are plotted in Figure 11.

K = 2 K = 3
In all these cases, the required computation time increases exponentially with the number of added modules. This is due to the fact that for each added task, n binary decision variables must be added to the MILP (with n the number of interval characterising the task). Remember that these binary variables are used to decide which job in the hyper-period of the task module is involved in the functional chain. Solving the MILP for the consistency is also more complex for two reasons:

1 Boolean variables have been added to linearise the objective function (c.f., Section 3.3).

2 The problem is highly symmetrical. This is known to be a burden for solvers.

We can also notice that solving the MILP for the experiments with a module similar to M 2 (flight management module) is more difficult. Let us consider again the scenario depicted in Figure 6. Finding the jobs of CRM and WP involved in the worst case depends on the time elapsed between the production of the query and the reception of the answer. Thus, inter-dependencies exist between part of the functional chain and the MILP becomes more complex to solve.

Extension to best case evaluation

One advantage of the proposed approach is that it is easily extendable. Indeed, it allows to evaluate other relevant properties like the minimum latency or consistency of the functional chains for instance. A classical example of system which must verify a minimum and a maximum latency is an airbag: upon a collision, it must neither be inflated too early nor too late. We first present the ideas for extending our approach to best case latency and compare the results to a local approach. Then we discuss the specificities of the best case consistency.

Best case latency

To apply our approach to the best case latency, we can reuse the MILP defined for the worst case. We just need to change the objective function to: minimise: t R -t 0 (21)

As for the worst case latency, we compare a local approach with a global approach. We use the same methodology as with the worst case to benchmark the global approach against the local one. Increasing the upper bound of communication delay does not impact the best case latency, thus only the lower bound (δ min i

) is varying. The results are plotted in Figure 12. The straight dashed line corresponds to the local approach. Because some steps of the results of the global approach are small, we do not plot the discontinuities (i.e., the jumps of the function). The function is left continuous. We can see the global approach improves the analysis of best case latency. The reasons are similar to the analysis of the worst case.

Best case consistency

The best case consistency is the minimum distance between the outputs of the functional chains: minimise: max The concept of best case consistency is less intuitive than the worst case, we give an illustration in Figure 13. The consistency of 3 functional chains F 1 , F 2 and F 3 are considered. The ranges of possible values for the outputs t 1,R1 , t 2,R2 and t 3,R3 are represented. These ranges depend on the constraints of the MILP. On this example, the best case consistency is obtained by giving the largest possible value to t 3,R3 , the smallest one to t 2,R2 and any value in between to t 1,R1 .

Figure 13 An example of best case consistency

F 1 F 2 F 3 BCC(F 1 , F 2 , F 3 ) t 1,R1
t 2,R2 t 3,R3 : : :

As for the maximisation problem, two technical variables maxT and minT are introduced to linearise the objective function: minimise: maxT -minT (23) and the same constraints (13), ( 14) and ( 15) are used. However, to apply our approach to the best case consistency, some constraints must be added to the MILP used for the worst case. Indeed, otherwise for the minimisation solver could assign the minimal value to maxT and the maximal one to minT. For each t c,Rc , c = 1..K, we add the two following constraints:

maxT ≥ t c,Rc minT ≤ t c,Rc (24)
thus, maxT is constrained to be equal to the largest t c,Rc and minT is constrained to be equal to the smallest t c,Rc . Applied to the case study, the best consistency is always null. This is because the problem is highly symmetrical and the data can be produced at the output at exactly the same time.

Conclusions and perspectives

The article presents an analysis method for end-to-end freshness and end-to-end temporal consistency properties on GALTT systems. This verification method is based on a MILP modelling. Worst case end-to-end properties are computed as optimal solutions of the MILP problem. An interesting feature of this approach is that it computes more accurate temporal bounds than local approaches as shown in Section 4.1. Another interesting point, as shown in Section 4.3, is its flexibility: one can easily compute best case end-to-end bounds by only modifying the objective function of the MILP form max to min. From a scalability point of view, the case study considered previously is composed of 11 tasks. This case study is representative from industrial systems (usually composed of 5 to 10 tasks). Our method applied to this case study does not take more than 1s for the latency or consistency analysis (with a non-optimised solver). We think that these results are promising.

In this article, we made however a strong hypothesis about the internal behaviour of the tasks. We implicitly considered that each job of each task does not induce a delay greater than its worst case response time, i.e., the end of its last time interval. Obviously it is not always the case in realistic systems. Some tasks can implement 'confirmation tests' waiting for a given amount of time (generally a multiple of its period) before producing a consolidated output. Obviously this internal latency impacts the global latency and the global consistency of the chain. Our next work is to extend our global method by tasks involving internal delays.

Figure 1

 1 Figure 1 Case-study: a flight management system (see online version for colours)

Figure 2

 2 Figure 2The functional chains of the flight management system

Figure 3

 3 Figure 3 Special case for the first job of the hyper-period (see online version for colours)

  For example, in the case study the task WaypointM is characterised by two jobs {[15, 31]} and {[60, 76]} (only one interval per job) on the hyper-period H 2 = 120. Thus, WCRT(WaypointM) = max{(31 -60) mod 120, (76 -15) mod 120} = max{91, 61} = 91

Figure 5 Figure 6

 56 Figure 5 Results of the case study -global vs. local approach (see online version for colours)

Figure 7 Figure 8

 78 Figure 7Results with module M ′ 3 added (see online version for colours)

Figure 9

 9 Figure 9 Worst case consistency -global vs. local approach (see online version for colours)

Figure 10 Figure 11

 1011 Figure 10 Worst case latency computation time, (a) number of NDB modules (b) number of FM modules (see online version for colours)

Figure 12

 12 Figure 12 Best case latency -global vs. local approach (see online version for colours)

Table 1

 1 Worst case response time of the case study tasks

	Task	KC	CockpitReqM	NDB	WayPointM	MFD
	WCRT	55	85	156	91	62

This work was supported by the French National Research Agency within the SATRIMMAP project.