Roberto Ferretti
email: ferretti@mat.uniroma3.it

Achille Sassi
email: ach.sassi@gmail.com

A semi-Lagrangian algorithm in policy space for hybrid optimal control problems *

Keywords: Hybrid control, Dynamic Programming, Semi-Lagrangian schemes, Policy iteration AMS Subject Classification 2010: 34A38, 49L20, 65B99, 65N06

The mathematical framework of hybrid system is a recent and general tool to treat control systems involving control action of heterogeneous nature. In this paper, we construct and test a semi-Lagrangian numerical scheme for solving the Dynamic Programming equation of an infinite horizon optimal control problem for hybrid systems. In order to speed up convergence, we also propose an acceleration technique based on policy iteration. Finally, we validate the approach via some numerical tests in low dimension.

Introduction

In the last two decades, the concept of hybrid control system has provided a sound mathematical framework for treating control systems in which continuous and discrete control actions mix together, and this framework has also been successfully adapted to optimal control problems. Among the various systems covered by this theory, we mention economic models with restocking, multigear and hybrid vehicles, and, more in general, systems with switchings in the dynamics and/or impulsive changes in the state. In this paper, we study efficient numerical methods for applying Dynamic Programming techniques to hybrid optimal control problems of infinite horizon type.

Among the various mathematical formulations of optimal control problems for hybrid systems, we will adopt here the one given in [START_REF] Branicky | A unified framework for hybrid control problem[END_REF][START_REF] Dharmatti | Hybrid control system and viscosity solutions[END_REF]. Let I = {1, . . . , m}, and consider the controlled system (X, Q) satisfying:

     Ẋ(t) = f X(t), Q(t), α(t) t ∈ (0, ∞) X(0) = x Q(0 + ) = q (1.1)
where x ∈ R d , and q ∈ I. Here, X and Q denote respectively the continuous and the discrete component of the state. The function f : R d × I × U → R d is the continuous dynamics, for a set of continuous controls given by:

U = α : (0, ∞) → U | α measurable, U compact .
The trajectory undergoes discrete transitions when it enters two predefined sets A (the autonomous jump set) and C (the controlled jump set), both of them subsets of R d × I. More precisely:

