Kevin Zagalo

Yasmina Abdeddaïm

Avner Bar-Hen

Liliana Cucu-Grosjean

Response Time Stochastic Analysis for Fixed-Priority Stable Real-time Systems

In this paper, we prove that a mean system utilization smaller than one is a necessary condition for the feasibility of real-time systems. Such systems are defined as stable. Stable systems have two distinct states: a transient state, followed by a steady-state where the same distribution of response times is repeated infinitely for each task. We prove that the Liu and Layland theorem holds for stable probabilistic real-time systems with implicit deadlines, we provide an analytical approximation of response times for each of those two states and a bound of the instant when a real-time system becomes steady.

Introduction

Probabilistic methods applied to real-time systems are motivated by the fact that worst-case execution times and minimal inter-arrival times are often rare events. The analysis induced by these methods consists in quantifying the probability that time requirements may not be satisfied. Their purpose is to soften real-time guarantees leading to an overestimated use of processors, or as we say, provide a pessimistic analysis [START_REF] Luis Díaz | Pessimism in the stochastic analysis of real-time systems: Concept and applications[END_REF]. Associating probabilities to worst-case values and to other possible cases may reduce the pessimism of a real-time analysis. Often, execution times have associated probabilities, but inter-arrival times may also be considered as probabilistic. In fact, many real-time systems have arrival times subject to noise and other interference from either operating systems etc. which can be modeled with probability distributions [START_REF] Wehaiba | Work in progress: Kdbench -towards open source benchmarks for measurement-based multicore wcet estimators[END_REF]. We call probabilistic real-time systems the real-time systems which have probabilistic parameters such as execution times and/or inter-arrival times.

Utilization-based schedulability conditions for single-core fixed-priority preemptive scheduling policies are well known [START_REF] Robert I Davis | A review of priority assignment in real-time systems[END_REF]. The seminal work of Liu and Layland [START_REF] Laung | Scheduling algorithms for multiprogramming in a hard-real-time environment[END_REF] introduces a sufficient condition for the feasibility of a real-time system using its maximal utilization. Nevertheless, a real-time system not satisfying this sufficient condition remains schedulable with a given probability (see [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF] in Section 3). Moreover, while probabilistic methods have been focused towards fitting in this sufficient condition by providing less pessimistic analyses, their domain of feasibility needs to be defined as well. In this paper, we build necessary feasibility conditions for fixed-priority scheduling policies based on the mean utilization of the real-time system. We demonstrate that a mean utilization smaller than 1 is mandatory for response times to be finite, c.f. Propositions 4 and 5. We call systems with a mean utilization smaller than 1, stable real-time systems.

Response times depend on properties of real-time systems such as the scheduling policy, the preemptiveness, etc.. Our motivation is to exploit those properties leading response times to the domain of a certain probability distribution. In this work, we emphasize the inverse Gaussian distribution as the appropriate distribution for an approximation of response times in the context of fixed-priority scheduling policies, using asymptotic results of queueing theory, c.f. Propositions 8. We propose two different approximations, a worst-case approximation before the system is in its steady-state (c.f. Proposition 6, see Definition 10) and another one when the system is steady, c.f. Propositions 1 and 7. We prove in Proposition 2 that a real-time system is in its steady-state at the first idle time.

The contributions of this paper are:

• a necessary condition for the feasibility of probabilistic real-time systems with implicit deadlines under preemptive fixed-priority scheduling policies,

• an approximation of the demand of probabilistic real-time systems,

• a first proof that the Liu and Layland theorem [START_REF] Laung | Scheduling algorithms for multiprogramming in a hard-real-time environment[END_REF] stands only for stable probabilistic real-time systems,

• a first proof that a real-time system is steady after the first idle time and an analytical formulation of its distribution,

• an analytical formulation of response time distributions of tasks before and after the first idle time,

• numerical evidence and comparison against simulated results and other existing estimations.

The paper is organized as follows. We present in Section 2 relevant existing results and in Section 3 the main definitions of the considered model. We introduce in Section 4 the Loynes theorem and the heavy-traffic theorem which are the central theorems used this work. These theorems provide analytical expressions of the steadiness of the system, through the stochastic analysis of backlogs. In Section 5 we provide the time demand analysis associated to the probabilistic real-time systems. Moreover, we demonstrate the importance of the first idle time within any feasibility analysis. We detail the experimental results associated to the proposed contribution in Section 6 and we discuss future work opened by our current contribution in Section 7.

Related work

Probabilistic methods for the analysis of response times have many applications in real-time systems [START_REF] Robert | A survey of probabilistic timing analysis techniques for real-time systems[END_REF]. Two main directions have been explored: static methods for the exact computation and approximation of response time distributions [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF][START_REF] Kim | An exact stochastic analysis of priority-driven periodic real-time systems and its approximations[END_REF][START_REF] Manolache | Memory and time-efficient schedulability analysis of task sets with stochastic execution time[END_REF][START_REF] Maxim | Response time analysis for fixedpriority tasks with multiple probabilistic parameters[END_REF] a priori, and the measurement-based application of the extreme value theory (EVT) [START_REF] Liu | An evt-based worst-case response time analysis of complex real-time systems[END_REF][START_REF] Lu | A statistical response-time analysis of real-time embedded systems[END_REF] approximating the distribution of the maximum values of response times a posteriori. Often, probabilities are considered for execution times and few papers consider probabilistic inter-arrival times and deadlines [START_REF] Gaujal | Dynamic speed scaling minimizing expected energy consumption for real-time tasks[END_REF][START_REF] John P Lehoczky | Real-time queueing theory[END_REF][START_REF] Maxim | Response time analysis for fixedpriority tasks with multiple probabilistic parameters[END_REF]. Moreover, the method introduced in [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF][START_REF] Kim | An exact stochastic analysis of priority-driven periodic real-time systems and its approximations[END_REF] requires a large amount of convolutions which have a high space and time complexity, and the analysis provided by Lehoczky [START_REF] John P Lehoczky | Real-time queueing theory[END_REF] is not suited to express response time distributions.

The contributions of this paper are based on queueing theory results [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF][START_REF] Chen | Fundamentals of queueing networks: Performance, asymptotics, and optimization[END_REF][START_REF] Sethuraman | Optimal stochastic scheduling in multiclass parallel queues[END_REF][START_REF] Panayotis | Optimal routing and scheduling of customers with deadlines[END_REF]. To the best of our knowledge, no result from the queueing theory is focused on general execution times and inter-arrival times, multi-class clients (i.e. different tasks) and the quantization of deadline misses of such systems. The only results based on queueing theory for real-time systems have been published within the thread of papers related to [START_REF] John P Lehoczky | Real-time queueing theory[END_REF], where the author approximates the number of simultaneously activated jobs by a reflected at the origin Brownian motion. However, the author makes strong hypotheses restricting the model. The exponential distribution of inter-arrivals and execution times suggested in [START_REF] John P Lehoczky | Real-time queueing theory[END_REF] is a strong hypothesis. Furthermore, his model has another important limitation as it considers systems of jobs of only one task, which does not allow a response time analysis.

Model

We consider a probabilistic real-time system Γ composed of a finite number of tasks ordered in the decreasing priority order, i.e. a task τ i has a higher priority than τ j if and only if i < j. The tasks are scheduled on a single-core processor, using a fixed-priority preemptive policy.

Probabilistic real-time tasks

A probabilistic task τ k is characterized by: 1. its execution time C k > 0 , with a distribution function F k , and finite mean

m k ≜ E[C k ] and non-negative variance s 2 k ≜ E[(C k -m k ) 2 ],
2. its inter-arrival time

T k > 0 between two consecutive instances of τ k , with a distribution function G k of rate λ k ≜ 1/E[T k ],
3. its relative implicit deadline D k > 0, with the same distribution function G k as its inter-arrival time.

Without loss of generality, we consider the functions F k and G k discrete. Also note that all the parameters written in capital letters are probabilistic.

The l-th instance of the task τ k is called a job and we denote it τ k,l . Its execution time is C k,l , with distribution function F k . T k,l is the inter-arrival time between the arrival times of the jobs τ k,l-1 and τ k,l , with distribution function G k , i.e. identically distributed as T k . We suppose all the rates of arrival λ k distinct.

Due to the fixed-priority policy, the response times of the task τ k depend only those of higher priority tasks. We call job of level k any job of a task of higher or equal priority than τ k , i.e. any job τ i,j , Let Γ be a real-time system with total mean utilization ū. Γ is said stable if and only if ū < 1.

1 ≤ i ≤ k, j ∈ N. Let u i ≜ E[C i ]/E[T i ] = λ i m i be the mean utilization of the task τ i , ūk ≜ k i=1 u i be the k-level mean utilization of Γ , ū ≜ τ k ∈Γ u k the total
The response time R k,l of a job τ k,l is the elapsed time between its activation (or arrival) A k,l ≜ l j=1 T k,j and the end of its execution. We consider implicit deadlines, i.e. the absolute deadline of the job τ k,l is A k,l+1 and A k,l+1 -A k,l = T k,l+1 ∼ D k is its relative deadline with distribution function G k . For the sake of comprehension, we refer to the relative deadline of τ k with the variable D k and to its inter-arrival time T k , even though they have the same distribution. At each arrival of jobs of a same task, previous jobs miss their deadline and are discarded.

Remark 1. Such discarding policy is the reason why response times of jobs of a same task are independent: response times are functions of their associated execution times and only the higher priority jobs execution times. As described in [START_REF] Huang | Control of patient flow in emergency departments, or multiclass queues with deadlines and feedback[END_REF], discarding jobs makes their backlog disappear from the computation of the response time of the next job of the same task. However, in the response time analysis in Section 5 we build response time approximations regardless of this discarding policy which makes them upper-bounds of response times, see [START_REF] Janssen | Renewal theory[END_REF].