• On hitting A, the trajectory jumps to a predefined destination set D ⊂ R d × I. The jump driven by a transition map g : R d × I × V → D, where V is a discrete finite control set. Denoting by τ i a time at which the trajectory hits A, the new state will be X(τ

+ i , Q(τ + i ) = g X(τ - i ), Q(τ - i ), w i , for a control w i ∈ V.
• Entering the set C, the controller can choose either to jump or not. If the controller chooses to jump, then the continuous trajectory is moved to a new point in D. Denoting by ξ k one such time of jump, we will have

X(ξ - k ), Q(ξ - k ) ∈ C and (x , q ) = X(ξ + k ), Q(ξ + k ) ∈ D.
The trajectory starting from x ∈ R d with discrete state q ∈ I is therefore composed of branches of continuous evolution given by (1.1) between two discrete jumps at the transition times τ i or ξ k . Now, considering an optimal control problem in the infinite horizon form, and including all control actions in a control strategy θ := α(•), {w i } i∈N , (ξ k , x k , q k ) k∈N we associate to θ a cost defined by: J(x, q; θ) := +∞ 0 X(t), Q(t), α(t) e -λt dt

+ ∞ i=0 c A X(τ - i ), Q(τ - i ), w i e -λτ i + ∞ k=0 c C X(ξ - k ), Q(ξ - k ), X(ξ + k ), Q(ξ + k ) e -λξ k (1.2)
where λ > 0 is the discount factor, : R d ×I×U → R + is the running cost, c A : A×V → R + is the autonomous transition cost and c C : C × D → R + is the controlled transition cost.

The value function V of the problem is then defined as:

V (x, q) := inf θ J(x, q; θ).

We point out that, in this generality, the problem requires strong assumptions to be mathematically well-posed. In particular, it should be ensured that the value function (1.3) is continuous, and that the so-called "Zeno executions" (i.e., the occurrence of an infinite number of transitions in a finite time interval) are avoided. We will give in the next section a precise set of assumptions, whereas in the examples we will apply the numerical technique under consideration in more general situations, showing that the recipe is robust enough to handle them.

To the best of our knowledge, the first rigorous theoretical study of the convergence of numerical schemes for the approximation of the value function of (1.1)-(1.2) has been given in [START_REF] Ferretti | Monotone numerical schemes and feedback construction for hybrid control systems[END_REF]. Here, solvability of the scheme by value iteration is proved, along with a convergence result based on the Barles-Souganidis theorem [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second-order equations[END_REF]. In spite of its robustness, however, value iteration is a relatively inefficient technique to compute the numerical solution, and an acceleration strategy would be highly desirable. From the very start of Dynamic Programming techniques [START_REF] Bellman | Dynamic Programming[END_REF][START_REF] Howard | Dynamic Programming and Markov processes[END_REF], policy iteration (PI) has been recognized as a viable, usually faster alternative to value iteration in computing the fixed point of the Bellman operator. Among the wide literature on policy iteration, we quote here the pioneering theoretical analysis of Puterman and Brumelle [START_REF] Puterman | On the convergence of policy iteration in stationary dynamic programming[END_REF], which have shown that the linearization procedure underlying policy iteration is equivalent to a Newton-type iterative solver. More recently, the abstract setting of [START_REF] Puterman | On the convergence of policy iteration in stationary dynamic programming[END_REF] has been adapted to computationally relevant cases [START_REF] Santos | Convergence properties of policy iteration[END_REF], proving superlinear (and, in some cases, quadratic) convergence of policy iteration. Moreover, we mention that an adaptation of policy iteration to large sparse problems has been proposed as "modified policy iteration" (MPI) in [START_REF] Puterman | Modified policy iteration algorithms for discounted Markov decision problems[END_REF], and has also become a classical tool.

In the present paper, we intend to study the construction and numerical validation of a SL scheme with PI/MPI sover for hybrid optimal control. To this end, we will recall the general algorithm, sketch some implementation details for the simple case of onedimensional dynamics, and test the scheme on some numerical examples in dimension d = 1, 2.

The outline of the paper is the following. In Section 2 we will review the main results about the Bellman equation characterizing the value function, and construct a semi-Lagrangian (SL) approximation for V in the form of value iteration. In Section 3 we will improve the algorithm by a policy iteration technique. Finally, section 4 will present some numerical examples of approximation of the value function and construction of the optimal control.

A Semi-Lagrangian scheme for hybrid control problems

First, we recall some basic analytical results about the value function (1.3). To this end, we start by making a precise set of assumptions on the problem.

Basic assumptions and analytical framework

In the product space R d × I, we consider sets (and in particular the sets A, C and D) of the form

S = {(x, q) ∈ R d × I : x ∈ S i , q = i}, (2.1)
in which S i represents the subset of S in which q = i. We assume that:

(A1) For each q ∈ I, A q , C q , and D q are closed subsets of R d , and D q is bounded. ∂A q and ∂C q are C 2 .

(A2) The function f is bounded. Moreover, it is Lipschitz continuous in the state variable x and uniformly continuous in the control variable α.

(A3) The map g : A × V → D is bounded and uniformly Lipschitz continuous with respect to x.

(A4) ∂A is a compact set, and for some γ > 0, the following transversality condition:

f (x, q, α) • η x,q ≤ -2γ
holds for all x ∈ ∂A q , and all α ∈ U , where η x,q denotes the unit outward normal to ∂A q at x. We also assume similar transversality conditions on ∂C.

(A5) We assume that, for all i ∈ I, (A8) c A (x, q, w) and c C (x, q, x , q ) are uniformly Lipschitz continuous in the variables x and x , and bounded with a strictly positive infimum. Moreover, for any x and q, the function c C satisfies (for some ∆ ≥ 0) the inequality c C (x, q, x , q ) < c C (x, q, x, q) + c C (x, q, x , q ) -∆ Via a suitable generalization of the Dynamic Programming Principle, it can be proved that the Bellman equation of the problem is in the form of a Quasi-Variational Inequality, and more precisely, once defined the Hamiltonian by H(x, q, p) := sup α∈U -(x, q, α) -f (x, q, α) • p and the transition operators M and N by: Mφ(x, q) := inf w∈V φ g(x, q, w) + c A (x, q, w) (x, q) ∈ A N φ(x, q) := inf (x ,q )∈D φ(x , q ) + c C (x, q, x , q ) (x, q) ∈ C

d(A i , C i ) ≥ β > 0 and d(A i , D i ) ≥ β > 0,
we have the following Theorem 1 ( [START_REF] Dharmatti | Hybrid control system and viscosity solutions[END_REF]) Assume (A1)-(A8). Then, the function V is the unique bounded and Hölder continuous viscosity solution of:

       λV (x, q) + H x, q, D x V (x, q) = 0 (x, q) ∈ (R d × I) \ (A ∪ C) max V (x, q) -N V (x, q), V (x, q) + H x, q, D x V (x, q) = 0 (x, q) ∈ C V (x, q) -MV (x, q) = 0 (x, q) ∈ A (2.2)
Note that uniqueness follows from a strong comparison principle, which also allows to use the Barles-Souganidis theorem [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second-order equations[END_REF] for proving convergence of stable and monotone schemes.

Numerical approximation

In order to set up a numerical approximation for (2.2), we construct a discrete grid of nodes (x j , q) in the state space and fix the discretization parameters ∆x and ∆t. In what follows, we will denote the discretization steps in compact form by δ := (∆t, ∆x) and the approximate value function by V δ .

Following [START_REF] Ferretti | Monotone numerical schemes and feedback construction for hybrid control systems[END_REF], we write the fixed point form of the scheme at (x i , q) as v

(q) i = V δ (x i , q) =      min N V δ (x i , q), Σ(x i , q, V δ ) (x i , q) ∈ C M V δ (x i , q) (x i , q) ∈ A Σ(x i , q, V δ ) else (2.3)
in which N , M and Σ are consistent and monotone numerical approximations for respectively the operators N , M and the Hamiltonian H. More compactly, (2.3) could be written as

V δ = T δ (V δ ).
We recall that, for λ > 0, under the basic assumption which ensure continuity of the value function, the right-hand side of (2.3) is a contraction [START_REF] Ferretti | Monotone numerical schemes and feedback construction for hybrid control systems[END_REF] and can therefore be solved by fixed-point iteration, also known as value iteration (VI):

V δ,j+1 = T δ (V δ,j ). (2.4)
To define more explicitly the scheme, as well as to extend the approximate value function to all x ∈ R d and q ∈ I, we use an interpolation I constructed on the node values, and denote by I[V δ ](x, q) the interpolated value of V δ computed at (x, q). With this notation, a natural definition of the discrete jump operators M and N is given by

M V δ (x, q) := min w∈V I[V δ ] g(x, q, w) + c A (x, q, w) (2.5) N V δ (x, q) := min (x ,q )∈D I[V δ ](x , q ) + c C (x, q, x , q ) (2.6)
On the other hand, a standard semi-Lagrangian discretization of the Hamiltonian related to continuous control is provided (see [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF]) by

Σ(x i , q, V δ ) := min α∈U ∆t (x i , q, α) + e -λ∆t I[V δ ] x i + ∆t f (x i , q, α), q .
(2.7)

In the SL form, the value iteration (2.4) might then be recast at a node (x i , q) as

v (q) i,j+1 =          min w∈V I[V δ,j ] g(x i , q, w) + c A (x i , q, w) (x i , q) ∈ A min min (x ,q )∈D I[V δ,j ](x , q ) + c C (x i , q, x , q ) , Σ(x i , q, V δ,j ) (x i , q) ∈ C Σ(x i , q, V δ,j ) else (2.8
) with Σ given by (2.7), and j denoting the iteration number.

Convergence of the scheme can be proved by using the arguments in [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF][START_REF] Ferretti | Monotone numerical schemes and feedback construction for hybrid control systems[END_REF]) if the interpolation I is monotone (e.g., a P 1 or Q 1 finite element interpolation):

Theorem 2 ( [START_REF] Ferretti | Monotone numerical schemes and feedback construction for hybrid control systems[END_REF]) Assume (A1)-(A8). Assume in addition that λ > 0, and that the interpolation I is monotone and invariant for the sum of constants. Then, V δ,j → V δ for j → ∞. Moreover, the approximate solution V δ converges to V locally uniformly in R d × I for ∆x, ∆t → 0.

Policy iteration algorithm

Following [START_REF] Sassi | Tecniche di Programmazione Dinamica nell'ottimizzazione di sistemi di controllo ibridi[END_REF], we give now an even more explicit form of the scheme, which is the one applied to the one-dimensional examples of Sec. 4. Once we set up a 1-D space grid of evenly spaced nodes x 1 , . . . , x n with space step ∆x, the discrete solution may be given the vector structure v := (v (1) , v (2) , . . . ,

v (m) ) ∈ R nm in which v (q) := (v (q) 1 , . . . , v (q)
n ) denotes the discretized value function associated to the q-th component of the state space. Within the vector v, the element v

(k) i appears with the index (k -1)n + i.
Keeping the same notation for all vectors, α ∈ U nm will denote the vector of controls of the system, α

(k)
i being the value of the control at the space node x i while the k-th dynamics is active. We also define the vector s ∈ I nm representing the switching strategy, so that s (k) i = l means that if the trajectory is in x i and the active dynamics is k, the system commutes from k to l. Note that, in the numerical examples of Sec. 4, discontinuous jumps will always appear only on the discrete component of the state space, so that, for example, we have x = x and this data need not be kept in memory (we will use the term switch to denote a state transition of this kind).

In the general case, we would also need to keep memory of the arrival point of the jump and/or of the discrete control w in the case of an autonomous jump. In general, the arrival point is not a grid point, so that we also need to perform an interpolation in (2.5)-(2.6). Therefore, the details for the general case can be recovered by mixing the basic arguments used in what follows.

The endpoint of this construction is to put the problem in the standard form used in policy iteration, min

(α,s)∈U nm ×I nm B(α, s)v -c(α, s) = 0, (3.1)
with explicitly defined matrix B and vector c. Note that, in (3.1), we have made clear the fact that a policy is composed of both a feedback control α and a switching strategy s. Define now the matrices D A , D C ∈ M nm {0, 1} as permutations of the array v. These matrices represent changes in the state due to the switching strategy: D A corresponds to autonomous jumps and D C to controlled jumps. Note that, in our case, the elements of D A and D C will be in {0, 1}, that there exists at most one nonzero element on each row, and that the two matrices cannot have a nonzero element in the same position. In order to determine the positions of the nonzero elements d Aa,b (s) = 1 in the matrix D A , we apply the following rule. For all (i, k) ∈ {1, . . . , n} × I, if the following conditions hold:

     (x i , k) ∈ A s (k) i = k (a switch occurs) (x i , s (k) i ) ∈ g(x i , k, V) (the switch is in the image of g), then, a = (k -1)n + i b = (s (k)
i -1)n + i. Similarly, the nonzero elements of the matrix D C , d C a,b (s) = 1, follow a slightly less strict rule. For all {1, . . . , n} × I, if the following conditions hold:

(x i , k) ∈ C s (k) i = k (a switch occurs), then, a = (k -1)n + i b = (s (k) i -1)n + i. Last, we define the matrix D(s) := D A (s) + D C (s),
which accounts for changes in the state related to the switching strategy, both autonomous and controlled.

We turn now to the continuous control part. First, we write Σ in vector form as

Σ(x, k, v) = min α (k) ∈U n ∆t (x, k, α (k) ) + e -λ∆t E(x, k, α (k) )v (k) ,
where the matrix E(x, k, α (k) ) ∈ M n (R) is defined so as to have

E(x, k, α (k) )v (k) = I[V δ ] x + ∆t f (x, k, α (k) ), k and I[V δ ] x+∆tf (x, k, α (k) ), k and (x, k, α (k) ) denote vectors which collect respectively all the values I[V δ ] x i + ∆t f (x i , k, α (k) i ), k and (x i , k, α (k)
i ). At internal points, using a monotone P 1 interpolation for the values of v results in a convex combination of node values. On the boundary of the domain, the well-posedness of the problem requires either to have an invariance condition (which implies that f (x i , k, α (k) i ) always points inwards) or to perform an autonomous jump or switch when the boundary is reached. Therefore, we should not care about defining a space reconstruction outside of the computational domain, although this could be accomplished by extrapolating the internal values.

The matrix E(α, s) ∈ M nm (R) is then constructed in the block diagonal form:

E(α, s) :=      
E (1) (α (1) , s (1) ) 0

• • • 0 0 E (2) (α (2) , s (2) ) . . . . . . . . . . . . . . . 0 0 • • • 0 E (m) (α (m) , s (m) )      
Assuming for simplicity that we work at Courant numbers below the unity (although this is not necessary for the stability of SL schemes), each block

E (k) (α (k) , s (k) ) ∈ M n (R) is a sparse matrix with non-zero elements e (k)
i,j determined so as to implement a P 1 space interpolation, in the following way: for every (i, k) ∈ {1, . . . , n} × I, define

h i,k := ∆t ∆x f (x i , k, α (k) i ) and if s (k) i = k h i,k < 0 then e (k) i,i-1 (α, s) = 1 + h i,k e (k) i,i (α, s) = -h i,k else, if s (k) i = k h i,k > 0 then e (k) i,i (α, s) = 1 -h i,k e (k) i,i+1 (α, s) = h i,k
Note that, if a switching strategy z ∈ I nm doesn't perform any switch (i.e. z (k) i = k for all (i, k) ∈ N × I), by definition of the matrix E(α, s) we would obtain, for all k ∈ I,

E (k) (α (k) , z (k) )v (k) = E(x, k, α (k) )v (k) ,
whereas, in the general case, if a switch occurs at x i , then the corresponding element of the matrix E (k) is zero. Finally, we define the vector c(α, s) ∈ R nm with a block structure of the form : c(α, s) = c (1) (α (1) , s (1) ), c (2) (α (2) , s (2) ), . . . ,

c (m) (α (m) , s (m) ) with c (k) (α (k) , s (k) ) ∈ R n such that, for every (i, k) in {1, . . . , n} × I, c (k) i (α (k) i , s (k) i ) = -∆t (x i , k, α (k) i ) s (k) i = k -ξ(k, s (k) i ) s (k) i = k
where ξ(k, l) denotes the switching cost (c A or c C ) from dynamics k to l.

With these notations, we can write the SL scheme (2.3) in vector form as

v = min (α,s)∈U nm ×I nm D(s) + e -λ∆t E(α, s) v -c(α, s) , (3.2)
or, defined the matrix B(α, s) := -I nm + D(s) + e -λ∆t E(α, s), as min

(α,s)∈U nm ×I nm B(α, s)v -c(α, s) = 0.
Once we have reformulated the Semi-Lagrangian scheme for the hybrid control problem in the standard form, we can solve it using Algorithm 1. The only difference with a standard PI algorithm is to include the switching strategy in the control policy. We remark that some theoretical result obtained in the "classical" setting is also true in the hybrid setting. In particular, convergence might still be obtained by monotonicity (see, e.g., [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF]) with minor changes in the proof, since the right-hand side of (2.8) is still in the form of a minimum: Theorem 3 Let v be the solution of (3.2), and v j be defined by Algorithm 1. If min α,s B(α, s)v 0 -c(α, s) ≤ 0, then the sequence v j is monotone decreasing, and v j → v.

Algorithm 1 Policy Iteration, 1D matrix form

j ← 0 STOP ← FALSE α 0 ∈ U nm s 0 ∈ I nm while STOP = FALSE do if [stopping criterion satisfied] then STOP ← TRUE else v j ← w solution of B(α j , s j )w = c(α j , s j ) (Policy Evaluation) (α j+1 , s j+1 ) ← arg min (a,σ)∈U nm ×I nm B(a, σ)v j -c(a, σ) (Policy Improvement) j ← j + 1 end if end while

Modified policy iteration

A different iterative solver for the numerical scheme has been first proposed and analysed in [START_REF] Puterman | Modified policy iteration algorithms for discounted Markov decision problems[END_REF], and is known as modified policy iteration. It consists in performing the minimization in (2.8) only once every N it iterations. In other terms, the policy evaluation step is replaced by N it iterations of linear advection (in which, however, the transport may occur among different components of the state space). For N it = 1 we obtain the value iteration, whereas for N it → ∞ the transport steps converge to an exact policy evaluation, and the algorithm coincides with the previous "exact" PI algorithm.

The pseudo-code in Algorithm 2 shows the MPI algorithm in one-dimensional matrix form, for a comparison with the exact algorithm (Algorithm 1).

Algorithm 2 Modified Policy Iteration, 1D matrix form

j ← 0 STOP ← FALSE α 0 ∈ U nm s 0 ∈ I nm while STOP = FALSE do if [stopping criterion satisfied] then STOP ← TRUE else if j = 0 (mod N it ) then (α j+1 , s j+1 ) ← arg min (a,σ)∈U nm ×I nm B(a, σ)v j -c(a, σ) (Policy Improvement) else (α j+1 , s j+1 ) ← (α j , s j ) end if v j+1 ← D(s j+1 ) + e -λ∆t E(α j+1 , s j+1 ) v j -c(α j+1 , s j+1 ) (Inexact Policy Evaluation) j ← j + 1 end if end while
Note that, in the numerical test section, the MPI algorithm has been applied to the two-dimensional examples. Although the formulation in dimension d = 2 could be accomplished by a suitable redefinition of the vectors and matrices, in practice the MPI algorithm does not need such a formalism.

Concerning convergence, the hybrid case can again be treated with the same theoretical tools of the original proof in [START_REF] Puterman | Modified policy iteration algorithms for discounted Markov decision problems[END_REF], which relies on the monotonicity of the (discretized) Bellman operator, as well as on giving an upper and a lower bound on the sequence v j by means of two converging sequences (one of which generated by value iteration). More precisely, the sequence considered in the convergence proof for the MPI is the sequence of approximations obtained after each policy improvement. In our notation, this is the subsequence v l corresponding to j = lN it + 1. We have therefore the following Theorem 4 Let v be the solution of (3.2), and v j be defined by Algorithm 2. If

min α,s B(α, s)v 0 -c(α, s) ≤ 0,
then, for any N it ≥ 1, the subsequence v l obtained for j = lN it + 1 is monotone decreasing, and v l → v for l → ∞.

Numerical tests

We give in this section some numerical examples in one and two space dimensions, comparing the performances of Value and Policy Iteration -exact PI algorithm in one dimension, and MPI in two dimensions. The comparison shows a substantial improvement in the convergence of the solver for the exact PI algorithm, whereas the MPI performs roughly the same number of iterations as the VI. Here, the bottleneck is apparently the contraction coefficient of the Bellman operator. Nevertheless, the MPI allows to avoid the minimization step in a large majority of the iterates, thus reducing the CPU time. Note that in both two-dimensional examples the control appears only as a switching strategy, and the complexity of policy evaluation steps is reduced by a factor 1/m. For more complex control actions the improvement in computing time would be even greater.

Stabilization of an unstable system

We now apply this technique to a stabilization problem: we consider a system with two dynamics: one "strong and expensive" and the other "weak and cheap". Only the former is able to keep the state of the system within the given set over time.

The state equation Ẋ(t) = f X(t), Q(t), α(t) is defined by

f (x, q, α) = x + d 1 α q = 1 x + d 2 α q = 2
where d 1 < d 2 and -1 ≤ α(t) ≤ 1 for every t in [0, +∞). Switching is mandatory only when the dynamics q = 1 is active and |X(t)| = 1, which implies that the state of the system belongs to the interval [-1, 1] for all t in [0, +∞).

Here and in what follows, c i,j denotes a constant switching cost from q = i to q = j, and the cost functional is defined as The values assigned to all the parameters are summed up in Table 1, whereas Table 2 reports the number of iterations required for given stopping tolerances. In the first three examples, the stopping criterion reads

(x, q, α) = x 2 + c 1 α 2 q = 1 x 2 + c 2 α 2 q = 2 d 1 d 2 c 1,2 c 2,1 c 1 c 2 λ t f 0.
v j -v j-1 ∞ < .
Note that, according to Table 2, squaring the tolerance makes the number of iteration N P of the PI algorithm increase linearly, which indicates roughly quadratic convergence, while the number N V for VI has a geometric behaviour as expected.

Figure 1 shows the optimal strategy obtained for ∆t = 0.0067, ∆x = ∆t||f || ∞ , X(0), Q(0) = (0.5, 1), and t ∈ [0, t f ]. This strategy consists in using the unstable dynamics q = 1 as long as the state belongs to the interval [-x, x] (where the value of x ∈ [-1, 1] depends on the given data). On the other hand, as soon as |X(t)| > x, the optimal choice is to switch from q = 1 to q = 2 in order to stabilize the system and force it back towards the origin, then switch again to the first dynamics which can be used at a lower cost.

m ρ 1 ρ 2 ρ 3 r c d τ ν 140 [kg] 0.06 0.09 0.12 0.2 [m] 0.3 10 [Nm] 6 • 10 3 [min -1 ] c α c x c i,j λ t f 1 0.5 0.1 i = j 0 i = j 1 10 [s]

Three-gear vehicle

In this test, we consider the optimal control of a vehicle equipped with a three-gear engine, focusing on the acceleration strategy and the commutation between gears. Physical parameters correspond to the italian scooter Piaggio Vespa 50 Special. The state equation for the speed of the vehicle is defined, for each gear q ∈ {1, 2, 3}, by

f (x, q, α) := 1 m T (β q x) rρ q α -c d x 2
where T (ω) := τ ω ν -ω ν 3 is the power band of the engine, τ and ν are respectively its maximum torque and r.p.m., β q := 60 rπρq is a conversion coefficient with ρ q := transmission shaft r.p.m. crankshaft r.p.m. , r is the radius of the wheel and c d the drag coefficient. The control α(t) ∈ [0, 1] represents the fraction of maximum torque used and the running cost is a linear combination of x and α: (x, α) = -c x x + c α α c x and c α are positive weights. Last, we define c i,j as the switching cost from from q = i to q = j.

The numerical results are obtained by assigning realistic values (Table 3) to the parameters defining the dynamics. The number of iterations is shown in Table 4 for various tolerances. The constant number of iterations obtained by PI might be due to the fact that optimal solutions (seem to) work with increasing values of q, this possibly meaning some sort of "causality" in the propagation of the value function.

Figure 2 shows the power band corresponding to our choice of τ and ν. Figure 3 shows the optimal solution obtained with ∆t = 0.027 [s], ∆x = ∆t||f || ∞ , (x, q) = (0.28 [m/s], 1) and t ∈ [0, t f ]. The optimal strategy is to reach the highest gear as fast as possible and then stabilize at a value α ≈ 0.5. A different scenario is shown in Fig. 4, in which we set the initial state to (x, q) = (14.58 [m/s], 1). Here, the control lets the vehicle slow down, then switches to the third gear in order to replicate the previous strategy.

Bang-Bang control of a chemotherapy model

In this test, we consider the control of a two-compartment model of tumor growth. For this model, and cost functionals of the kind we will consider below, optimal controls are known to be of bang-bang type (see [START_REF] Ledzewicz | Optimal bang-bang control for a two-compartment model in cancer chemotherapy[END_REF]). In this case, we can recast the problem in hybrid form, by considering an evolution in lack of chemotherapy (Q = 1):

Ẋ1 (t) = -a 1 X 1 (t) + 2a 2 X 2 (t) Ẋ2 (t) = a 1 X 1 (t) -a 2 X 2 (t) (4.1)
and a different evolution at full-dose chemotherapy (Q = 2):

Ẋ1 (t) = -a 1 X 1 (t) Ẋ2 (t) = a 1 X 1 (t) -a 2 X 2 (t). (4.2)
Here, the two compartments represent the number of cells at different stages of their lives, and the chemotherapy acts by preventing the generation of new tumor cells in the first compartment by inhibiting the mitosis of cells in the second compartment. The cost functional is defined as

J(x, q; θ) = ∞ 0 r 1 Ẋ1 (t) + r 2 Ẋ2 (t) + Q(t) -1 e -λt dt, (4.3)
in which Ẋ1 (t) and Ẋ2 (t) are given by (4.1)-(4.2) for respectively Q(t) = 1 and Q(t) = 2, and we have to minimize a combination between the growth of the tumor mass and the toxic effect of the drug on healthy cells (note that this latter term appears only when Q(t) = 2). Due to the geometric properties of the problem, Zeno executions cannot occur, and we can avoid to introduce a switching cost, which would have no practical meaning. Setting the switching cost to zero also causes the two value functions to coincide, i.e., V (x, 1) ≡ V (x, 2), and in this case a switch can occur at t = 0 + . While the general theory usually rules out this situation, no particular problems arise in this specific case. The values of the parameters are assigned as in Table 5, according to the current literature (see [START_REF] Ledzewicz | Optimal bang-bang control for a two-compartment model in cancer chemotherapy[END_REF]). Figg. 5-7 show the value function(s) of the problem, the optimal switching with respect to time and space and a sample trajectory starting from the initial state (x 1 , x 2 ) = (2, 1). The value function has been computed with a 100 × 100 grid on the domain [0, 2] 2 , and ∆t = 0.1. Note that there exists a clear discontinuity for the gradient of the value function, which corresponds to the switching curve in Fig. 7, which separates the black region, in which the optimal solution is Q(t) = 1, from the white region, in which the optimal solution is Q(t) = 2. The approximate optimal control shows a limit cycle in which a quasi-periodic switching between the two dynamics takes place.

Table 6 compares the two (VI and MPI) numerical solvers. Here and in the following test, the MPI algorithm has been implemented with N it = 10, and an initial block of 10 value iterations has been performed at the very start in order to provide a better initial guess. As remarked above, the Modified Policy Iteration algorithm performs essentially the same number of iterations than the value iteration algorithm, but at a lower cost.

DC/AC inverter

The last test presents a single-phase DC/AC inverter, whose conceptual structure is sketched in Fig. 8. In this device, a DC source generates an AC output by means of a suitable operation of the switches S 1 , . . . , S 4 , as well as a suitable choice of the three components (R, L and C) which appear in series in the RLC load. Following [START_REF] Chai | Hybrid feedback control methods for robust and global power conversion[END_REF], we consider as state variables X 1 = i L (the current through the inductor L, i.e., through the load) and X 2 = v C (the voltage across the capacitor C), the state equations being

     Ẋ1 (t) = V DC L (Q(t) -2) - R L X 1 (t) - 1 L X 2 (t) Ẋ2 (t) = 1 C X 1 (t). (4.4)
The physical meaning of the discrete state variable depends on the state of the switches S 1 , . . . , S 4 , and more precisely

Q(t) =      1 if S 1 , S 3 = OF F and S 2 , S 4 = ON 2 if S 1 , S 4 = OF F and S 2 , S 3 = ON 3 if S 2 , S 4 = OF F and S 1 , S 3 = ON .
The cost functional is defined so as to force the system to evolve (approximately) along an ellipse of the state space (see [START_REF] Chai | Hybrid feedback control methods for robust and global power conversion[END_REF]), namely With the parameters chosen (see Table 7), a ≈ 0.84, b ≈ 1.34 and the required output of the system would be given by two sinusoids of amplitude respectively 126 A for X 1 and 200 V for X 2 , both at the frequency of 1 Hz. The approximate solution has been computed on a 100 × 100 grid on the domain [-250, 250] 2 , with ∆t = 0.01, and state constraint boundary conditions have been treated by penalization, assigning a stopping cost of 5 • 10 8 on the boundary. The effect of the lack of full controllability is apparent in Fig. 9, which shows one component of the value function (they are practically undistinguishable from one another) and the optimal switching with respect to time. Fig. 10 shows the output (X 1 (t), X 2 (t)) of the controlled system, whereas, as an example, Fig. 11 reports the switching map of the second component of the state space. Here, the optimal solution is to keep Q(t) = 2 in grey regions, commute to Q(t) = 1 in black regions and to Q(t) = 3 in white regions.

Finally, Table 8 compares the two numerical solvers (VI and MPI). In this last test, the stopping condition has been computed on the relative l 1 update,

v j -v j-1 1 v j 1 < ,
to avoid problems with both high values of the solution and the occurrence of a discontinuity caused by the lack of controllability.

Conclusions

We have constructed and validated a Semi-Lagrangian scheme for hybrid Dynamic Programming problems in infinite horizon form. The numerical scheme has been made more efficient by a Policy Iteration type solver. Numerical tests performed on examples of varying complexity show that the scheme is robust and that the approximate optimal control policy obtained is stable and accurate, although the complexity remains critical with respect to the dimension of the state space. This validation suggests that this could be a feasible method to design optimizationbased static controllers in low dimension.

Figure 1 :

 1 Figure 1: Value function, trajectory and optimal control, weak-strong test

Figure 2 :

 2 Figure 2: Power band of the engine.

Figure 3 :

 3 Figure 3: Value functions, trajectory and optimal control, three-gear vehicle, first case.

Figure 4 :

 4 Figure 4: Value functions, trajectory and optimal control, three-gear vehicle, second case.

Figure 5 :

 5 Figure 5: Value function and optimal switching for the chemotherapy test.

Figure 8 :

 8 Figure 8: Abstract structure of the single-phase DC/AC inverter.

b 2 2 1 a 2 + x 2 2 b 2 -c 2 .

 222222 = c, in which the constants a and b are defined in terms of the physical parameters R, L, C and of the desired pulsation ω. This makes it natural to define the running cost as(x, q, α) = x ] 0.7 [Ω] 0.1 [H] 0.1 [F] 2π [s -1 ] 22500 1 0 200

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Value function and optimal switching for the DC/AC inverter.

Table 1 :

 1 Choice of parameters, weak-strong test

	5 2 0.2	0	0.25 4 1 20
		N V		N P
	10 -3	456	8
	10 -6 1147 10
	10 -12 2786 12

Table 2 :

 2 Number of iterations (VI and PI) for a given tolerance , weak-strong test

Table 3 :

 3 Choice of parameters, three-gear vehicle test

		N V	N P
	10 -3	337	6
	10 -6	586	6
	10 -12 1084	6

Table 4 :

 4 Number of iterations (VI and PI) for a given tolerance , three-gear vehicle test.

Table 5 :

 5 Choice of parameters, chemotherapy test

	2
	1.8
	1.6
	1.4
	1.2
	1
	0.8
	0.6
	0.4
	0.2
	0
	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	Figure 7: Optimal switching map for the chemotherapy test.

Figure 6: Trajectories X 1 (t) and X 2 (t) for the chemotherapy test.

Table 6 :

 6 Number of iterations for a given tolerance , chemotherapy test.

Table 7 :

 7 Choice of parameters, inverter

		N V	N P
	10 -3	469	469
	10 -6 13481 13491

Table 8 :

 8 Number of iterations for a given tolerance , DC/AC inverter test.

* This work is partially supported by the EU under the 7th Framework Programme Marie Curie Initial Training Network "FP7-PEOPLE-2010-ITN", SADCO project, GA number 264735-SADCO. For the second author, also by iCODE Institute project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02.