Example 1 (Utilization and probabilistic deadlines). Let us consider the task set {τ 1 , τ 2 }, with [START_REF] Wehaiba | Work in progress: Kdbench -towards open source benchmarks for measurement-based multicore wcet estimators[END_REF]. Furthermore, because deadlines are implicit, if T 2,2 = 3.1, the deadline of the job τ 2,1 is 3.1, if T 2,2 = 4 its deadline is 4. The deadline of any job τ 1,j is 2 because T 1,j is always 2 for all j ∈ N.

P(C 1 = 1) = 1, P(T 1 = 2) = 1, P(C 2 = 1) = 1/2, P(C 2 = 2) = 1/2, P(T 2 = 3.1) = 1/2, P(T 2 = 4) = 1/2. Then ū = 1/2 + ((1 + 2)/2)/((3.1 + 4/2) ≈ 0.92, u max = 1/2 + 2/3.1 ≈ 1.
A job τ k,l respects its deadline if its response time R k,l is smaller or equal to its deadline, i.e. R k,l ≤ D k . In order to be schedulable, the worst-case deadline failure probability [START_REF] Navet | Worst-case deadline failure probability in real-time applications distributed over controller area network[END_REF][START_REF] Georg Von Der Brüggen | Efficiently approximating the worst-case deadline failure probability under edf[END_REF] of a task τ k should be smaller than its permitted failure rate α k ∈ (0, 1) [START_REF] Guo | Edf schedulability analysis on mixed-criticality systems with permitted failure probability[END_REF], i.e.

P sup l∈N R k,l > D k = P sup l∈N R k,l > t dG k (t) ≤ α k (1)
where α k is a parameter given a priori.

Definition 2. Let Γ be a real-time system. Γ is said feasible in the fixedpriority domain if and only if there exists at least one fixed-priority preemptive scheduling policy such that all tasks τ k ∈ Γ satisfy (1).

Remark 2. For all tasks τ k ∈ Γ and all jobs τ k,l , l ∈ N the execution times C k,l are independent. We discuss in Section 7 how the model may be extended in the dependent case, i.e. when the execution times C k,l , τ k ∈ Γ, l ∈ N have a dependence structure. The inter-arrival times T k,l , l ∈ N of a given task τ k are independent, however inter-arrival times of two different tasks do not need to be independent for the provided results.

For a given job τ k,l , its execution time C k,l and the inter-arrival time T k,l are not necessarily independent, meaning that the analysis we provide holds if execution times are correlated to the time separating a job to the previous instance of the same task.

In the rest of this paper, we refer to probabilistic real-time systems as simply real-time systems when there is no ambiguity.

Elements of probability

Let us define some elements of probability that we use in this work. Let X and Y be two non-negative variables.

Definition 3 (Distribution function, probability function and distribution).

The distribution function of X is the function defined by F X (x) = P(X ≤ x). The probability function of X is the function p X (x) = (dF X /dx)(x). We denote the distribution of X by dF X (x) = p X (x)dx. Definition 4 (Identically distributed). X and Y are identically distributed and we write X ∼ Y , or X ∼ dF Y , if and only if F X (x) = F Y (x) for all x > 0.

Definition 5 (Independence). X and Y are independent if and only if P(X ≤ x, Y ≤ y) = F X (x) • F Y (y) for all x, y > 0.

Definition 6 (Conditional distribution). The conditional distribution function of X given Y = y is written

P(X ≤ x | Y = y) = 1 p Y (y) • d dy P(X ≤ x, Y ≤ y)
for all x, y > 0, and is such that

P(X ≤ x) = P(X ≤ x|Y = y)dF Y (y) (2)
for all x > 0.

Definition 7 (Convolutions). Let X and Y be independent and Z = X + Y .

Then the distribution function of Z is the convolution

F Z = F X * F Y defined by F Z (z) = F X (z -y)dF Y (y) = F Y (z -x)dF X (x)
Definition 8 (Brownian motions). A standard Brownian motion is a process B such that B(0) = 0 and with independent increments B(t) -B(s) of distribution function Φ •/ √ t -s for any 0 ≤ s < t, where

dΦ(x) = exp -x 2 /2 / √ 2πdx
is the standard normal distribution.

A Brownian motion of drift u and deviation v > 0 is a process W such that there exists a standard Brownian motion B such that

W (t) = ut + vB(t) (3)
Definition 9 (Inverse Gaussian distribution). Let θ > 0 and γ > 0. The inverse Gaussian distribution of parameters (θ, γ) has a probability function ψ defined by

ψ (t; θ, γ) = γ 2πt 3 exp - γ(t -θ) 2 2θ 2 t (4)
Remark 3. Note that deterministic1 real-time systems form a subclass of probabilistic real-time systems, as the distributions of the parameters of Γ can describe a unique value if needed. Hence, we propose results covering also deterministic parameters, would they be execution times, inter-arrival times or deadlines. For instance, if the execution time C k of a task τ k ∈ Γ is deterministic and equal to 5, then its distribution function is

F k (x) = 1 [5,∞) (x)
, where

1 A (x) = 1 if x ∈ A 0 otherwise

Backlog stochastic process

Let us define the following stochastic processes:

1. N k (t) ≜ ∞ l=1 1 [0,t] (A k,l ) as the number of jobs of τ k released before t > 0, of mean E[N k (t)] = λ k t [1, (1.1.22)]. N k is right-continuous.
2. the k-level demand W k (t) as the accumulation of the execution times required by jobs of priority higher or equal than τ k , regardless of potential deadline misses, released before the instant t to complete,

W k (t) ≜ k i=1 Ni(t) j=1 C i,j of mean E[W k (t)] = k i=1 E[N i (t)]E[C i ] = ūk t, see [14, (6.53)]. W k is right-continuous.
3. the k-level backlog β k (t) as the remaining demand at t ≥ 0, defined by β k (0) ≥ 0 and

β k (t) ≜ β k (0) + W k (t-) - t 0 1 {β k (s)>0} ds (5)
where W k (t-) = lim ε→0,ε>0 W k (t -ε) is the left limit of W k at the instant t, and

t 0 1 {β k (s)>0} ds is the total busy time of level k [18, p. 2] before t > 0.
Example 2 (Backlog). Consider the task set of Example 1. Suppose both tasks τ 1 and τ 2 are activated at time t = 0 and T

2,2 = 3.1, C 2,1 = 2, C 2,2 = 1. Then, β 2 (0) = 0, β 2 (2) = β 2 (0) + W 2 (2-) - 2 0 1 {β2(s)>0} ds = 0 + (1 + 2) -2 = 1 β 2 (3) = β 2 (0) + W 2 (3-) - 3 0 1 {β2(s)>0} ds = 0 + (1 + 2) -3 = 0 β 2 (3.1) = β 2 (0) + W 2 (3.1-) - 3.1 0 1 {β2(s)>0} ds = 0 + (1 + 2) -3 = 0 β 2 (3.2) = β 2 (0) + W 2 (3.2-) - 3.2 0 1 {β2(s)>0} ds = 0 + (1 + 2 + 1) -3.1 = 0.9
At the instant t = 3.2, the backlog of level 2 is 0.9.

The process β k describes the remaining demand without considering deadline misses, i.e. while jobs are discarded their demand remains in the backlog analysis. Thus, the process β k is an upper-bound of the blocking time, see [START_REF] Janssen | Renewal theory[END_REF], i.e. the response time analysis that we provide in Section 5 is pessimistic as defined in [START_REF] Luis Díaz | Pessimism in the stochastic analysis of real-time systems: Concept and applications[END_REF].

In stable real-time systems, there exist two main states of the backlog [8, Section 4.2].

Definition 10 (Transient and steady states). For a stable real-time system Γ , the steady-state backlog is defined by

βk ≜ lim t→∞ β k (t)
We say that Γ is steady at the instant t > 0 when its backlog is in steady-state, i.e. β k (t) ∼ βk , and transient before β k reaches its steady-state. : Level 3 demand of 1 000 instances of the DK model system, the heavy-traffic demand process and the classical deterministic worst-case analysis (WCET) considering only the maximal execution time for each task, for ū = 0.838 and u max = 1.208 and hyper-period T , for Γ defined in Table 1.

In the following sections, we define the domain of this steadiness, provide an expression of the instant when the system goes from transient to steady, and prove that this time instant exists and is finite under some conditions.

Heavy-traffic theory

In this section, we present an analytical approximation of the demand of probabilistic real-time systems. Let us denote π k (x) ≜ P βk ≤ x the distribution function of the steady-state backlog of level k. We use background results of queueing theory presented in 4.1 and 4.2, provide the exact formulation of the steady-state backlog distribution π k in 4.3 and illustrate this result in the deterministic inter-arrival case already studied by Diaz et al. [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF] in Section 4.4.

The Loynes theorem

As the demand and backlog processes W k and β k are well studied in queueing theory [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF][START_REF] Chen | Fundamentals of queueing networks: Performance, asymptotics, and optimization[END_REF][START_REF] Jeanblanc | Mathematical methods for financial markets[END_REF], we provide the formula of the steady-state of the system, by adapting the Loynes theorem

[1, Theorem 2.1.1]. Theorem 1 (Loynes theorem). Let A k l = inf{t > 0 : k i=1 N i (t) = l}
C I k l ,l -(A k l+1 -A k l ) + (6)
where x + = max(0, x), c.f. [ Our goal is to find the distribution of the steady-state backlog βk . However, given the generality of this model, we cannot provide an exact description of the process β k . This is why we use the heavy-traffic assumptions to find an exact formula for the distribution of the steady-state backlog βk when the system utilization gets close to 1, c.f. Figure 2, i.e. we build a process βk such that its steady-state is also βk .

The heavy-traffic theorem

A first step in the approximation of the backlog process β k is the approximation of the demand process W k . Theorem 2 (Heavy traffic theorem, Section 6.5, [START_REF] Chen | Fundamentals of queueing networks: Performance, asymptotics, and optimization[END_REF] ). The level k heavy-traffic demand process

Ŵk (t) ≜ lim n→∞ W k (nt)/ √ n is a Brownian motion of drift ūk and deviation v k , where v 2 k ≜ k i=1 λ i • s 2 i . See Figure 1.
In [START_REF] John P Lehoczky | Real-time queueing theory[END_REF], the author uses the heavy-traffic approximation providing the distribution of the lateness of jobs in a system with exponential inter-arrival times. To illustrate the heavy-traffic approximation, one can think of water continuously flowing into a sink at a rate λ i and the execution times as the rate 1/m i the water leaves the sink. It is usually understood as true when the system is at full processor utilization, because the theorems of heavy-traffic theory are exact when ū → 1. However we use it as a way to build an upper-bound of the backlog process. Indeed, the heavy-traffic assumption should be seen as a bound, or more specifically, a way to suppose that the system utilization is at its maximum (i.e. ū = 1), providing upper-bounds that are exact when the processor utilization at 100%. Figure 3b illustrates this upper-bound becoming exact in Figure 3c and Figure 3d.

In this section, we use the heavy-traffic assumption to find the steady-state backlog βk as an approximation upper-bounding the blocking time of the system in its steady-state, see [START_REF] Janssen | Renewal theory[END_REF]. The approximation in Theorem 2 leads to the standard Brownian motions which is continuous. It means that instead of looking at the demand for a large amount of time, we consider a re-scaled version of the demand in order to build a good approximation. Theorem 2 can be written as

Ŵk (t) = ūk t + v k B(t) (7)
where B is a standard Brownian motion, as defined in (3).

Steady-state backlog distribution

In order to find the steady-state backlog βk , we work with the heavy-traffic demand Ŵk . We define the heavy-traffic backlog process by

βk (t) = β k (0) + Ŵk (t) - t 0 1 { βk (s)>0} ds (8)
This heavy-traffic backlog process βk is a Brownian motion reflected at the origin of drift ūk -1 and deviation v k , and its limit is also βk , c.f. [4, Section 6.5]. Its stationary distribution π k being the same for the initial backlog process β k , i.e. lim t→∞ β k (t) = βk = lim t→∞ βk (t) we get its explicit expression in the following proposition. Proposition 1. Let ε > 0 and ū = 1 -ε. Consider the parameter

η k ≜ 2 (1 -u k ) / λ k (γ 2 a,k + γ 2 e,k ) (9)
where γ a,k and γ e,k denote the coefficients of variation2 of respectively interarrival and execution times of τ k . Then, when ε → 0, the distribution function π k of the steady-state backlog βk is defined for x > 0 by

π k (x) = k i=1 (1 -exp (-η i x)) k j=1 j̸ =i (1 -η i /η j ) -1 (10)
Proof. From [4, Section 6.5, p. 144], we know that when ε → 0, βkβk-1 is an exponential variable of parameter η k . As execution times are independent, βiβi-1 and βjβj-1 are independent for any i ̸ = j. Furthermore, βk = β1 + k i=2 βiβi-1 , thus we know from [START_REF] Sheldon | Probability models for computer science[END_REF] that the distribution of βk is

dπ k (x) = k i=1 η i k i=1 exp (-η i x) k j=1,j̸ =i (η j -η i )
dx for x > 0, which is the convolution of k exponential distributions of parameters η 1 , . . . , η k . See Figure 3 for an illustration of the probability function of π k .

|β k (s) -βk (s)| = O (t log log t) 1/4 (log t) 1/2 (11)
Note that [START_REF] Friebe | Identification and validation of markov models with continuous emission distributions for execution times[END_REF] applies to the general case, hence the deterministic case, which means that the approximation of the stationary distribution in the periodic case is given in an explicit formula, without solving any system of equations neither computing a large number of convolutions.

The Diaz and Kim model: periodic workloads

For the case where the tasks of the real-time system Γ are periodic and have deterministic deadlines, i.e.

T k = λ -1 k ∈ N and G k (x) = 1 [λ -1 k ,∞) (x)
for all τ k ∈ Γ and x > 0, also known as the Diaz and Kim (DK) model [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF][START_REF] Luis Díaz | Pessimism in the stochastic analysis of real-time systems: Concept and applications[END_REF][START_REF] Kim | An exact stochastic analysis of priority-driven periodic real-time systems and its approximations[END_REF], the authors approximate the distribution function π k , resolving linear system equations and compute response times distributions with the help of convolutions.

The fastest computational complexity of convolutions is O(N log N ) when N is the number of values that a probabilistic variable can take. Computing the exact values of π k quickly becomes an expensive operation when the number of tasks or the number of possible execution times gets larger, even with methods providing tools to soften those computations like Markovic et al. [START_REF] Marković | On the Convolution Efficiency for Probabilistic Analysis of Real-Time Systems[END_REF], Milutinovic et al. [START_REF] Milutinovic | Speeding up static probabilistic timing analysis[END_REF] or Maxim et al. [START_REF] Maxim | Re-sampling for statistical timing analysis of real-time systems[END_REF] for example. Moreover, the computation of response times has the same problem, as the number of possible values of βk quickly becomes large. With [START_REF] Robert I Davis | A review of priority assignment in real-time systems[END_REF] we have an explicit formula of the distribution of βk , with [START_REF] Gaujal | Dynamic speed scaling minimizing expected energy consumption for real-time tasks[END_REF] we know the error of the heavy-traffic approximation, and with Theorem 1 we have an analytical expression of the backlog in the deterministic case.

Indeed, let Tk = lcm λ -1 1 , . . . , λ -1 k be the hyper-period of level k of Γ . In the DK model [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF][START_REF] Kim | An exact stochastic analysis of priority-driven periodic real-time systems and its approximations[END_REF], the authors consider the k-level backlog β k (t Tk ), t ∈ N, i.e. the remaining demand of level k at the beginning of the t-th hyper-period. Diaz et al. [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF] have proven that the sequence β k t Tk t∈N is a stationary Markov chain when ūk < 1. The sequence β k (t Tk ) t∈N is defined by β k (0) ≥ 0 and

β k (t + 1) Tk = β k (t Tk ) + W k ( Tk ) -Tk + for t ∈ N.
In order to apply [START_REF] Robert I Davis | A review of priority assignment in real-time systems[END_REF] to the DK model, we consider a system with only one task with an execution time equal to the total demand of one hyper-period W k ( Tk ). Indeed, the demand of one hyper-period is repeated at each new hyperperiod [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF]. Let W k,t , t ∈ N be independent and identically distributed variables with the same distribution than W k ( Tk ). Then from [START_REF] Robert I Davis | A review of priority assignment in real-time systems[END_REF] we have

βk ∼ sup n≥1 n t=1 W k,t -n Tk + (12)
The representation in ( 12) is an efficient method to approximate the stationary distribution π k of the backlog process β k . Indeed, let us take an integer n > 0, and generate a sample (W k,t ) n t=1 independent and identically distributed sequence with the distribution of W k Tk . Eq. ( 12) provides the variable of distribution π k found in Diaz et al. [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF] and Kim et al. [START_REF] Kim | An exact stochastic analysis of priority-driven periodic real-time systems and its approximations[END_REF]. It also means that the variable max 1≤j≤n j t=1 W k,t -Tk + is an approximation of βk when n is large enough. This method is not expensive in complexity as it requires only to build the distribution function of W k Tk once. Remark 4. In the DK model, the coefficient of variation γ a,k is 0, thus

η k = 2m 2 k 1 -u k λ k s 2 k

Time demand analysis

In this section, we express the response times distribution of a task τ k using the heavy-traffic analysis of Section 4. The Markovian property of Brownian motions allows to approximate the distribution of any response time R k,l in terms of the backlog βk-1 (A k,l ) and the first response time R k,1 , thus response times will be conditioned to backlogs and execution times, and represented as idle times following the inverse Gaussian distribution. In a second part, we provide an analytical expression of the worst-case heavy-traffic response time distribution, and in a third part we do the same for the steady-state heavy-traffic response time distribution. Finally, we explain how to simulate heavy-traffic response times of a task τ k from the distribution functions F k and G k .

The response time R k,l of a job τ k,l is the smallest instant after its arrival time A k,l smaller than the time it is given to run the level k -1 demand and its execution time, i.e. the smallest t > 0 such that

B k-1 (A k,l ) + C k,l + W k-1 (t + A k,l ) -W k-1 (A k,l ) ≤ t ( 13
)
where

B k (t) = k i=1 min (β i (t) -β i-1 (t), C i (t))
is the blocking time of level k at the instant t > 0, β 0 (t) = 0, and C i (t) is the execution time of the most recent job of τ i released before t, for all t > 0. This relation holds for any fixed-priority preemptive single-core model with the discarding policy discussed in Section 2. By definition, the blocking time is such that

B k (t) ≤ min β k (t), k i=1 C i (t) (14)
for all t > 0, which makes our response time analysis build upper-bounds of response times. Furthermore, the heavy-traffic demand Ŵk-1 is a Brownian motion, thus Ŵk-1 (t + A k,l ) -Ŵk-1 (A k,l ) ∼ Ŵk-1 (t), and it is continuous. Then we define the heavy-traffic response time of a job τ k,l as

Rk,l ∼ inf t>0 βk-1 (A k,l ) + C k,l + Ŵk-1 (t) = t (15)
In the following sections we consider response times as idle times properly conditioned, and from this analysis provide the worst-case response time distribution and the steady-state response time distribution.

First idle time

In the following, we consider the conditional probability that the system starts with a level k -1 backlog x ≥ 0 and the job τ k,1 ,

P x k (•) = P(• | β k-1 (0) = x, A k,1 = 0)
Whenever we need to suppose β k-1 (0) = x, A k,1 = 0, we say that we work under the probability P x k . As the steady state of backlogs repeats itself infinitely, the analysis of response time stands for all future response times. In the transient state though, each response time has its own distribution function. Hence, in Section 5.3 we A question rises at this point: when does the system get steady ? Let us look closely at the dynamics of the heavy-traffic backlog process. Let β(t) ≜ max τ k ∈Γ βk (t), Ŵ (t) ≜ max τ k ∈Γ Ŵk (t) and v ≜ max τ k ∈Γ v k . As we work with fixed-priority scheduling policies, β, Ŵ and v refer to the lowest priority level.

Proposition 2. Let ū < 1, x > 0 and β(0) = x. The steady-state backlog is reached at the first idle time

I x ≜ inf t > 0 : β(t) = 0
Proof. For any x > 0, at the instant I x , the next idle time has the distribution of I 0 , so is the third, the fourth, etc. This is sufficient to say that the system is steady for all t > I x . Moreover, one can see in [4, Theorem 6.8] that when ū < 1 the steady-state backlog β is 0.

Proposition 2 is coherent with the proven results in the deterministic case [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF] and the simulations shown in Figure 2. Let us now refine the analysis for each task τ k ∈ Γ . Let I x k be the first idle time of level k, i.e. when the initial backlog of level k is equal to x > 0.

I x k = inf t > 0 : βk (t) = 0 (16)
Lemma 1. The distribution of I x k is an inverse Gaussian distribution with probability function

ψ x k (t) ≜ ψ t; x 1 -ūk , x 2 v 2 k , t > 0
Proof. The idle time I x k defined in ( 16) is a quantity called first-passage time of a Brownian motion [28, Eq. ( 28)]. Until the first idle time, we know that βk (t) = x + Ŵk (t) -t, because 1 β k (s)>0 = 1 for 0 ≤ s < t ≤ I x k . Then from (7) we have

I x k ∼ inf {t > 0 : B(t) = t(1 -ūk )/v k -x/v k } where B is a standard Brownian motion. When β k (0) = x, the distribution of I x
k is an inverse Gaussian distribution of mean x/(1 -ūk ) and shape x 2 /v 2 k . See Figure 5.1 and [15, p. 146], for more details.

The first idle time is also a probabilistic variable, but we can build an instant t idle (ε, x) such that for any ε > 0 and all t > t idle (ε, x), the probability that the first idle time of the lowest priority level is greater than t is less than ε, i.e. P (I x > t) ≤ ε. We call first ε-idle time the instant t idle (ε, x) from which we can guarantee the system is steady with probability 1 -ε and initial backlog x > 0.

We know that β is a Brownian motion reflected at the origin. This means that until it reaches zero, i.e. an idle time, it is a simple Brownian motion. In fact, between two consecutive idle times the backlog process has the same dynamic. The following proposition bounds with a given probability the first idle time and the instant the system gets steady. Proposition 3. Let ε ∈ (0, 1) and x ≥ 0. If ū < 1, the system of initial demand x is steady with probability at least 1 -ε at the instant

t idle (ε, x) =   q 1-ε + q 2 1-ε + 4 1-ū v x 2 1-ū v   2 ( 17
)
where

q 1-ε = Φ -1 (1-ε) is the (1-ε)-quantile of a standard normal distribution.
Proof. Let us consider the lowest priority first idle time I x and remark that

P(I x > t) = P inf s∈[0,t] Ŵ (s) -s > -x (18)
≤ P Ŵ (t) > t -x = P B(t) > t 1 - ū v - x v = 1 -Φ (1 -ū)t -x v √ t
Hence, let us define t idle (ε, x) such that

1 -Φ (1 -ū)t idle (ε, x) -x v t idle (ε, x) ≤ ε
which implies that for t > t idle (ε, x) we should have

t > ūt + Φ -1 (1 -ε)v √ t + vx (19)
see Figure 5.1. This is a second order polynomial equation, which admits no solution when ū > 1, and, when ū ≤ 1 and ε ∈ (0, 1), the smallest solution of ( 19) is [START_REF] Kim | An exact stochastic analysis of priority-driven periodic real-time systems and its approximations[END_REF].

Note that for all x ≥ 0, lim ε→0 t idle (ε, x) = ∞, lim ε→1 t idle (ε, x) = 0, and

• ū < 1 =⇒ t idle (ε, x) < ∞, • ū = 1 =⇒ t idle (ε, x) = ∞, • ū > 1 =⇒ t idle (ε, x) is not defined.
Remark 5. We provide the analysis for the lowest priority level, but in fact each priority level has its own first ε-idle time.

Finally we consider the maximum first ε-idle time t max idle (ε) by considering a synchronous activation for the distribution µ max of the initial demand, i.e. β(0) = τ k ∈Γ C k . See Proposition 5 for a more detailed explanation. The maximum first ε-idle time is then

t max idle (ε) = t idle (ε, x)dµ max (x) (20)

Conditioning response times

From [START_REF] Jeanblanc | Mathematical methods for financial markets[END_REF], any job τ k,l has a heavy-traffic response time distribution that can be expressed from the response time Rk,1 ,

P Rk,l > t | βk-1 (A k,l ) = x = P x k Rk,1 > t (21)
Let us condition this probability for specific values of execution times.

Lemma 2. Under P x k the heavy-traffic response time Rk,1 given C k,1 = y is identically distributed as I x+y k-1 . Proof. As stated in [START_REF] Liu | An evt-based worst-case response time analysis of complex real-time systems[END_REF], the proper conditioning on backlogs provides the distribution of the response time of τ k,1 . Furthermore, when the backlog βk-1 (A k,l ) = x and the execution time C k,l = y, the response time Rk,l is the time it takes for all level k -1 jobs to finish plus the time it takes for level k -1 to stay idle for x + y instants. It means that we can artificially set the initial backlog to x + y and look at I x+y k-1 , the first idle time of level k -1, as represented in Figure 6. In other words,

P x k Rk,1 > t | C k,1 = y = P inf s∈[0,t] Ŵk-1 (s) -s > -x -y = P I x+y k-1 > t
which is sufficient to conclude.

Example 3 (Response times as idle times). Consider a task set {τ 1 , τ 2 }, and that τ 2 is activated at t = 0, i.e. A 2,1 = 0, and

C 2,1 = y = 3. Let P(C 1 = 1) = 1/2, P(C 1 = 2) = 1/2, P(T 1 = 2) = 1/2, P(T 1 = 4) = 1/2, and suppose β 1 (0) = x + y = 8, C 1,1 = 1, C 1,2 = 2, C 1,3 = 2, and A 1,1 = 2, A 1,2 = 6, A 1,3 = 10 and A 1,4 = 14.
Then the response time R 2,1 is the first instant t when β 1 (t) = 0: Hence R 2,1 = 13. This is illustrated in Figure 6. To build all possible values of R 2,1 , one must do the same for all combinations of possible values of (A

β 1 (1) = β 1 (0) = x + y β 1 (3) = β 1 (0) + C 1,1 -3 = x + y -2 β 1 (7) = β 1 (0) + C 1,1 + C 1,2 -7 = x + y -4 β 1 (11) = β 1 (0) + C 1,1 + C 1,2 + C 1,3 -11 = x + y -6 β 1 (13) = β 1 (0) + C 1,1 + C 1,2 + C 1,3 -13 = x + y -8 = 0 t β 1 (t) x + y A 1,1 A 1,2 C 1,2 A 1,3 R 2,1 = I x+y 1
. Lemma 2 provides an analytical expression of this result for the heavy-traffic response time R2,1 .

From ( 15) and Lemma 2 we establish a necessary condition of the feasibility of Γ . Proposition 4. A non-stable real-time system with implicit deadlines and independent execution times is not feasible under a fixed-priority scheduling policy.

Proof. Let τ k ∈ Γ be such that ūk > 1. The Loynes theorem 1 states that there is a finite number of idle times of level k, which implies with (15) that heavy-traffic response times get infinite at some point, i.e. for all t > 0, P sup l∈N Rk,l > t = 1. In other words, there is no permitted failure rate α k ∈ (0, 1) such that τ k is schedulable as defined in [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF]. As we consider a fixedpriority scheduling policy, then the lowest level backlog is larger than all k-level backlogs β k (t) at any time t > 0.

With the same reasoning we get the explicit formula of the distribution function of response times.

Let us denote the conditional probability that the distribution of the initial backlog β k-1 (0) is the probability dµ and the first job released is τ k,1 ,

P µ k (•) ≜ P(• | β k-1 (0) ∼ dµ, A k,1 = 0) = P x k (•)dµ(x)
Lemma 3. Under P µ k , the probability function of the heavy-traffic response time Rk,1 is

h µ k (t) = ψ x k-1 (t)d(µ * F k )(x)
Proof. We conclude with Lemmas 1 and 2 and the definition of conditional probabilities, c.f. (2).

In Sections 5.3 and 5.4, we prove that the proper initialization of the system puts the system in two specific cases: the worst-case and the steady-state.

Worst-case response time

Before the system reaches its steady-state, we say it is transient. In that case we cannot provide the exact distribution of the backlog in an analytical formulation. However, we can bound it by using the worst-case backlog.

Let

µ k = F 1 * • • • * F k be the distribution function of the demand of level k of a synchronous activation k i=1 C i . Proposition 5. If ūk-1 < 1,
the worst-case heavy-traffic response time distribution of τ k is defined by

P sup l∈N Rk,l > t = P µ k-1 k Rk,1 > t (22)
Proof. First, let us consider ūk-1 < 1, as stated in Theorem 1 the backlog process converges to βk which is finite. Jobs are discarded if they miss their deadlines, and as we consider implicit deadlines, there can be at most one job per task activated simultaneously. This leads into considering

β k-1 (0) = k-1
i=1 C i as the maximum backlog of level k -1, c.f. [START_REF] Janssen | Renewal theory[END_REF], and use the property stated in [START_REF] Liu | An evt-based worst-case response time analysis of complex real-time systems[END_REF].

When ūk-1 = 1, idle times of level k may or may not be finite, thus we cannot conclude anything on the distribution of response times.

When ūk-1 > 1, the largest response time does not come from the synchronous activation and is in fact ∞. As we have already demonstrated in the proof of Proposition 4, response times of the task τ k increase to ∞, due to the absence of idle times of level k -1. Then we conclude that the worst-case heavy-traffic response time is ∞.

See Figure 3.

The worst-case response time being identified as the synchronous activation case, we provide an explicit expression its distribution function. Proposition 6. Let τ k ∈ Γ and ūk < 1. The probability function of the worstcase heavy-traffic response time of τ k is

h max k (t) ≜ ∞ 0 ψ x k-1 (t)dµ k (x) (23)
Proof. From Proposition 5 we know that when ūk < 1, the worst-case response time of τ k ∈ Γ is Rk,1 under P

µ k-1 k
. Furthermore, from Lemma 3 we have

h max k (t) = ψ z k-1 (t) d(µ k-1 * F k )(z).
In addition, we have the relation µ k-1 * F k = µ k .

Steady-state response time

The backlog process of level k being stationary and with a stationary distribution π k , in the steady-state [START_REF] Liu | An evt-based worst-case response time analysis of complex real-time systems[END_REF] becomes

P π k-1 k Rk,l > t = P π k-1 k Rk,1 > t (24)
for any l ∈ N. Then if it holds for any l ∈ N, it holds for all the response times after the convergence of the backlog process. This is why in the steady-state, the distribution of heavy-traffic response times is unique. Consider the steady-state response time

Rk ≜ inf t > 0 : βk-1 + C k + Ŵk-1 (t) = t (25)
As in Section 5.3, we get the distribution function of steady-state response times.

Proposition 7. Let τ k ∈ Γ , ūk < 1, hk be the probability function of the steady-state response time of τ k , and for i ∈ {1, . . . , k} let η i be as defined in Proposition 1 and

ξ i,k = k j=1,j̸ =i (1 -η i /η j ) -1 . Then for all t > 0, hk (t) = k-1 i=1 ξ i,k-1 ∞ 0 ψ z k-1 (t)dν i,k (z)
where dν i,k (z) = η i z 0 e -ηi(z-x) dF k (x)dz is the distribution of C k + βiβi-1 . Proof. We know from Proposition 1 that the steady-state level (k -1) backlog distribution is π k-1 . We know from Lemma 2 that the stead-state response time Rk is the first idle time with an initial backlog of βk-1 + C k . Finally from Lemma 3, we know that the probability function of Rk is hk ≜ h

π k-1 k
. Then, let us compute the convolution π k-1 * F k which is the distribution function of the sum βk-1 + C k . For any z > 0, we have

d(π k-1 * F k )(z) = dπ k-1
dz * F k (z)dz and

dπ k-1 dz * F k (z)dz = dz ∞ 0 dπ k-1 dz (z -x)dF k (x) =dz z 0 η 1 . . . η k-1 k-1 i=1 exp (-η i (z -x)) k-1 j=1,j̸ =i (η j -η i ) dF k (x) =dz z 0 k-1 i=1 η i e -ηi(z-x) ξ i,k-1 dF k (x) = k-1 i=1 ξ i,k-1 dν i,k (z)
Finally by definition of conditional probabilities, c.f. (2), we have hk (t) =

ψ z k-1 (t) d(π k-1 * F k )(z).
Remark 6. In the periodic case, Diaz [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF] has proven that the system is always steady when u max k ≤ 1, e.g. Figure 3a, i.e. the steady-state backlog is βk = 0, and that is it reached at t = 0.

hk (t) = ψ x k-1 (t)dF k (x) ( 26
)
This steadiness is illustrated in Figure 2.

Simulations

In Lemma 2 we prove that response times R k,l can be simulated from a sample of idle times of level k -1, a sample of execution times and the proper initialization backlog sample, the steady-case or the worst-case.

Generating heavy-traffic response times

The procedure is as follows:

• Generate a backlog b with the distribution function µ k-1 (resp. π k-1 ),

• Generate an execution time c with the distribution function F k ,

• Generate the response time R k = I b+c k-1 , with the probability function ψ b+c k-1 .

See Figure 6 and Figure 7a. Proof. It is a consequence of Lemma 3. Alternatively, we can see that by construction ( 23) and ( 7) come from conditional probabilities. The representation in Lemma 2 indicates that the distribution of response times is conditioned by the values of the backlog and the execution time.

Experimental results

The purpose of this section is to illustrate that the response times generated from Proposition 8 provide a good approximation, by comparing distribution functions of simulations and its associated EVT estimation, and generated heavytraffic response times distributions. Closer are the curves, better is the estimation. Those results are not exhaustive and are an illustration, we do not cover in this work the sensitivity of the model for different values of ū. In order to do so, and to compare the distribution of the generated response times with the ones provided by SimSo, we modified the SimSo tool [START_REF] Chéramy | Simso: A simulation tool to evaluate real-time multiprocessor scheduling algorithms[END_REF] to simulate probabilistic systems and compare the scheduled (SimSo) simulations and heavy-traffic response times simulations. We use the data generated by SimSo to apply EVT on response times (using the Scipy framework3 on Python). Finally we compare our results with SimSo simulations and EVT estimations. For simplicity the periods are deterministic in these simulations.

In order to illustrate the stability described in this paper, we use the task set Γ provided in Table 1, with C k being the set of possible values for the execution time C k and p k = dF k /dx its probability function. The level 2 maximum utilization is smaller than 1, hence backlogs of level 2 converge quickly to 0, see Figure 3a. The level 3 maximum utilization is greater than 1 and the level 3 mean utilization is smaller than 1. Thus, τ 3 is the task of interest, see Figure 2b and Figure 3b. The level 4 mean utilization is close to 1, hence Figure 2c and Figure 3c illustrate the behavior of backlogs when the utilization approaches its phase transition. The level 5 mean utilization is greater than 1, which illustrates the explosion (infinite response times) of the system, see Figure 2d and Figure 3d.

In Figure 2 we observe that the demand W (t) follows the line ūt (its mean). In Figure 2 and3 we see what happens when the system mean utilization is smaller, close and greater than 1 : for smaller values of ū the system stays with zero backlog at some point, see Figure 2b, but for values greater than 1 the system explodes, see Figure 2d and3d. In Figure 2c and 3c we see that even for ū close to 1, the system always admits finite idle times. Most importantly, we see in Figure 2b and 3b that when ū < 1 and u max > 1, the analysis holds and provides quantifiable response times.

In Figure 5.1, we see to what corresponds idle times and ε-idle times graphically. In Figure 5.1 the average idle time corresponds to the point where the line t and ūt meet, and the distribution of the idle times correspond to the frequency the demand process meets the line t for each instant t > 0.

Finally, in Figure 7a, we can see the simulations presented in Section 6.1 and a comparison with response times simulated via SimSo [START_REF] Chéramy | Simso: A simulation tool to evaluate real-time multiprocessor scheduling algorithms[END_REF]. The two ground truth samples are the two subsets SimSo-transient and SimSo-steady, which are composed respectively of simulated response time released before and after the maximum ε-idle time t max idle (ε) = 154 where ε = 10 -6 . The maximum idle time t max idle (ε) is computed via Monte-Carlo approximations with the representation in [START_REF] Laung | Scheduling algorithms for multiprogramming in a hard-real-time environment[END_REF]. The EVT estimation (green curve) from these SimSo simulations and the steady-state and worst-case response times suggested in this paper are compared via their distribution functions in Figure 7b. In this case, where ū is not too close to 1 and u max greater than 1, we can see that the Worst-case response times, the EVT estimation and the steady-state response times are greater (in the stochastic sense [START_REF] Luis Díaz | Pessimism in the stochastic analysis of real-time systems: Concept and applications[END_REF]) than the true response times simulated with SimSo. The heavy-traffic Worst-case response time seems to be an upper-bound in practice, and the Steady-state response time is quite accurate.

Closer is the Steady-state to the SimSo-steady curve, more accurate the approximation is. We can see a big difference between the Worst-case curve and the Steady-state curve. However the proposed analysis does not permit to

Future work

The results of this paper can be generalized if some assumptions of the model are relaxed. Indeed, we have seen that the steady-state only depends on interarrivals only through their mean value in the heavy-traffic analysis. This means that one could suppose inter-arrival independent, also known as Poisson arrivals, and use the wide list of results on Poisson arrivals to extend the results of this paper.

Extension to dynamic-priority policies

The first step into dynamic scheduling, as for example Earliest Deadline First (EDF) is to study the processes β and Ŵ introduced in Section 5.1. Indeed, for fixed-priority policies, those variables are simply the backlog and demand of the lowest priority level. However, for EDF, levels of priority need to be defined not only for tasks but for jobs. In [START_REF] Luis Díaz | Stochastic analysis of periodic real-time systems[END_REF]Section 4.4], authors use the concept of ground jobs which are jobs released at an instant where the system is idle, and as shown in the Loynes theorem, when ū < 1, idle times are finite. This means that an analysis mixing the concept of ground jobs and idle times as defined in this paper can provide an extension of our results for dynamic-priority scheduling.

Extension to dependent execution times

Considering time-driven execution times, i.e. execution times with distribution functions changing over time, is possible. We used drifted Brownian motions to describe the demand. However, we used a model with constant coefficients. A more general class of drifted Brownian motion involve drift and deviation function, depending not only on the values of the demand but also on the the time interval they are being considered. This is a natural extension of the model presented in our work and leads to a more general type of real-time system that needs to be defined. Considering dependence between tasks is also possible. It would lead the analysis to multivariate Brownian motion which requires some new concepts but uses the same arguments we used in this work. The demand should describe not only priority levels but task levels, and cannot be superposed as we have done. Indeed, as a some properties are quite instinctive when independence is assumed (Wald's lemma and the law of large numbers for example), they need some workaround when this assumption is relaxed.

Extension to multi-core processors

The Loynes theorem also exists for the multi-core case [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF], which means that there is a way to apply the results we provide in this paper to a multi-core model, for a global fixed-priority scheduling policy, and a global dynamic-scheduling policy. Furthermore, for independent inter-arrivals, optimal algorithm already exist for priority assignment and processor assignment, c.f. [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF][START_REF] Onno | Workloads and waiting times in single-server systems with multiple customer classes[END_REF][START_REF] Chen | Fundamentals of queueing networks: Performance, asymptotics, and optimization[END_REF][START_REF] Huang | Control of patient flow in emergency departments, or multiclass queues with deadlines and feedback[END_REF][START_REF] Van Mieghem | Dynamic scheduling with convex delay costs: The generalized c-mu rule[END_REF].

Conclusion

We have seen that real-time systems can reach steadiness over a finite and quantifiable amount of time, and that a necessary condition to assure this stability is that the mean utilization is smaller than 1. In practice, systems with a mean utilization ū > 0.9 have many deadline misses, which requires from system designers to quantify deadline miss probabilities carefully by using the distributions of response times provided in this paper. No schedulability tests considering the steadiness of response times exist. This is a natural step in our opinion for the use of the analysis provided in this work. For example, the analysis we have presented in this paper is well suited for an application of a Monte-Carlo response time analysis [START_REF] Bozhko | Monte carlo response-time analysis[END_REF] which has recently been proven efficient.

We have expressed the probability function of response times in a specific family of distributions, the inverse Gaussian distribution, and provided a method to generate them. However, the distribution functions of execution time need to be known, but they are usually unknown. The methods built in this paper could be used in empirical and measurement-based methods, for example using clustering methods [START_REF] Friebe | Identification and validation of markov models with continuous emission distributions for execution times[END_REF][START_REF] Zagalo | Identification of execution modes for real-time systems using cluster analysis[END_REF].

(a) Histogram of the generated response times of τ 3 in Table 1, as defined in Section 6.1.

(b) Response time distribution functions of τ 3 in Table 1.

  mean utilization of Γ , and, the maximum utilization of level k, u max k possible value of C i and t min i is the minimum possible value of T i , and the total maximum utilization u max ≜ τ k ∈Γ c max k

Figure 1

 1 Figure1: Level 3 demand of 1 000 instances of the DK model system, the heavy-traffic demand process and the classical deterministic worst-case analysis (WCET) considering only the maximal execution time for each task, for ū = 0.838 and u max = 1.208 and hyper-period T , for Γ defined in Table1.

  be the activation time of the l-th job of level k, and I k l the index of the task of the l-th job of level k. If ūk < 1, then the steady-state backlog βk exists, is finite and we obtain βk = sup n≥0 n l=1

Figure 2 :

 2 Figure 2: The heavy-traffic backlog process of systems with different mean utilization, initialized with β(0) = τi∈Γ C i .

Figure 3 :

 3 Figure 3: Steady-state backlog simulations with different mean utilizations in the DK model. In grey the histogram of simulations of β(n)/ √ n for n = 10000, in black the probability function of β.

Figure 4 :

 4 Figure 4: Trajectories of Ŵ and average of the first idle time for ū = 0.838.

Figure 5 :

 5 Figure 5: First maximum ε-idle time for ε = 10 -6 . for ū = 0.838.

Figure 6 :

 6 Figure 6: Representation of the response time R 2,1 as an idle time when A 2,1 = 0, C 2,1 = y and β 1 (0) = x, in Example 3.

Proposition 8 .

 8 The probability function of R k is h max k (resp. hk ).

Figure 7 :

 7 Figure 7: Simulations of a 10 000 instances for the SimSo simulations, steady and transient, EVT estimation of the worst-case response time of the SimSo simulations and simulations of heavy-traffic worst-case response time and steady-state response-time when ū = 0.838.

Table 1 :

 1 Task set used in the simulations of the experimental results.

	task τ k T k	C k	p k	ūk	u max k
	τ 1	4	(1, 2)	(0.5, 0.5)	0.375	0.5
	τ 2	6	(1, 2)	(0.5, 0.5)	0.625 0.833
	τ 3	8	(1, 2, 3)	(0.5, 0.3, 0.2)	0.838 1.208
	τ 4	10	(1, 2, 3)	(0.6, 0.2, 0.2)	0.998 1.508
	τ 5	12 (1, 2, 3, 4) (0.5, 0.3, 0.1, 0.1) 1.148 1.841
	quantify analytically this difference.

A real-time system is deterministic if all task parameters are described by unique values.

The ratio of the standard deviation and the mean.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gereme.html


Hakim Elchaoui Elghor
email: hakim.elchaoui@ibisc.univ-evry.fr

David Roussel

Fakhr-Eddine Ababsa

Houssine El

Bouyakhf

El-Houssine Bouyakhf

3D planar RGB-D SLAM system

Keywords: RGB-D SLAM, Planar features, 3D Plane-based maps

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In recent years, mobile robots have received an increasing interest especially for indoor environment applications. In these environments RGB-D sensors offer an opportunity to significantly develop robotic navigation and interaction capabilities. Since they combine the advantages of RGB cameras with the ability to obtain geometric information, many works tend to exploit the potential of these novel sensors. For applications such as Simultaneous Localization and Mapping (SLAM), introducing RGB-D cameras allows to create three-dimensional maps in real time. Several RGB-D SLAM systems have been proposed. They can be placed in two large family: Sparse and Dense SLAM systems. Although the purpose is the same, the two approaches diverge in the modeling and processing. Sparse SLAM approaches are based on visual odometry. They use visual features correspondences with registration algorithms, as RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]orICP [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF], to estimate and refine transformations between poses. The algorithm developed by Henry et al. [START_REF] Henry | RGB-D mapping: Using kinect-style depth cameras for dense 3D modeling of indoor environments[END_REF] was one of the first RGB-D systems using features points to estimate camera poses and represent the environment by surfels [START_REF] Pfister | Surfels: surface elements as rendering primitives[END_REF]. It creates and optimizes a Graph-Based SLAM. This modeling [START_REF] Grisetti | A tutorial on graph-based SLAM[END_REF] consists in constructing a graph which nodes are sensor poses and where an edge between two nodes represents the transformation (egomotion) between these poses. Such formulation enables the graph optimization step which aims to find the best nodes configuration that produces a correct topological trajectory and easier loop-closures detection when revisiting the same areas. Endres et al. [START_REF] Endres | 3D mapping with an RGB-D camera[END_REF] followed the same path and proposed a graph-based RGB-D SLAM which became very popular among Robotic Operating System (ROS) users due to its availability. The implementation and optimization of the pose-graph is performed by the G 2 o framework [START_REF] Kümmerle | G2o: a general framework for graph optimization[END_REF]. To represent the environment, they used 3D occupancy grid maps generated by the OctoMapping approach [START_REF] Hornung | Octomap: an efficient probabilistic 3d mapping framework based on octrees[END_REF]. This RGB-D SLAM system offers a good trade-off between the quality of pose estimation and computational cost. Indeed, sparse SLAM approaches are typically fast due to the sensor's egomotion estimation based on sparse points. In addition, such a lightweight implementation can be embedded easily on mobile robots and small devices. However, the reconstruction quality is limited to a sparse set of 3D points which leads to many redundant and repeated points in the map and lack of semantic description of the environment.

Instead, dense SLAM methods enable good pose estimation and high quality scene representation. However, they tend to drift over time and fail against scenes with poor geometric structure. To overcome high computational costs, these approaches use sophisticated equipments such as high performance graphics hardware GPU which may constrain the used platform. Dense RGB-D SLAM systems were introduced by Newcombe et al. in the well known Kinect Fusion [START_REF] Izadi | Kinectfusion: real-time 3d reconstruction and interaction using a moving depth camera[END_REF][START_REF] Newcombe | Kinectfusion: real-time dense surface mapping and tracking[END_REF]. It is a dense system integrating real-time depth measurements into a volumetric data structure to create highly detailed voxel-based maps. Camera poses are estimated by tracking only live depth data with a dense iterative closest point (ICP) algorithm. Despite the high quality maps, the algorithm fails in environments with poor structure and is restricted to small workspaces due to high memory consumption. Whelan et al. proposed a moving volume method [START_REF] Whelan | Real-time large-scale dense rgb-d slam with volumetric fusion[END_REF] to overcome the restricted area problem. By moving the voxel grid with each current camera pose, real-time detailed mapping of unbounded areas became possible. Unlike voxel-based reconstruction, Keller et al. [START_REF] Keller | Real-time 3d reconstruction in dynamic scenes using point-based fusion[END_REF] proposed a more efficient solution. They proposed a point-based fusion representation supporting spatially extended reconstructions with a fused surfel-based model of the environment.

Recently, perceiving the geometry surrounding robots has become a research field of great interest in computer vision. For both robot planing and augmented reality applications as [START_REF] Ababsa | Robust camera pose tracking for augmented reality using particle filtering framework[END_REF], using some geometric assumptions is a crucial prerequisite. Likewise, in current RGB-D SLAM systems researchers begin to pay a significant interest to geometric primitives in order to build three-dimensional (3D) structures. The observed geometry can be a good solution to better constrain the problem and help improve 3D reconstructions. For indoor environments, 3D planes can be relevant as they are extremely common and easily deduced from point clouds. Thus, they were introduced in several recent works. One of the earliest RGB-D SLAM approaches incorporating planes has been developed by Trevor et al. [START_REF] Trevor | Planar surface slam with 3d and 2d sensors[END_REF]. They combined a Kinect sensor with a large 2D planar laser scanner to generate both lines and planes as features in a graph based representation. Data association is performed by evaluating the joint probability over a set of interpretation trees of the measurements seen by the robot at one pose. Taguchi et al. [START_REF] Taguchi | Point-plane SLAM for handheld 3D sensors[END_REF] presented a bundle adjustment system combining both 3D point-to-point and 3D plane-to-plane correspondences. Their system shows a compact representation but a slow camera tracking. This work was extended by Ataer-Cansizoglu et al. [START_REF] Ataer-Cansizoglu | Tracking an rgbd camera using points and planes[END_REF] to find point and plane correspondences using camera motion prediction. However, the constant velocity assumption used to predict the pose seems to be difficult to satisfy when using handheld camera. The RGB-D SLAM system [START_REF] Keller | Real-time 3d reconstruction in dynamic scenes using point-based fusion[END_REF] was extended by Salas-Moreno et al. [START_REF] Salas-Moreno | Dense planar slam[END_REF] to enforce planarity on the dense reconstruction with application to augmented reality. In a recent work, Xiang and Zhang [START_REF] Gao | Robust rgb-d simultaneous localization and mapping using planar point features[END_REF] proposed an RGB-D SLAM system based on planar features. From each detected 3D plane they generate a 2D image and try to extract its 2D visual features. These extracted features are back-projected on the depth image to generate 3D feature points used to estimate the egomotion with ICP. More recently, Whelan et al. [START_REF] Whelan | Incremental and batch planar simplification of dense point cloud maps[END_REF] performed incremental planar segmentations on point clouds to generate a global mesh model consisting of planar and non-planar triangulated surfaces. In [START_REF] Kaess | Simultaneous localization and mapping with infinite planes[END_REF], the full representation of infinite planes is reduced to a point representation in the unit sphere S3. This allowed to parameterize the plane as a unit quaternion and formulate the problem as a least-squares optimization of a graph of infinite planes.

In our works, we are also focused on searching alternative 3D primitives in structured indoor scenes. Especially, we use 3D planar regions to generate a reduced significant representation of the environment in sparse RGB-D SLAM systems. Indeed, an RGB-D point cloud contains 307200 points and requires 3.4 Megabytes in memory. Together with the sensor's characteristic noise, the large number of points lead to significantly redundant 3D maps when assembling several point clouds. Thus, using planar assumptions on the observed geometry can deal with sensor noise and redundant representation of the environment. This paper, based on our previous work [START_REF] Elghor | Planes detection for robust localization and mapping in RGB-D SLAM systems[END_REF], describes an RGB-D SLAM system benefiting from planes segmentation to generate lightweight 3D plane-based maps. Unlike previous work, we propose a new method for building low-cost 3D maps based on reliable estimations. As human living environments are mostly composed of planar features, such technique is suitable to overcome the sensor's weakness without using a dense approach.

Our contributions in this paper are three-fold: i. We detail our 3D plane-based SLAM system implementation. ii. We propose a 3D compact planar modeling for indoor environment representations. iii. We introduce a new protocol to evaluate the system performance over real world scenes.

In the reminder of this paper we give an overview of our system in the following section. In Sect. 3 we provide an evaluation of the qualitative and quantitative performance of our system using a handheld RGB-D sensor. Finally, we summarize and report future works in Sect. 4.

System Description

The schematic representation of our system is shown in Fig. 1. Our starting point was inspired by the RGB-D SLAM system introduced by Endres et al. [START_REF] Endres | 3D mapping with an RGB-D camera[END_REF]. Inputs are color images and depth maps (RGB-D data) provided by the Kinect sensor. Our system introduces 3D planes to estimate camera's transformations and generate a planar reconstruction of the environment by merging planes into a global map. We begin by extracting local 3D planar regions (planes models with their bounding boxes) π l from Depth maps (See next subsection) and 2D feature points from RGB data concurrently. 2D features points belonging to detected planar regions are projected onto these planes as 3D planar feature points. These planar feature points are used to estimate the local transformation T l between two camera poses and refine it with an iterative RANSAC. Let the global coordinates system be the first registered camera pose. In each pose addition we store the global transformation T w leading to this coordinates system using current and previous local transformations. Then, global detected planar regions π w can be obtained. The 3D planar map is updated either by merging detected planar regions to the existing ones (Planes set π), or by adding new planes.

Planar Regions Detection and Representation

Planar regions detection is performed using the Point Cloud Library (PCL) [START_REF] Rusu | 3d is here: point cloud library (PCL)[END_REF]. We extract co-planar points sets into the full point cloud. To find the k most prominent planar regions, which contains the maximum number of 3D points (inlier points), we proceed to a multi-plane segmentation with an iterative RANSAC. At each stage, we detect the main planar region and remove its inliers from the point cloud. Then, we extract the next main planar region while the number N of its inliers is still significant (at least 700 points).

A detected planar region π is parameterized by its plane model (n x ,n y ,n z ,d)

⊤ where n ⊤ (n x ,n y ,n z ) is the normal vector and d the distance from the origin. A 3D point p(x, y, z) lying on this plane model satisfies the familiar equation n x x + n y y + n z z = -d, otherwise written n ⊤ p = -d. This equation will be used to generate our 3D planar feature points detailed in the next subsection.

To represent a plane in the map, we also need its bounding box B. When a planar region is detected in the local frame, the point cloud of its inlier points is stored and then projected into the world coordinates after the egomotion estimation. Then, we generate the matrix of the second central moments of this point cloud. Afterwards, we proceed to a Singular Value Decomposition (SVD) of this matrix in order to find the main axis vectors of the point cloud and therefore the bounding box according to these vectors.

In the sequel, each detected planar region will be formerly defined by π[n, d, N, B] encompassing the plane model (n, d), its inliers population N,a n d its bounding box B. These parameters will be useful during the global planar map construction (Sect. 2.3).

3D Planar Feature Points

Studies conducted in [START_REF] Gao | Robust rgb-d simultaneous localization and mapping using planar point features[END_REF][START_REF] Taguchi | Point-plane SLAM for handheld 3D sensors[END_REF] showed that planar primitives are safer and robust to noise which leads to more accurate pose estimation while making processing faster. Then, it is relevant to introduce planar feature points in our system since depth maps provided by the Kinect sensor are noisy and contain missing depth values. The aim is to minimize 3D points measurement errors by benefiting from 3D plane fitting. Thus, we define our 3D planar feature points as 3D feature points satisfying a detected 3D planar model. The generation of these 3D planar feature points begins with 3D planar regions detection from point clouds and visual 2D feature points extraction from RGB images concurrently. For each 2D feature point we retrieve the depth value, generate its 3D feature point and check it against detected planar regions. 3D feature points belonging to a planar region are kept and others are rejected. In addition, we perform a regularization step by projecting all remaining 3D planar feature points into their respective planar regions.

Planar-Maps Building

As mentioned before, our goal is to produce lightweight 3D planar maps which can be useful for indoor robot navigation and augmented reality applications. For low-cost applications such as small robots, this choice can be very efficient to avoid 3D point clouds representations which are highly redundant and require memory resources. Based on the detected planar regions and the estimated transformation in each added pose, our system adds and grows planar structures in the map over time. To construct the global map, planes must be represented in the same 3D coordinates system. As the planes detection is performed in local frames, registered planar regions can be transformed from the camera to world coordinates systems using the global egomotion estimation T w . If the matrix M (Rotation R and Translation t) represents this global transformation, a 3D point p w into the world coordinates can be found easily using its correspondent point in the camera p c by the well known equation: p w =R p c + t and conversely p c = R ⊤ p w -R ⊤ t. Then, a plane in the world coordinates system is defined by its normal vector and distance π w (n w , d w ) with n w = Rn c and d w = d c -n ⊤ c R ⊤ t. Once these parameters are defined, we proceed to a plane-to-plane comparison in order to update the global map. Whenever a new planar region matches to a registered one, we update the plane by merging their plane models. Correspondence between planar regions is examined in two stages (Fig. 2). Firstly, angle between two compared plane models must not exceed a threshold set to 10 • and the distance limit between them is set to 5 cm. When this is satisfied, we check that the two planar regions overlap or are spatially close at least. This condition is added to avoid merging planes which have the same model but don't represent a continuous plane in the real world. An example of such situation can be a wall with an open door in the middle: the two parts of the wall have the same equation. Thus, using the bounding boxes of concerned planar regions we obtain the overlap between them if it exists. In case there is no overlap, we compare the minimum distance between the two bounding boxes to a distance threshold of 50 cm. When two planar regions match, they are merged together in a new resulting plane according to their inlier points populations. Their bounding boxes are compared and extrema are assigned to the merged planar region. If no correspondence can be found, the detected region is added to the map as a new plane. Even more than planes themselves, our map contains theoretical intersections between these planes. Planes intersections are generated using an adjacency criterion. We represent this intersection by lines and points when two or three planes intersect. This makes our map more significant and workable for other applications, and represents a first step towards a more semantic map.

Evaluation and Results

In this section we evaluate our work with a series of experiments using RGB-D data acquired with a Kinect. We present real-time results obtained on a PC with an Intel Core i5-2400 CPU at 3.10 GHz × 4. Planes detection was performed within depth images with a plane thickness threshold set to 1.5 cm.W e rejected all points exceeding this threshold. Also, only the first three main planes were considered during planar regions detection. Since planes detection produces larger planes first (in terms of population), adding further planes may introduce a lot of irrelevant small planar patches. These thresholds were chosen with preliminary experiments to provide accurate results. Experimental protocol and results concerning planes quality, map building and poses trajectory are shown in the following subsections. Unfortunately similar systems do not offer public softwares which makes comparisons impossible as we propose novel evaluations.

Planes Detection Process

In order to determine the overall performance of our system, we evaluated the influence of 3D data size on the detected planes quality. By default the Kinect images resolution is 640 × 480. Hence, a subsampling step seems essential to limit the plane detection process time. Also, this allows online data acquisition and processing with the Kinect 30 Hz update rate. Otherwise, the point cloud will contain 307200 points, which requires several seconds for planes detection. This is often used by sparse systems to overcome computational time. So we used a subsampling factor when creating a point cloud. This subsampling factor reduces the depth image dimensions. Then, we studied the impact of subsampled data on the quality and speed of planes detection over different distances. We performed specific experiments on a plane placed in front of the Kinect. For each distance we changed the subsampling factor and observed the detected planes. Table 1 (page 8) details results of these experiments. For each experiment, we considered planes detection runtimes, estimated distances and the number of erroneous detected extra planes. From these three columns we can deduce that the subsampling factor 4 gives a good trade-off between runtime and planes estimation quality. We can notice that planes detection quality decreases with the distance as the depth noise becomes greater than the planar region thickness threshold. As the depth noise increases with distance, planes detection over 3.5 m cannot be performed. Analysis of these results allowed us to choose subsampling data with factor 4 as it doesn't degrade the detected planes quality and enables real time processing. We also noticed that angle errors between detected planes from the same distance does not exceed 6 • in any case.

3D Planar Maps

Here we show results of our planar SLAM system for a real time captured office scene. The Kinect was mounted on a movable support that enables it to be placed everywhere within the scene. The office scene contains several planes with various sizes set in different locations and features several parallel and intersecting planes. Figure 3 shows an example of our 3D plane-based map with computed intersections. Each plane in the map was generated progressively by the merging process described above. Theoretical intersections between adjacent planes was also updated whenever a plane model changed. Benefiting from this minimal representation, we generated a lightweight 3D global plane-based map unlike the usual heavyweight point-based maps.

Poses Estimation Error. In our previous work [START_REF] Elghor | Planes detection for robust localization and mapping in RGB-D SLAM systems[END_REF] we evaluated generated trajectories and runtimes of our approach against benchmark data (Sect. 4.1 Benchmark datasets). Evaluations have shown that we significantly reduced the egomotion run-time while keeping a good accuracy of the estimated trajectories over the compared approaches. Here we complet this study by introducing our novel evaluation protocol consisting in revisiting particular locations in the scene several times and comparing stored poses. As we don't have an external tracking system, the idea is to consider the first passage at a location as the ground truth and to compare further visits poses. We chose five sparse locations in the scene (crosses in Fig. 4 left) and performed a series of experiments where we moved through these locations randomly. Such experiments are very appreciated as by evaluating the revisited areas we evaluate loop closures and thus the system's accuracy. Table 2 summarises the Absolute Trajectory Error (ATE) and the Relative Pose Error (RPE) proposed by Sturm et al. [START_REF] Sturm | A benchmark for the evaluation of RGB-D SLAM systems[END_REF]. ATE computes the absolute difference between the ground truth and estimated poses after alignment and RPE measures the local accuracy of the trajectory over a fixed time interval ∆ (equals 1 here). Results show poses estimation accuracy and consequently the robustness of our system. This confirms the 3D planar feature points reliability already shown in our previous work.

Evaluations of the Map. As our goal is to build planar maps which can be used in several contexts such as mobile navigation and augmented reality, it is very important to assess our map against real world scenes. Our evaluation, first time proposed, consists in comparing measurements of the generated maps with those of real scenes, namely distances and angles between planes. Tables 3 and4 represent distances and angles estimations between numbered planes in Fig. 4 compared to the real ones in the office scene. Evaluations show that errors between real and estimated measurements do not exceeds 6 % for distances and 1 % for angles which is suitable for application requiring good accuracy. Hence, this lightweight map presents a good tradeoff between quality and memory consumption which is required by mobile robots and small devices.

3D Representations Comparison. Figure 5 shows a very close view of the second plane in Fig. 4 with two different 3D representations. The top row presents a point-based representation (Raw point cloud) and at the bottom we show our plane-based map which is a 3D planar representation. The point-based representation of this single plane already contains 340 thousand 3D points including many overlapping points due to assembling several point clouds. Such a representation would require a heavyweight process for the visualization of the entire map which leads to memory inefficiency. Our lightweight plane-based map is ready and much more usable for semantic labeling than raw 3D points based maps.

Conclusion

In this paper we proposed a simple 3D planar maps representation for RGB-D SLAM systems. We generated a 3D map based on detected 3D planes on the scene instead of the heavyweight point-based representation. During the reconstruction process, discovered planar areas are appended to the plane-based map. Moreover, all matched planes are progressively merged together in this map to give a compact representation and more information about the scene. Thus, we generated dense planar maps with significantly reduced sizes without relying on a dense approach. Besides, evaluations show that our system is able to build good quality maps while generating accurate trajectory. Also, the best tradeoff between precision and process time for planes detection have been evaluated. Once our plane-based map is available, the next step will consist in building the semantics associated to these planes such as floors or walls and planes structures such as desks or rooms which represents a first step towards semantic maps.

Fig. 1 .

 1 Fig. 1. Main system pipeline.

Fig. 2 .

 2 Fig. 2. Planar map building.

Fig. 3 .

 3 Fig. 3. Example of a 3D planar map resulting from our system. (left) Generated planes and their corresponding point cloud. (right) Our planar map with planes intersections.

Fig. 4 .

 4 Fig. 4. Overhead view of (right) Point-based map with real trajectory inside and (left) Plane-based map with evaluation locations (crosses).

Table 2 .

 2 ATE and RPE poses estimation errors. ATE RMSE (m) 0.006 ± 0.003 Relative pose error (RPE) RMSE -Translation (m) 0.014 ± 0.002 Relative pose error (RPE) RMSE -Rotation (deg) 1.30 • ± 0.54 •

Fig. 5 .

 5 Fig. 5. Close overhead view of a planar region in the scene: (top) Raw point cloud representation and (bottom) 3D Plane-based map.

Table 1 .

 1 Impact of subsampled data on planes detection.

	Real	Subsampling Estimation Inliers Estimated	Number of erroneous
	distance (m) factor	runtime (s) points distance (m) detected planes
	1.5	1	2.37	234813 1.49	0
		2	0.23	58797 1.48	0
		4	0.07	14545 1.49	0
	2.5	1	7.26	125107 2.44	1
		2	1.78	32320 2.41	2
		4	0.4	8619 2.41	0
	3.5	1	7.57	62209 3.32	2
		2	1.72	14637 3.37	1
		4	0.4	1904 3.36	0

Table 3 .

 3 Real and estimated distances between planes.Plane i Plane j Real distance (Avg. ± Std. Dev.) Estimation (Avg. ± Std. Dev.)

	2	3	1.50 m ± 0.02 m	1.45 m ± 0.06 m
	2	5	2.82 m ± 0.02 m	2.70 m ± 0.08 m
	2	6	3.42 m ± 0.02 m	3.40 m ± 0.06 m
	2	7	3.20 m ± 0.02 m	3.00 m ± 0.10 m
	2	9	1.95 m ± 0.02 m	1.81 m ± 0.10 m
	4	8	2.95 m ± 0.02 m	2.88 m ± 0.05 m

Table 4 .

 4 Real and estimated angles between planes. Plane i Plane j Real angle (Avg. ± Std. Dev.) Estimation (Avg. ± Std. Dev.)

	1	2	56 • ± 2 •	56.80 • ± 1.5 •
	3	4	90

• ± 2 • 92.37 • ± 2.1 • 4 5 90 • ± 2 • 90.12 • ± 0.7 • 7 8 90 • ± 2 • 92.52 • ± 1.4 • 8 9 90 • ± 2 • 91.75 • ± 0.9 •