Corentin Damman

Gregory Edison

Fabrice Guet

Eric Noulard

Luca Santinelli

Jerome Hugues

Architectural Performance Analysis of FPGA Synthesized LEON Processors

Keywords: Performance Analysis, Embedded Systems, Benchmarking, Processor Synthesis, FPGA, LEON, Probabilistic Worst-Case Execution Time

Modern embedded processors have gone through multiple internal optimization to speed-up the average execution time e.g., caches, pipelines, branch prediction. Besides, internal communication mechanisms and shared resources like caches or buses have a significant impact on Worst-Case Execution Times (WCETs). Having an accurate estimate of a WCET is now a challenge. Probabilistic approaches provide a viable alternative to single WCET estimation. They consider WCET as a probabilistic distribution associated to uncertainty or risk.

In this paper, we present synthetic benchmarks and associated analysis for several LEON3 configurations on FPGA targets. Benchmarking exposes key parameters to execution time variability allowing for accurate probabilistic modeling of system dynamics. We analyze the impact of architecture-level configurations on average and worst-case behaviors.

I. INTRODUCTION

Timing constraints contribute defining the correctness of most of embedded systems. For example, the processor that manages the injection of your engine or the electronic flight controls of an aircraft have to provide results within welldefined time windows. Such systems are called real-time systems. In particular, hard real-time systems if timing constraints misses are not tolerated and soft-real-time systems if some timing constraints misses can be tolerated.

Recent improvements and optimization of processors induce huge complexity onto real-time systems. We are now facing multi-core architectures with internal communications, cache memories, interruptions, multi-functionalities, etc. Each of which has to be properly analyzed in order to guarantee both system correctness and predictability. Classical single value Worst-Case Execution Time (WCET), based on modeling internal architectures behavior, may have reached its limit. Indeed, complexity in micro-architecture and/or costs in modeling system interactions, may lead to a large pessimism for the worst-case scenarios.

Probabilistic approaches now emerge as alternatives to single value WCETs. The probabilistic real-time modeling assumes the task WCET as worst-case distribution, the probabilistic WCET (pWCET), able to upper-bound any possible task execution behavior. The pWCETs allow accounting for the probability of occurrence of worst-case conditions which could be vanishingly small [START_REF] Cucu-Grosjean | Measurement-based probabilistic timing analysis for multi-path programs[END_REF]. Hence, pWCETs may lead to important reduction of computing capability over-provisioning since they cope better with tasks actual behavior.

Contributions :

In this paper we investigate the impacts of CPU architecture elements on tasks execution behaviors. For that purpose we synthesize LEON3 processors with different architectures on two Xilinx FPGAs. We modify some core configuration parameters like cache policy, the branch prediction and the floating-point unit to determine their impact. We run some benchmarks selected from the Mälardalen WCET project [START_REF]Mälardalen WCET benchmark programs[END_REF] on the synthesized architectures. Finally, we analyze the results with the use of a measurement-based probabilistic timing analysis tool called DIAGXTRM [START_REF] Guet | On the reliability of the probabilistic worst-case execution time estimates[END_REF], to perform average and worst-case performance analysis with probabilistic models. We will publish a data pack with all the models for the final version of the paper, at the following URL: https:// forge.onera.fr/ projects/ syntheticbench Organization of the paper: In Section II and Section III we present the architectural parameters applied for synthesizing LEON3 processors. Section IV describes the benchmarks and the task execution conditions applied. In Section V we introduce the DIAGXTRM framework applied for computing the pWCETs and for timing performance analysis. In Section VI we present the results obtained in terms of both average performance and worst-case probabilistic models. Section VII is for conclusions and future work.

A. Related Work

Performance analysis of multi-core platforms is a quite established research field, [START_REF] Mhadhbi | Performance Evaluation of FPGA Soft Cores Configurations Case of Xilinx MicroBlaze[END_REF][START_REF] Mattsson | Evaluation of synthesizable cpu cores[END_REF]. The approaches proposed so far refer to measurements for evaluating the impact of architectural elements on system average performance, [START_REF] Jiang | Measurement-based research on cryptographic algorithms for embedded real-time systems[END_REF][START_REF] Njuguna | A survey of fpga benchmarks[END_REF].

The probabilistic timing analysis estimates pWCETs and can be either static-based or measurements-based. Static Probabilistic Timing Analysis (SPTA) [START_REF] Davis | Analysis of probabilistic cache related pre-emption delays[END_REF][START_REF] Altmeyer | Static probabilistic timing analysis for real-time systems using random replacement caches[END_REF] requires an exact model of the system in order to infer the probability law, thus computing the pWCETs. Measurements-Based Probabilistic Timing Analysis (MBPTA) relies on measurements and on the Extreme Value Theory (EVT) in order to provide pWCET distributions, [1-3-10]. The MBPTA does not require models of the system, but only measurements of its behavior.

II. MICROPROCESSOR DESIGN SPACE

Recent processors architectures rely on multiple building blocks. In the following, we consider those with an impact on the pWCET and present their key parameters.

• Cache memories have a direct impact on performance through instruction/data prefetching. If there is no cache, the system is forced to load each instruction from the memory and to do multiple loads and stores for each variable directly. The consequent saturation of the memory buses will have an impact on the instruction latencies.

• Cache replacement policies control how cache misses are handled, and how new pages are loaded. Different policies may either speed up or slow down a software.

• Branch prediction controls how instructions are prefetched in the pipeline. Such policies usually rely on heuristics.

• Math co-processor enables faster mathematical operations compared to software-emulated ones.

For this study, we implement LEON3 processors, which general architecture is presented in Figure 1, [START_REF]LEON3 32-bit processor core[END_REF]. The LEON3 is a 32-bit SPARC V8 processor, [START_REF] Gaisler | LEON3 Product Sheet[END_REF]. We synthesize each LEON3 platform configuration using the GAISLER GRLIB VHDL IP LIBRARY written in VHDL, [START_REF] Gaisler | GRLIB VHDL IP Core Library[END_REF], for the GAISLER GR-XC6S-LX75 and the XILINX ML605, which are two boards with XILINX FPGAs, DDR RAM, Flash, Ethernet, USB, UART, etc. The two FPGAs are flashed with various configurations, according to the four key parameters considered. The architecture and its configuration space are presented in the next subsections.

A. Architecture

The architecture we consider is based on two LEON3MP processors at 50MHz. We make the following design choices.

The integer unit is configured with SPARC V8 multiply and divide instructions. Emulated floating-point operations also benefit from this option. The integer multiplier is based on a 32x32 pipelined multiplier with a two-cycles latency.

For the cache system, LEON3 processors feature two separate instruction and data caches (Harvard architecture) with snooping and locking. The Instruction Cache (IC) is implemented as a multi-way cache with 4 ways, 4 kbytes/way and 32 bytes/line. Regarding the Data Cache (DC), it is implemented with 4 ways, 4 kbytes/way and 16 bytes/line. Four cache replacement algorithms are proposed: random, direct, Least-Recently-Replaced (LRR) and Least-Recently-Used (LRU). Local instruction and data RAM are disabled. We also choose to disable the Memory Management Unit (MMU) by default.

We implement a hardware module called LEON3 STATIS-TICAL UNIT (L3STAT). Eight counters are enabled, four per processor: the execution time counter (counted in CPU ticks at a frequency of 50MHz), the data and instruction cache miss counters and the branch prediction miss counter. In this work we make use of the execution time counter for evaluating each each configuration and its impact on the task executions.

1) Cache Memory: Cache memories are key elements of embedded systems since, depending on the state of the cache, task execution time could change. Both the worst-case timing analysis and the accuracy of the pWCET estimates would be impacted by the cache. We implement two cases: with caches (caches) and without caches (nocaches). In the first case the cache hierarchy is fully active, while in the second we disabled the twos caches for each processor. We expect a penalty on the execution time from the absence of cache.

2) Cache Replacement Policies: Cache replacement policies impact system performance i.e. task execution time. The way cache lines are replaced would affect next accesses; thus the latency of retrieving the information in cache changes. The replacement policies are critical to system variability especially with benchmarks which saturate the cache. The different replacement algorithms implemented are: the Random algorithm, which selects an item randomly and evicts it; the Direct algorithm, where the address of the new item is directly used to calculate its location in the cache; the LRR policy: it evicts the item least recently replaced; the LRU policy which evicts the item least recently accessed.

3) Branch Prediction (BP): The LEON3 is an advanced 7stage pipelined processor, which implies the use of an efficient BP. According to the GRLIB Configuration Help, the BP option would improve performance with up to 20%, depending on application. In order to verify that gain of performance and measure the impact of this option on the worst-case behavior, we run the benchmarks on a platform with branch prediction (BP) and without branch prediction (noBP) for comparison.

4) Floating-Point Unit (FPU): The FPU is a highperformance, fully pipelined IEEE-754 FPU. The two FPUs provided are the GAISLER RESEARCH'S GRFPU and GRFPU-LITE. The GRFPU is a high-performance pipelined FPU with high area requirements. GRFPU-LITE provides a balanced option with high acceleration of floating-point computations combined with lower area requirements compared to GRFPU, [START_REF] Gaisler | GRFPU -High Performance Floating Point Unit[END_REF]. They support all SPARC FPU instructions. If the FPU is disabled, a simulated, software FPU can be used. The option without FPU is identified noFPU, while that with FPU are identified as GRFPU or GRFPU-lite.

B. Architectural Configurations

A platform configuration is made of a set of architectural elements i.e. caches, cache replacement policies, BP and FPU.

We define the reference configuration as Reference = (caches, Random (IC), Random (DC), BP, noFPU). In the following, we make use of Reference for indicating also the choice on the element composing the reference configuration. For example, when talking about replacement policies, by referring to Reference we mean the random replacement policy. All the other configurations are compared with Reference.

Other configurations considered are: Direct = (caches, Direct (IC), Direct (DC), BP, noFPU), LRR = (caches, LRR (IC), LRR (DC), BP, noFPU), LRU = (caches, LRU (IC), LRU (DC), BP, noFPU), noBP = (caches, Random (IC), Random (DC), noBP, noFPU), nocaches = (nocaches, N/A, N/A, BP, noFPU), GRFPU = (caches, Random (IC), Random (DC), BP, GRFPU), GRFPU-lite = (caches, Random (IC), Random (DC), BP, GRFPU-lite)

The configuration of a typical System on a Chip (SoC) architecture like the GAISLER GR712RC is GR712RC = (caches, LRU (IC), LRU (DC), BP, GRFPU), [START_REF] Gaisler | GR712RC User's Manual[END_REF]. We will consider each choice of this architecture.

III. MICROPROCESSOR SYNTHESIS

The first step of this work is to synthesize the different architectural configurations presented in Section II. Figure 2 shows the entire methodology with the steps considered, respectively for synthesizing, benchmarking and analyzing. We use the GR-XC6S board for the architectural configurations Reference, Direct, LRR, LRU, noBP and nocaches; we use the ML605 board instead for the configurations GRFPU and GRFPU-lite to overcome limits in synthetizing FPUs. The GR-XC6S is a SPARTAN-6 FPGA, the ML605 board feature a Virtex-6 FPGA [START_REF]Spartan-6 Family Overview[END_REF][START_REF]Virtex-6 Family Overview[END_REF].

Figure 4 shows the occupations of the FPGAs in terms of slice registers (mostly Flip Flops), slice Look-Up Tables (LUTs, mostly used as logic) and the global numbers of occupied slices. The Table I presents detailed values.

We compare the different architectures related to Reference; as a reminder, Reference is the architecture with both data and instruction caches using random policy, without MMU, with branch prediction and without FPU.

• Removing the data and instruction caches reduce slices by 6%, but deteriorate run-time performances, see Section VI.

• We can see that the architectures using the Direct or the LRR replacement algorithm uses almost the same slices as the Random policy. However, the LRU policy occupies 15% more space. This area overhead is due to the addition of 5 flip-flops per line (for a 4-way LRU) to store the accesses.

• Disabling the branch prediction only saves 3% of the number of occupied slices, and reduces the performance.

• Regarding the implementation of the FPU on the ML605 board, we can see that the addition of the GRFPU uses almost twice the slices required without FPU. The GRFPU-LITE is, on the other hand, less demanding in terms of occupied slices.

IV. EMBEDDED SYSTEM BENCHMARKING

The main objective of the benchmarking is to reveal the impact of each architectural change on both average and worst-case execution time of tasks. We chose the Mälardalen University WCET project benchmarks to implement tasks and to push the processor to its limits.

• cnt and matmult are used to fill entirely the data cache with a lot of matrix operations. cnt counts non-negative numbers in a matrix with a parametrized size of 1000 × 1000. matmult performs a matrix multiplication. As for cnt, the size is parameterizable and we used 200 × 200 matrix. Those benchmarks are applied to measure the performance of the cache replacement policy and the impact of the lack of caches if they were removed.

• nsichneu is similar to the first two benchmarks, but this time for the instruction cache. It uses a large amounts of ifstatements (more than 250) to simulate an extended Petri Net. The benchmark is provided in three different versions, the nsichneu inner is considered.

• jfdctint is a typical application used in embedded real-time systems. By typical application we mean an application that does not saturate neither the data cache nor the instruction cache. jfdctint performs a discrete-cosine transformation on a 8×8 pixel block. The code is made of long calculation sequences (i.e., long basic blocks) and single-nested loops.

• lms is used for FPU performances comparison. The benchmark realizes a Least Mean Squares (LMS) adaptive signal enhancement on a sine wave with added white noise.

All benchmarks are single path tasks, and so no functional variability (input vector) is considered. Also there is no need to look for the worst case path.

We use GRMON to run all benchmarks, [START_REF] Gaisler | GRMON2 User's Manual[END_REF]. This monitor allows a non-intrusive debugging of LEON3 systems. It also performs communication with the modules implemented in the architecture, like the L3STAT unit. Finally, it is used to download and run applications on FPGAs.

A. Execution Conditions

The benchmarks are executed on all architectural configurations, with medium optimization -O2 and SPARC V8 multiply -divide instructions. We investigate two execution conditions:

The first execution condition is to run the task alone, without any other interaction. This execution condition is the baseline, as it is exempt of any interference tasks by construction. All benchmarks were compiled for this execution condition, and the flag -msoft-float was used to emulate floating-point operations when no FPU was present. This execution condition is named ELF after its compiler name.

The second execution condition consider using the RTEMS Real-Time Operating System (RTOS), [START_REF]RTEMS, An Open Real-Time Operating System[END_REF]. This execution condition is more representative of an actual real-time embedded system. It is compiled in Asymmetric Multi Processor mode, to exercise perturbations between cores. This execution condition is identified as RTEMS.

The RTEMS Multi-Processor application is running in the pre-emptive mode on two cores. In the first core, two tasks are launched: the benchmarking task with priority 2 and a second interference task, doing infinite loops, with priority 1 (higher priority). The second core contains another task made of an infinite loop. This minimal intra-and inter-core interference is used to determine the effect of interactions and is representative of a not heavily loaded system. The methodology could apply to any kind of interference, and it will be considered in future work.

B. Initial State Assumptions

The benchmarking execution diagram is presented in Figure 3. The trace of measurements is made of 5000 consecutive executions of the same benchmark. There exists a trace per benchmark, per configuration and per execution condition.

We have to ensure that the processor is in the same initial state at each benchmark execution. Before each run, the hardware counters from the L3STAT unit are reset, the caches are flushed and the pipeline is filled with nop() instructions. Then, the counters are activated, the benchmark is executed and finally the eight countersare disabled and printed for that execution. The resulting trace T is a matrix of size 8 × 5000 (eight counters, measured for each of the 5000 executions). The execution time is measured in CPU ticks and the trace duration spans from 4 minutes to 30 hours, depending on the benchmark to run as well as the architectural configuration. This allows including multiple systemic effects in the traces which could happen late in the executions.

V. PROBABILISTIC PERFORMANCE EVALUATION

The traces obtained from the benchmarking traces are applied to estimate the pWCETs, Figure 2.

Measurements-Based Probabilistic Time Analysis (MBPTA) allows defining tasks execution behaviors from runtime measurements. We use the MBPTA approach called DIAGXTRM to evaluate the impacts of architectural elements on both average and worst-case behavior of tasks. Only the execution time counter from the L3STAT traces is used.

The Extreme Value Theory (EVT) is the statistical tool composing any MBPTA that produces continuous distributions which are safe estimation of the task worst-case behavior 1 : the pWCETs. For this specific paper, DIAGXTRM applies the EVT in its peak over threshold version identifying the pWCETs from the Generalized Pareto Distribution (GPD) family. The distribution shape parameter ξ identifies the distribution that better cope with the measurements within this family. In particular, ξ < 0 is for a Weibull distribution, ξ = 0 is for a Gumbel distribution and ξ > 0 is for a Frechet distribution.

Hence, the pWCET is a worst-case thresholds with a probability associated to the risk to go past this value.

The EVT, thus the MBPTA, is applicable if the measurements are 1. stationary, 2. independent or 3. extremal independent and 4. they match a specified theoretical model, Figure 5. These situations correspond to four hypothesis that can be asserted from measurements. DIAGXTRM verifies all those hypotheses with the corresponding statistical tests. The confidence on the hypotheses validation is related to the quality of the pWCET model produced, named reliability of the pWCET. Only if all the tests succeed, the pWCET obtained is reliable and safe. To note that the independence and the extremal independence are partially overlapping hypotheses; if one is verified, there is no need to verify the other, Figure 5.

The robust statistics [START_REF] Guet | On the reliability of the probabilistic worst-case execution time estimates[END_REF][START_REF] Mhadhbi | Performance Evaluation of FPGA Soft Cores Configurations Case of Xilinx MicroBlaze[END_REF][START_REF] Mattsson | Evaluation of synthesizable cpu cores[END_REF][START_REF] Jiang | Measurement-based research on cryptographic algorithms for embedded real-time systems[END_REF][START_REF] Njuguna | A survey of fpga benchmarks[END_REF][START_REF] Davis | Analysis of probabilistic cache related pre-emption delays[END_REF][START_REF] Altmeyer | Static probabilistic timing analysis for real-time systems using random replacement caches[END_REF][START_REF] Hansen | Statistical-Based WCET Estimation and Validation[END_REF][START_REF]LEON3 32-bit processor core[END_REF][START_REF] Gaisler | LEON3 Product Sheet[END_REF][START_REF] Gaisler | GRLIB VHDL IP Core Library[END_REF][START_REF] Gaisler | GRFPU -High Performance Floating Point Unit[END_REF][START_REF] Gaisler | GR712RC User's Manual[END_REF][START_REF]Spartan-6 Family Overview[END_REF][START_REF]Virtex-6 Family Overview[END_REF][START_REF] Gaisler | GRMON2 User's Manual[END_REF][START_REF]RTEMS, An Open Real-Time Operating System[END_REF][START_REF] Buckley | Fuzzy statistics: hypothesis testing[END_REF] is applied by DIAGXTRM for quantifying the uncertainties from statistical tests. It defines the confidence cl on the hypothesis verified. 1 By safe estimation we mean a pWCET which is larger than or equal to any possible task execution time, [START_REF] Davis | Analysis of probabilistic cache related pre-emption delays[END_REF].

a) Stationarity Hypothesis: The analysis of trace can show the relationship between consecutive measurements and evaluate the impact that previous (in time) measurements on future ones, i.e. the independence between measurements. The EVT applicability relates to stationary traces. An established test for verifying this hypothesis is the Kwiatowski Phillips Schmidt Shin test, [START_REF] Kwiatkowski | Testing the null hypothesis of stationarity against the alternative of a unit root : How sure are we that economic time series have a unit root?[END_REF]. The confidence associated is cl 1 .

b) Independence Hypothesis: The statistical dependence is checked through the existence of correlated patterns of measurements, using the Brock Dechert Scheinkman statistical test [START_REF] Brock | A Test for Independence based on the Correlation Dimension[END_REF]. The confidence on the independence is cl 2.1 .

c) Extremal Independence Hypothesis: When overall independence does not hold, another way is to look for the independence of extreme measurements 2 . The extremal index θ ∈ [0, 1] is an indicator of dependence degree between extreme measurements, [START_REF] Northrop | An efficient semiparametric maxima estimator of the extremal index[END_REF]. We denote this parameter cl 2.2 .

d) Matching Hypothesis: The matching test is based on a quadratic statistic which measures the square distance between the pWCET model estimated and an hypothetical exact pWCET model following the GPD. The Cramer Von Mises criterion is applied for verifying the degree of matching between the two distributions, leading to parameter cl 3 .

Finally, DIAGXTRM computes pWCET parameters such as the threshold u and the pWCET distribution shape ξ that guarantee the best pWCET model among the possible EVT pWCETs, see [START_REF] Guet | On the reliability of the probabilistic worst-case execution time estimates[END_REF] for details.

We stress the fact that with MBPTA approaches, the resulting pWCET model represents the worst-case for the execution conditions accounted for by the measurements. It is not an absolute worst-case. Thus safety, confidence and reliability relates to the specific measurements considered.

VI. TRACE ANALYSIS

Given a trace of measurements, DIAGXTRM profiles both the average and the worst-case behavior of tasks. We use that tool for each architectural configuration and present the results in the following subsections.

A. Average Performance

At first we focus on the average execution time to compare the architectural configuration effects.

1) Execution Time Variability: Figure 6 illustrate the execution time variability of any system. More precisely, the trace is the result of 5000 executions of the cnt benchmark on the reference architecture without any RTOS. That proves that even when an application is run in a standalone mode, variability is present due to systemic and functional changes happening at runtime.

2) Impact of the Cache Replacement Policy: We analyze the traces of the two benchmarks matmult and nsichneu which saturate respectively the data and the instruction caches.

Figure 7a presents the box plots of execution times for the matmult benchmark. We can see that in order to improve the average performance of the data cache, the best cache 2 By extremal measurements we allude to observations relatively far from the average values as well as observations well separated in time.

replacement policy is the LRU algorithm. Then comes the random, the direct and finally the LRR algorithm.

As a matter of fact, programs usually compute results thanks to a restricted set of variables. It is then natural that the LRU algorithm is the best choice for the data cache replacement policy. This result is consistent with the fact that multiple SoC architectures implement the LRU algorithm for their data cache replacement policy, like the GR712RC configuration.

Figure 7b illustrates the execution times as box plots for the benchmark nsichneu. We can see that the LRU algorithm is not the best replacement policy for the instruction cache. The random algorithm is the first, then the direct, the LRU and finally the worst is the LRR.

For instruction cache the random replacement policy has better performance on average than the other policies, [START_REF] Altmeyer | Static probabilistic timing analysis for real-time systems using random replacement caches[END_REF]. Unlike data variables accesses, sequences of instructions of any program are not really regular. That explains why the random algorithm is the best instruction cache replacement policy. However, the choice of another algorithm is not critical as the drawback induced on the execution time is not significant, comparatively to the bigger differences of execution times related to the data cache replacement policies, Figure 7a (to note that the scales are not the same in the two figures). 3) Disabling Key Architectural Elements: In order to evaluate the gain of performance provided by the branch prediction or the cache memories, we disabled them from reference architecture. The Table II summarizes the execution time distributions and the performance drawbacks for different benchmarks. We chose matmult and nsichneu because they overuse the cache memories when present. Thus, those benchmarks are prone to reveal the catastrophic effects of disabling the instruction and the data caches. The third benchmark analyzed is jfdctint which is applied for investigating the impact on the execution time for an application that does not saturate caches. As described in Section II, the branch prediction deteriorates performances up to 20% if disabled. A maximal drawback of 11.19% is measured, Table IIa.

The cache memories impact on the execution time distribution is huge. For the same application, a factor of 8.78 on the execution time can be observed (Table IIb), where a factor of at least 2-3 was expected according to the GRLIB Configuration Help. For a typical application, the performance drawbacks are not negligible, Table IIc.

The results demonstrate that branch prediction and cache memories are key elements for embedded systems. Although they clearly improve average performances, they increase execution time variability. The system predictability or analyzability are affected by those elements.

4) Impact of The FPU Architecture: We use the benchmark lms that computes a lot of floating-point operations. Figure 8 presents the first 1000 executions of the resulting traces. We note that there are different execution modes, each with some variability, depending on the FPU architecture. The only element which is not reset between each benchmark execution is the FPU branch prediction, as the integer unit branch prediction and the cache memories are flushed each time. The multiple modes are then probably due to a non-zero prediction distance as detailed in [START_REF] Tyson | The effects of predicated execution on branch prediction[END_REF]. In the remaining 4000 measurements only one execution mode exists, revealing a convergence to a mode after some executions. Its mean value is used for comparison between the 3 FPUs. The emulated FPU has a mean level of 17 826 717 CPU ticks, the GRFPU converges around 522 888 CPU ticks and the GRFPU-LITE at 923 719 CPU ticks. This corresponds to improvements of respectively 33× and 18× faster.

The advantages of using a hardware FPU instead of an emulated FPU for floating-point applications are comforted. The GRFPU is up to 34 times faster than a software FPU. The GRFPU-LITE is a good compromise between speed improvement (19 times faster) and occupied slices (see Section III).

A side comment regards the capability that DIAGXTRM has of modeling architectural elements without the need for an exact element model. For example, with lms and the FPU architecture, although not knowing the internal behavior of the FPUs applied, DIAGXTRM is able to characterize the different modes just from measurements.

Table III summarizes the average execution time of each architectural change. This is done for the five benchmarks of interest and for each execution condition.

Two additional configurations were synthesized for each FPGA: theyare identified as Best. For the GR-XC6S it is Best1 = (caches, Random (IC), LRU (DC), BP, noFPU) and for the ML605 it is Best2 = (caches, Random (IC), LRU (DC), BP, GRFPU). They are the best configuration since they collect the choices which improves the most the average performances, respectively for each FPGA architectures.

On a global scale, the Best configurations improve the performance more than they deteriorate it, Table III. The architectures without BP or without caches are clearly the worst. The LRU algorithm is not always the faster, as we can see with the nsichneu benchmark, but its choice rather than the random replacement policy for the data cache is profitable. As we can see, Best configurations allows an additional gain of 1.01% on the GR-XC6S and of 8.36% on the ML605 in comparison with the second best architectures, respectively Reference and GRFPU.

Finally, the effect of the RTOS can be observed: the mean execution times undergo a penalty of 2.52% on average with RTEMS as compared to standalone ELF.

B. Worst-Case Performance

This section focuses on the worst-case performance and the model used for estimating the pWCETs. The methodology applied as well as the quality of the pWCET models (reliability and hypothesis confidence) are detailed. Two cases are investigated: the impact of different cache replacement policies and the impact of the RTOS on worst-case behaviors. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq qq q qq q q q q q q q q q q q q qq q q q q qq q q qqq q q q q q qq q q q q q 139350000 139400000 139450000 1) Probabilistic Distribution Fitting: Figure 9 shows the measurements and the pWCET derived from the EVT. The inverse cumulative distributions are represented. u is the threshold parameter that DIAGXTRM computes for obtaining the best pWCET model; k is the number of measurements above the threshold selected. We see how the pWCET model is able to perfectly cope with the input measurements and then safely infers the rare events (risk probability ≤ 10 -9 ).

2) Quality of the Probabilistic WCET Model: As described in Section V, the reliability of the pWCET model is based on the confidence levels of the four hypotheses {cl 1 , cl 2,1 , cl 2,2 , cl 3 }. Each hypothesis is scored out of 4.0, and a minimum of 1.0/4.0 is required in order to apply the method. Hypotheses cl 2.1 and cl 2.2 are covering each other, so that a score of 1.0/4.0 in one of the two is sufficient. The spider net in Figure 10 represents the reliability of the pWCET model from these confidence levels. 3 For each benchmark, a score is determined as executionT ime max(executionT imes) . Then, the sum of the scores are compared relatively to those of the Reference configurations.

Hence, all the assumptions are satisfied, the model is reliable. Stationarity, independence, extremal independence and matching are verified with almost the maximal confidence. The EVT is applicable with realistic architectures as there is enough variability from the measurements. 3) Comparisons of Worst-Case Behaviors: From the pWCET distributions we can extract the WCET thresholds with a 10 -9 risk probability i.e. confidence. This allows us to identify the impact of each architectural change on the worst-case behaviors with single values. Table IV indicates the pWCET distribution shape parameter ξ for each pWCET model as well as the theoretical estimations i.e. the WCET thresholds at 10 -9 . DIAGXTRM automatically selects ξ to best fit the input measurements. As it can be seen, the measurements have small variance hence they are better approximated with a Weibull GPD distribution, ξ < 0.

The accuracy of the pWCET is: accuracy IV. The small accuracy indicates that the measured values embed already the worst-case conditions. Thus, the pWCET as Weibull distribution is able to perfectly cope with the measurements and the task average behavior. The case with accuracy equal to 1.32E-01, Table IV ELF, indicates that the measurements have large variability. This means that the worst-case conditions could have a significant impact on the task execution times. The resulting Weibull, in order to account for that and foreseen the worst-cases, has to differ more from the measurements. Nonetheless, the pWCET is equally reliable due to the maximum confidences achieved.

The impact of the cache replacement policy on the worstcase behavior is presented in Table IVa. The results are consistent with those presented in Section VI-A2.

The benchmark lms was used to characterize the impact of the RTOS. We noted that the behavior of the FPU is not stable in the first 1000 runs of the benchmark, thus we use the last 4000 iterations to determine the pWCET distribution (i.e. achieve the hypotheses scores presented before in Figure 10b). The Table IVb gives the results. They confirm that the use of a RTOS like RTEMS induces a small penalty of 2.85% on the task performance: minimal interference were added and then accounted by a more pessimistic pWCET. VII. CONCLUSION AND FUTURE WORK In this work, we review i) configuration space of LEON processors on FPGA and ii) benchmarking LEON3 processors elements. Some architectural choices have been made, described and then synthesized on two FPGA platforms. Tasks benchmarks have been selected for their ability of revealing the impact of each architectural configuration on average and worst-case behaviors. Some execution conditions have been chosen in order to expose key parameters in execution time variability in both standalone and RTOS environments.

The results verify that time execution variability is present in any embedded processor. Average performance analyses show that the best data cache replacement policy is the LRU algorithm, while the random algorithm is the best choice for the instruction cache. The results confirm that cache memories, branch predictions and high-performance math co-processor are key elements of modern processors. Both contribute in improving average performance at the cost of increased variability in the system. We highlight configurations that are best compromise: they improve performance in almost every case, and provide guidelines for predictable embedded systems.

Probabilistic models have been estimated in order to see the impact of different cache replacement policies and the impact of the RTOS on the worst-case task behaviors. We validated the reliability and the accuracy of the resulting pWCETs.

Our approach could be applied with other benchmarks or other architectural configurations. Future work will study interference among tasks in RTOS, and evaluate the impact of scheduler or resource sharing at software and hardware levels (cache, MMU, etc.). This will improve our understanding of modern processors so as to guarantee time constraints.

Fig. 1 .

 1 Fig. 1. LEON3 Architectural Blocks.

Fig. 2 .

 2 Fig. 2. Methodology applied.

Fig. 4 .

 4 Fig. 4. Device utilization summary.

Fig. 5 .

 5 Fig. 5. Decision diagram of DIAGXTRM with tests and action applied.

  Histogram of execution times. Fig.6. Time trace and histogram of execution times (benchmark cnt on the GR-XC6S reference architecture without RTEMS).

  Reference architecture with software FPU (emulated).

  Architecture with hardware FPU: GRFPU-LITE. Fig. 8. Time traces of execution times for different Floating-Point architectures (benchmark lms on the ML605 without RTEMS).

Fig. 9 .

 9 Fig. 9. Measurements and model for pWCET estimation (benchmark matmult on the GR-XC6S without RTEMS).

  .Reference xc6s.elf.Direct xc6s.elf.LRR xc6s.elf.LRU (a) Different cache replacement policies (bench. matmult on the GR-XC6S without RTEMS).

  .GRFPU ml605.rtems.GRFPU (b) Impact of the RTOS: without and with RTEMS (benchmark lms on the ML605).

Fig. 10 .

 10 Fig. 10. Reliability of the pWCET estimations (scores out of 4.0 for each hypothesis).

  close the pWCET (Theoretical estimation) is to the maximum measured execution time (Max measured). As the models are all Weibull distributions, thus bounded to the right, they are quite close to the maximum measured values, Table

TABLE I

 I

		DEVICE UTILIZATION SUMMARY.
	GR-XC6S Occupied slices	ML605 Occupied slices
	architecture	out of 11662	architecture	out of 37680
	Reference	5219 ⇒ 44%	Reference	10667 ⇒ 28%
	Direct	4957 ⇒ 42%	GRFPU	20198 ⇒ 53%
	LRR	4791 ⇒ 41%	GRFPU-lite	13724 ⇒ 36%
	LRU	6895 ⇒ 59%
	noBP	4868 ⇒ 41%
	nocaches	4446 ⇒ 38%

TABLE III MEAN

 III EXECUTION TIME IN THOUSANDS OF CPU TICKS OF EACH BENCHMARK FOR DIFFERENT ARCHITECTURAL CONFIGURATIONS.

	Embedded system			GR-XC6S Architectures				ML605 Architectures
	benchmarking	Cache Replacement Policies		Key Elements	Best		FPU		Best
	Env.	Bench.	Reference	Direct	LRR	LRU	noBP	nocaches	Best1	Reference	GRFPU GRFPU-lite	Best2
		cnt	77.082	77.082	77.082	77.082	81.132	436.601	77.082	77.276	77.276	77.276	77.276
		matmult	139.292 142.168	149.393	132.077	154.869	1095.468	132.077	129.721	129.721	129.721	129.721
	ELF	nsichneu	12.624	18.337	29.122	21.806	13.517	110.600	12.624	21.638	21.791	21.785	12.627
		jfdctint	2.469	2.469	2.469	2.469	2.547	16.328	2.469	2.501	2.485	2.485	2.501
		lms	17.738	17.737	18.312	17.737	19.108	178.854	17.738	17.827	0.523	0.924	0.524
	RTEMS	jfdctint lms	2.496 17.980	2.498 18.032	2.503 18.304	2.495 17.977	2.575 19.387	18.268 207.048	2.496 17.979	2.530 18.078	2.532 0.535	2.532 0.939	2.532 0.535
	Global time scores 3	-	+8.32%	+24.58%	+11.67%	+6.62%	+663.63%	-1.01%	--19.35%	-18.91%	-27.71%

TABLE IV PWCET

 IV DISTRIBUTION PARAMETERS AND ESTIMATIONS FOR DIFFERENT BENCHMARKS AND ARCHITECTURES(a) Impact of the cache replacement policy (benchmark matmult on the GR-XC6S without RTEMS). Impact of the RTOS (benchmark lms on the ML605).

	GR-XC6S	Distribution	DIAGXTRM results
	architecture	shape ξ	Theoretical	Measured	Accuracy
	Reference	-0.249	139447252.58 139433822	9.63E-03
	Direct	-0.219	142170799.12 142170390	2.88E-04
	LRR	-0.095	149397501.05 149396595	6.06E-04
	LRU	-0.317	132076829.57 132076819	8.00E-06
	(b) ML605	Distribution	DIAGXTRM results
	with GRFPU	shape ξ	Theoretical	Measured	Accuracy
	ELF	-0.149	527244.18	526547	1.32E-01
	RTEMS	-0.596	542279.08	542228	9.42E-03