Ahlem Mifdaoui
email: ahlem.mifdaoui@isae.fr

Thierry Leydier
email: thierry.leydier@virtualitereelle.com

Beyond the Accuracy-Complexity Tradeoffs of Compositional Analyses using Network Calculus for Complex Networks

Keywords: Network Calculus, Pay Multiplex Only Once (PMOO), Performance analysis, Delay bounds, Accuracy-complexity trade-off, Total Flow Analysis (TFA)

Achieving the accuracy-complexity tradeoffs for compositional timing analyses using Network Calculus is still a hot research topic. In this specific area, we propose in this paper an improved version of the Total Flow Analysis (TFA) algorithm, called TFA++, when taking into account the impact of the finite transmission capacity of the network links on the input and output traffic models at each network node. First, we review the existing analysis algorithms by identifying their main limitations in terms of accuracy and complexity, through a simple but representative network example. Afterwards, we define the TFA++ algorithm and we detail the main steps of the followed methodology to compute the delay upper bounds. Moreover, we conduct comparative analyses of the derived delay bounds and analysis times with the different algorithms, with respect to the network size and load. In doing this, we highlight noticeable enhancements of both metrics under TFA++, in comparison to the existing algorithms; thus the high accuracy and low complexity of TFA++. Finally, this statement has been asserted through a representative avionics case.

INTRODUCTION

With the inherent complexity and the growing number of components, in addition to the use of advanced communication technologies and multi/many core processors, many challenges have emerged from designing new generation Cyber-Physical Systems (CPS). Hence, innovative cost-effective methods and tools are needed to design, analyze and verify the hardware and software architecture components of such systems, to guarantee predictability and reliability requirements while minimizing the development and maintenance costs. Particularly, for new generation CPS, the communication networks are considered as a bottleneck for performance and timing predictability; thus an appropriate performance analysis to provide safe guarantees has to be considered.

Most performance analysis approaches of complex networks can broadly be categorized under two main classes, simulation-based and analytical-based. The former enables the performance analysis of large-scale networks, but it is not appropriate to draw firm conclusions since it does not cover all corner cases, i.e., worst-case scenario; whereas, the latter is based on high-level analytical models that can be easily injected within optimization procedures, to speed up the design space exploration. These facts make the use of analytical-based approaches outwardly growing for early verification of networks in CPS, and one of the most relevant approaches in this specific area is the Network Calculus [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF]. The Network Calculus is a compositional algebraic framework to derive maximum bounds on system performance, such as delays and backlogs. The high modularity and scalability of such a framework make it particularly efficient to conduct performance analysis of complex communication networks [START_REF] Perathoner | Influence of Different Abstractions on the Performance Analysis of Distributed Hard Real-Time Systems[END_REF], such as Switched Ethernet [START_REF] Loeser | Low latency hard real-time communication over switched Ethernet[END_REF], the AFDX [START_REF] Grieu | Analyse et evaluation de techniques de commutation Ethernet pour l'interconnexion de systemes avioniques[END_REF], Networks on Chip [START_REF] Qian | Analysis of worst-case delay bounds for on-chip packet-switching networks[END_REF] and networks with cyclic dependencies [START_REF] Amari | Worst-case Timing Analysis of Ring Networks with Cyclic Dependencies using Network Calculus[END_REF].

However, one of the main challenging issues for Network Calculus is computing accurate performance bounds with a reasonable time-effort. Particularly, the lack of accuracy may lead to resource over-dimensioning of the networks under design; thus increasing development costs. Most of the related work in this area is focusing on improving the timing analysis algorithms [START_REF] Lenzini | A Methodology for Computing End-to-end Delay Bounds in FIFO-multiplexing Tandems[END_REF] [START_REF] Schmitt | Delay bounds under arbitrary multiplexing: When network calculus leaves you in the lurch[END_REF] [START_REF] Bouillard | Exact Worst-case Delay in FIFO-multiplexing Feedforward Networks[END_REF], to cope with the accuracy-complexity tradeoffs. Such research efforts lead to two main classes of approaches using Network Calculus: algebraic and optimization-based methods. In this paper, our main focus is the algebraic class, since it is the one keeping the composition property of Network Calculus. Particularly, there are mainly three existing timing analysis algorithms [START_REF] Schmitt | The DISCO Network Calculator -A Toolbox for Worst Case Analysis[END_REF]: Total Flow Analysis (TFA), Separated Flow analysis (SFA) and Pay Multiplex Only Once (PMOO). Each one of these aforementioned algorithms has its pros and cons in terms of accuracy and complexity, but there is no best algorithm in terms of both metrics.

To overcome these limitations, we propose in this paper an improved TFA algorithm (TFA++) to enhance accuracy, while keeping low complexity. The main idea of such an improved algorithm consists in considering the link transmission capacity impact on the input and output traffic models of each crossed node in the network.

Hence, the main contributions in this paper are as follows: (i) First, we review the existing analysis algorithms by identifying their main limitations in terms of accuracy and complexity, through a simple but representative network example. In doing this, we consolidate the idea of the non-existence of a best algorithm;

(ii) Second, we define the TFA++ algorithm (Algorithm 1) and we detail the main steps of the followed methodology to compute the delay upper bounds; (iii) Third, we conduct comparative analyses of the derived delay bounds and analysis times with the different algorithms, with respect to the network size and load. In doing this, we highlight noticeable enhancements of both metrics under TFA++, in comparison to the existing algorithms; thus the high accuracy and low complexity of TFA++. Finally, this statement has been asserted through a representative avionics case.

The rest of the paper is organized as follows: we start with presenting the main concepts of the Network Calculus framework in Section 2. Afterwards, we report the main related work and the identified limitations in Section 3. Then, we present the proposed TFA++ algorithm and report evaluation results in Sections 4 and 5, respectively. Finally, we draw the main conclusions and future work in Section 6.

THE NETWORK CALCULUS AS A COMPOSITIONAL ALGEBRAIC THEORY

The Network Calculus framework has been founded by the seminal work of Cruz in [START_REF] Cruz | A calculus for network delay. I. Network elements in isolation. Information Theory[END_REF][START_REF] Cruz | A calculus for network delay. II. Network Analysis. Information Theory[END_REF], and then extended with min-plus Algebra operations in [START_REF] Chang | Performance Guarantees in Communication Networks[END_REF] and [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF]. The latter extension is based on the idea of modeling the communication nodes as in conventional system theory, with an input function, a transfer function and an output function, where addition and multiplication are replaced by minimum and addition, respectively. Particularly, we will detail in this section how this algebraic extension has enabled the composition property of Network Calculus. We will answer herein some primordial questions when applying Network Calculus to conduct performance analysis of a realistic network: how to model the input traffic? how to model the node specifications? how to deal with a network of nodes to compute end-to-end performance?

Traffic Model

Network Calculus describes data flows by means of cumulative functions, R(t), defined as the number of transmitted bits during the time interval [0, t]. Function R(t) is always a non-decreasing function of time with R(0) = 0. Consider a system S receiving input data with cumulative function R(t), called input function. After the data processing, the output data is described using another cumulative function R * (t), called output function. The horizontal distance d(t) between the input and output functions represents the experienced delay of an input data received at time t in the system. The vertical distance x(t) between input and output functions represents the backlog, i.e., total number of bits present in the system at instant t.

One of the fundamental concepts in Network Calculus is the notion of maximum arrival curve, which provides an upper bound on the number of events, e.g., bits or packets, observed during any interval of time, as shown in Figure 1. Consequently, this curve covers all possible scenarios of traffic arrivals, observed at any instant, i.e., if α is a maximum arrival curve, then for any interval duration ∆, there will be at most α(∆) events. This concept allows modeling a large panel of event arrival patterns, such as periodic, sporadic, with or without jitter, and bursty or not. The definition of the arrival curve is as following: Definition 2.1. (Arrival Curve) [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF] A function α is an arrival curve for a data flow with an input cumulative function R, such that R(t) is the number of bits received until time t, iff: The considered arrival patterns necessary to define the arrival curve can be obtained from traffic traces if any, or formal specification. The latter is more common for real-time communication systems under design. The network designer generally specifies a traffic contract for each application, that is enforced using a controller, e.g., policers or shapers. The most common controller follows a leaky bucket algorithm, which guarantees for the controlled traffic a maximum burst b and a maximum rate r on the communication medium, i.e., the traffic flow is (b, r )-constrained. In this case, the arrival curve is an affine curve, defined as λ b,r (t) = b + r .t for t > 0. Furthermore, there is the following interesting result concerning the arrival curves. Lemma 2.2. [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF] If α 1 and α 2 are arrival curves for a flow R, then so is

∀t, s ≥ 0, s ≤ t, R(t) -R(s) ≤ α(t -s) (1)
α 1 ⊗ α 2 , where (f ⊗ д)(t) = inf 0≤s ≤t { f (t -s) + д(s)}.

Node Model

To conduct worst-case performance analysis, we need to put constraints on the input traffic through the maximum arrival curve notion. In return, we need to guarantee a minimum offered service within crossed nodes to cover the worst-case behavior and infer upper bounds on performance metrics, e.g., backlog and delay. This is done through the concept of minimum service curve, which is defined as following.

iff: ∀t ≥ 0, R * (t) ≥ inf s ≤t (R(t) + β(t -s)) (2)
A very useful and common model of service curve is the ratelatency curve β R,T , with R the minimum guaranteed rate and T the maximum latency before starting the service. This rate-latency function is defined as follows:

β R,T (t) = 0 if t ≤ T and R(t -T ) otherwise.
This service curve is easy to define in the case of one input/ output node serving one or many traffic flows coming from the same source and going to the same destination. However, to handle more realistic scenarios with a network of nodes, implementing aggregate scheduling which multiplexes the crossing flows at the input and demultiplexes them at the output, we need to define the left-over service curve guaranteed to each traffic flow within each crossed node, considering the impact of the other traffic flows in contention, to infer the offered guarantees for each flow. The computation of such a left-over service curve depends on the implemented scheduling policy within each crossed node, and the most common ones are Blind Multiplexing, FIFO and Fixed Priority (FP). It is worth noting that this derivation needs strict service curve property in the general case, except for FIFO and Constant bit rate nodes. A minimum strict service curve is defined as follows.

]s, t], R * (t) -R * (s) ≥ β(t -s).
The main result concerning the left over service curves computation is as follows: Theorem 2.5 (Residual service curve -Blind Multiplex). [START_REF] Bouillard | Service curves in Network Calculus: dos and don'ts[END_REF] let f 1 and f 2 be two flows crossing a server that offers a strict service curve β such that f 1 is α 1 -constrained, then the residual service curve offered to f 2 is:

β 2 = (β -α 1 ) ↑
where f ↑ (t) = max{0, sup 0≤s ≤t f (s)}

Performance Analysis

Knowing the arrival and service curves, one may compute the performance bounds for data flows, e.g., delay and backlog. In the case of single node with one input/ output, these bounds are computed according to the following theorem.

Theorem 2.6 (Performance Bounds). [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF] Consider a flow constrained by an arrival curve α crossing a system S that offers a minimum service curve β and a maximum service2 curve γ . The performance bounds obtained at any time t are given by: Output arrival curve 3

: α * (t) = α ⊘ β(t) Tight Output arrival curve: α * (t) = (γ ⊗ α) ⊘ β(t) Backlog 4 : ∀ t : q(t) ≤ v(α, β) Delay 5 : ∀ t : d(t) ≤ h(α, β)
The calculus of these bounds is greatly simplified in the case of leaky bucket arrival curve and the Rate-Latency service curve. In this case, the delay is bounded by b R + T , the backlog bound is b + r * T , and the output arrival curve is b + r (T + t). These results are illustrated in Figure 2.

To extend this result to a network of nodes, one of the strongest result in the Network Calculus framework is the computation of an end-to-end service curve for a tandem of nodes crossed by the 0 data (bits) time same flows. This curve is computed as the convolution of residual service curves in each node, computed according to Theorem 2.6. This result is described in the following theorem.

( ) t α ( ) t β ( , ) h α β ( , ) v α β
Theorem 2.7 (Concatenation-Pay Bursts Only Once).

[9] Assume a flow crossing two servers with respective service curves β 1 and β 2 . The system composed of the concatenation of the two servers offers a minimal service curve β 1 ⊗ β 2 to the flow.

Furthermore, this result infer an interesting property known as "Pay bursts Only Once Phenomena". Indeed, the end-to-end delay bound for a data flow, computed using the end-to-end service curve obtained with Theorem 2.7, clearly outperforms the sum of delay bound per node, computed iteratively using Theorem 2.6 and denoted as additive delay bound. The computation of these two bounds show the appearance of the burst term many times in the additive delay bound, as opposed to only once for the other.

CONVENTIONAL TIMING ANALYSIS ALGORITHMS AND LIMITATIONS

In the research community, there has been a growing interest in the subject of accuracy-complexity tradeoffs of performance analysis approaches using Network Calculus and we particularly focus herein on the compositional algebraic approaches. There are three main algorithms to compute the end-to-end delay bound for each flow of interest (f.o.i). All of them need to start the procedure by defining the initial input arrival curve within the network, i.e., within the source nodes, and the service curve guaranteed within each crossed node to the aggregate flows. Then, the main difference between the different algorithms relies on how to use these traffic and node models to derive the end-to-end delay bounds.

Total Flow Analysis (TFA) [START_REF] Schmitt | The DISCO Network Calculator -A Toolbox for Worst Case Analysis[END_REF] This algorithm consists in computing iteratively the delay upper bound in each crossed node for the aggregate flow, then the sum results in end-to-end delay bounds [START_REF] Cruz | A calculus for network delay. II. Network Analysis. Information Theory[END_REF]. The main steps are: (i) summing up all the arrival curves of individual flows at the input of each crossed node; (ii) computing the delay bound based on Theorem 2.6, when considering the arrival curve of the aggregate flows and the global service curve guaranteed within the node;

(iii) computing the output arrival curve of each flow using Theorem 2.6, when considering the left-over service curve based on Theorem 2.5.

Separated Flow Analysis (SFA) [START_REF] Schmitt | The DISCO Network Calculator -A Toolbox for Worst Case Analysis[END_REF] This algorithm consists in considering the end-to-end service curve, through the concatenation of the left over service curves, guaranteed for the f.o.i within each crossed node (applying Theorem2.7). This method is also based on three steps: (i) summing up the input arrival curves of flows in contention with the f.o.i; (ii) computing the left over service curve guaranteed to the f.o.i, based on Theorem 2.5; (iii) computing the output arrival curve of the f.o.i.

These steps need to be conducted for each flow within each crossed node, to infer all the input arrival curves for the next node. This procedure has to be continued along each flow path to have all the left over service curves; thus the end-to-end service curve.

Pay Multiplex Only Once (PMOO) [START_REF] Schmitt | Improving Performance Bounds in Feed-Forward Networks by Paying Multiplexing Only Once[END_REF] The main idea of this algorithm is based on accounting the flow serialization phenomena along the flow path to compute tighter end-to-end delay bounds. This fact consists in defining a smart application order of Theorem 2.7 and Theorem 2.5 to start concatenating as much as possible the service curves of nodes, crossed by the same group of interfering flows with the f.o.i (applying The-orem2.7); and then to compute the left-over service curve of the concatenated system guaranteed to the f.o.i (applying Theorem 2.5).

In doing this, we pay only once the bursts of interfering flows with the f.o.i along its path.

Limitations and Discussion

To better evaluate these algorithms, we conduct a brief comparative analysis of derived delay bounds and analysis times under the different algorithms, when increasing the network load for a simple but representative network example, shown in Figure 3. The derived delay bounds and analysis times under the different algorithms regarding the network load are illustrated in Figure 4 and Table 1, respectively. As we can notice, the delay bounds increase linearly with the network load. Furthermore, the PMOO guarantees the lowest delay bounds, which confirms the high accuracy of PMOO in comparison with TFA and SFA algorithms. On the other hand, the TFA offers the lowest maximum analysis time per flow, which shows the low complexity of TFA in comparison with SFA and PMOO. These results consolidate our statement that there is no best algorithm in terms of both metrics.

TFA SFA PMOO Max. Analysis Time per flow (us) 897 1683 953 Table 1: Maximum Analysis Times per flow of the different algorithms for the Network example Therefore, there is a need of a new solution to bridge the gap between these existing algorithms and guarantee high accuracy and low complexity. To achieve this aim, we make the choice of improving the TFA algorithm to favor the time-effort metric, and this extension is detailed in the next section.

IMPROVED TFA ALGORITHM

In this section, we present the main idea of our proposed algorithm, called TFA++, to compute accurate delay bounds with low complexity. Afterwards, we detail the main steps of the computation methodology and explain the algorithm. Finally, we express the delay bounds for an illustrative network example under TFA++ and PMOO, to show the expected enhancement in terms of accuracy, while keeping a low complexity.

Main Idea

Without loss of generality, we model the network as a direct graph where the nodes are the output ports of the different network components and the edges the network links with finite transmission capacities. It is worth noting that the output ports are generally the only network parts inducing unknown multiplexing delays; whereas the other kinds of ports or blocks generally infer maximum constant delays that we can easily add at the end of the delay bound computation.

The main idea of the TFA++ is to take into account the impact of the link capacity on the traffic model to enhance the end-to-end delay bound tightness. Hence, we present herein two improvements due to this fact.

Tighter input arrival curve of an aggregate traffic

To refine the input arrival curve of an aggregate traffic at a given node with multiple inputs, where each input l has a finite transmission capacity C l , we define the following notations:

• i ∋ k the set of flows crossing the node k;

• Each flow i has an input arrival curve at each crossed node k, α k ⊖1 i (t); • Pred N is the set of nodes transmitting flows on the input of node N .

The aggregate traffic sent by each node l ∈ Pred N on the input of node N can be modeled through two possible arrival curves:

(1) the first one is simply the sum of individual flow arrival curves crossing l, i ∋l α N ⊖1 i ; (2) knowing the finite transmission capacity of the link from node l to node N , C l , we have a second arrival curve of this same aggregate traffic, λ 0,C l .

Hence, based on Lemma 2.2, the aggregate traffic flow at the input of node N has the following arrival curve:

l ∈P r ed N (λ 0,C l ⊗ i ∋l α N ⊖1 i (t)) (3)
According to Theorem 3.1.6 in [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF], for functions passing through the origin, the convolution of both functions is at most equal to their minimum. Moreover, for the particular case of concave functions, e.g., leaky bucket curves, the convolution is equal to the minimum.

Hence, the refined input arrival curve in Eq. 3 used in TFA++ is necessarily tighter than the sum of the individual arrival curves, used in the original TFA algorithm:

l ∈P r ed N (λ 0,C l ⊗ i ∋l α N ⊖1 i (t)) ≤ l ∈P r ed N i ∋l α N ⊖1 i (t) (4)
Tighter output arrival curve of an individual flow The finite transmission capacity at the output of a node N , C N , infers a maximum service curve γ N (t) = C N .t.

Hence, through applying Theorem 2.6, we can compute the following output arrival curve from node N for each flow i, α N i , as follows:

α N i (t) = (α N ⊖1 i ⊗ γ N ) ⊘ β N i (t) (5)
where β N i (t) is the left-over service curve of the flow i at node N computed through applying Theorem 2.5.

The refined output arrival curve in Eq. ( 5) used in TFA++ is tighter than the one used in the original TFA algorithm:

(α N ⊖1 i ⊗ γ N ) ⊘ β N i (t) ≤ α N ⊖1 i ⊘ β N i (t) (6)
In the next section, we explain the use of both improvements to compute the end-to-end delay bounds through the Algorithm 1.

Computation Methodology

Algorithm 1 Improved Total Flow Analysis (N etwork, F )

1: for each flow of interest f ∈ F do 2:
for each sink ∈ sinks(f ) do 3: for each flow k ∈ K do 12:

D EED (f , sink) ← 0 4: for each node N ∈ path sr c(f )=>sink (f ) do 5: MaxDelay ← 0 6: K ← set-of-flows-at-input (N , f ) 7: Pred ← Predecessors(N ) 8: MaxInput ← ComputeInput(N , f
α N k (t) ← ComputeOutput(N , f ) 13:
end for 14:

D EED (f , sink) ← D EED (f , sink) + D 15:
end for 16:

D EED (f ) ← Vector(D EED (f , sink)) 17:
end for 18:

D EED ← Matrix(D EED (f ))
19: end for 20: return D EED

The Delay bounds analysis based on the proposed TFA++ is described in Algorithm 1 for feed-forward networks, i.e., no cyclic dependencies, and supporting multicast communications.

For each flow of interest f (line 1) and for each one of its paths path sr c(f )=>sink (f ) (line 2), the delay bound is computed within each crossed node N , accounting the set of flows in contention with f , K (line 6) and the the set of predecessors (line 7). Finally, the end-to-end delay is simply equal to the sum of these individual delays along the flow path (lines 4-15).

First, we compute the input arrival curve of the node N , considering the impact of the link capacity, as explained in Eq. (3) (line 8). Then, we compute the total service curve of the crossed node (line 9), to enable the computation of the delay bound within the crossed node through applying Theorem 2.6 (line 10). Afterwards, for each flow in the interference set, we need to compute the output arrival curve when taking into account the impact of the link capacity, as explained in (Eq. 5) (lines 11-13). Finally, the end-to-end delay bound is the sum of all the crossed node delays until reaching the corresponding sink (line 16). We consider herein the illustrative example in Figure 5 to express the derived delay bounds under TFA++ and PMOO for the flow of interest f 1 with the arrival curve α 1 . We detail herein the delay bounds computation when considering identical flows, i.e., f 1 and f 2 , with leaky-bucket arrival curves, λ b,r , and nodes with ratelatency service curves, β C,T .

Illustrative Example

Under PMOO, as we can notice, both flows are crossing the same nodes 3 and 4 until the final destination. Hence, we can start by concatenating the nodes 3 and 4 through applying Theorem 2.7, then computing the residual service curve offered to the f.o.i f 1 through Theorem 2.5. This combination gives the following end-toend service curve under PMOO for f 1 :

β 1 ⊗ (β 3 ⊗ β 4 -α 3⊖1 2 ) ↑ = β C-r,T +(2T .C+b+r .T )/C-r (7)
Consequently, the end-to-end delay bound under PMOO when applying Theorem 2.6 is D P MOO = T + (2T .C + 2b + r .T )/(C -r ).

Under TFA++, we compute the end-to-end delay bound as the sum of local delays in nodes 1, 3 and 4 crossed by the f.o.i. f 1 . First, the delay bound in node 1 is simply equal to the maximum horizontal distance between the arrival and service curves; thus As we can notice, D T F A++ < D P MOO for this illustrative example. This fact shows the promising accuracy of TFA++ in comparison to PMOO, and we will consolidate this statement in the next section under different network configurations.

EVALUATION

In this section, we first report the computed delay bounds and time analysis with respect to network size and load, to benchmark the existing timing analysis algorithms against our proposed TFA++ in terms of accuracy and complexity. Afterwards, we present the delay bounds for a representative avionics study to consolidate our evaluation. All the computed metrics in this section are based on the WoPANets tool [START_REF]WoPANets: Worst-Case Performance Analysis of embedded Networks[END_REF].

Comparative Analyses

To conduct the comparative analyses, we consider the case study with the following assumptions: (i) The topology is a mesh connecting M nodes, as the one illustrated in Figure 6; (ii) All nodes guarantee a rate-latency service curve β R,T with R = 100Mb/s and T = 40us;

(iii) Each node generates leaky-bucket constrained flows with a burst σ and a rate ρ; (iv) The network has a maximum bottleneck utilisation rate U .

Furthermore, we define various network configurations, where each network configuration is defined with the tuple (σ , ρ, M, U ). We vary the network size M or load U of this tuple at a time, to highlight its impact on the derived delay bounds and time analysis. As we can notice, TFA++ implies the lowest delay bounds for both scenarios, in comparison to existing algorithms. This fact shows its high accuracy. Moreover, when increasing the network size, i.e., the number of nodes, or the network load, i.e., the number of flows, after a given threshold value, PMOO leads to worse delay bounds than SFA, since it becomes more and more difficult to find a smart application order of Theorem 2.7 and Theorem 2.5 under these conditions; thus a more pessimistic bounds.

On the other hand, Figures 9 and 10 show the impact of the network size and load on the time analysis under the different timing algorithms, respectively. As illustrated, TFA++ infers almost the same time analysis than the original TFA, which still is the lowest one in comparison to existing algorithms. In addition, we can notice that the PMOO analysis time becomes the highest when increasing the network size due to the high call number of convolution operator (⊗); whereas when increasing the number of flows, the SFA leads to the highest analysis time due to the high call number of deconvolution operator (⊘).

Hence, these results confirm our first conclusions in Section 3 concerning the pros and cons of the existing algorithms, where the PMOO is considered in general as the most accurate one but also the most complex one. Moreover, they show the performance of our proposed algorithm in terms of high accuracy and low complexity, in comparison to existing solutions.

Avionics Case Study

We report herein the derived delay bounds under the different timing algorithms for a representative AFDX network of A350, illustrated in Figure 11. This network consists of 7 switches and more than 60 end-systems. The latter generate more than 1100 multicast Virtual Links, inferring more than 8600 flows. The maximum 12 shows the delay bounds for the different flows following the ascending order under the different timing algorithms. As we can notice, both TFA++ and PMOO imply the most accurate delay bounds for most of the flows. However, for this realistic case study, there is no strict order relation between both algorithms. Hence, to derive the tightest delay bounds, we can consider the minimum of the computed delay bounds for each flow using PMOO and TFA++. It is worth noting that the TFA++ still has the lowest time analysis for this avionics case study. This fact shows that our proposal still is a good value for this realistic case study, in comparison to the existing algorithms in terms of accuracy and complexity.

CONCLUSIONS

A new timing analysis algorithm based on Network Calculus, called TFA++, has been proposed in this paper to cope with the known accuracy-complexity tradeoffs in the literature. The main idea of TFA++ is to take into account the impact of the finite transmission capacity of the links on the input and output traffic models of each crossed node. Moreover, this algorithm outperforms the main existing solutions, in terms of high accuracy and low complexity for some network configurations.

However, the avionics case has shown that there is no strict order relation between PMOO and TFA in terms of accuracy. This fact encourages us to lead further investigations on the comparative analysis of both algorithms, when considering various network topologies and flow parameters, i.e., the burst, the rate or the path length.

Figure 1 :

 1 Figure 1: Arrival curve

Definition 2 . 3 .

 23 [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF] (Simple Minimum Service Curve) The function β is the simple service curve for a data flow with an input cumulative function R and output cumulative function R *

Definition 2 . 4 .

 24 [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF] (Strict service curve) The function β is the strict service curve for a data flow with an input cumulative function R and output cumulative function R * , if for any backlogged period1

Figure 2 :

 2 Figure 2: Backlog and Delay Bounds

Figure 3 :

 3 Figure 3: Illustrative Network Example

Figure 4 :

 4 Figure 4: Delay bounds vs Network load of the different algorithms for the Network example

  )

Figure 5 :

 5 Figure 5: Illustrative Example

D 1 =

 1 T +b/C. Afterwards, to compute the delay in node 3, we define the input arrival curve of node 3 as in Eq. (3), 2 * min(C.t, b +r (t +T )), and consequently D 3 = T + b/(Cr ). Finally, we compute the output arrival curve of node 3 as in Eq. (5), min(C.t, α 4⊖1 1 + α 4⊖1 2 ); thus D 4 = T . Hence, the end-to-end delay bound under TFA++ is D T F A++ = 3T + b/(Cr ) + b/C.

Figure 6 :

 6 Figure 6: Considered Network topology for the comparative analysis

Figures 7

 7 Figures7 and 8show the impact of the network size and load on the delay bounds under the different timing algorithms, respectively. As we can notice, TFA++ implies the lowest delay bounds for both scenarios, in comparison to existing algorithms. This fact shows its high accuracy. Moreover, when increasing the network size, i.e., the number of nodes, or the network load, i.e., the number of flows, after a given threshold value, PMOO leads to worse delay bounds than SFA, since it becomes more and more difficult to find a smart application order of Theorem 2.7 and Theorem 2.5 under these conditions; thus a more pessimistic bounds.On the other hand, Figures9 and 10show the impact of the network size and load on the time analysis under the different timing algorithms, respectively. As illustrated, TFA++ infers almost the same time analysis than the original TFA, which still is the lowest one in comparison to existing algorithms.

Figure 7 :Figure 8 :

 78 Figure 7: The impact of the network size on the delay bounds for (σ = 400bytes, ρ = 20kbps, M ∈ [2x2 : 1x1 : 16x16], U = 30%)

Figure 9 :Figure 10 :

 910 Figure 9: The impact of the network size on the analysis time for (σ = 400bytes, ρ = 20kbps, M ∈ [2x2 : 1x1 : 16x16], U = 30%)

Figure 11 :

 11 Figure 11: Representative Avionics Case Study

Figure

  Figure12shows the delay bounds for the different flows following the ascending order under the different timing algorithms. As we can notice, both TFA++ and PMOO imply the most accurate delay bounds for most of the flows. However, for this realistic case study, there is no strict order relation between both algorithms. Hence, to derive the tightest delay bounds, we can consider the minimum of the computed delay bounds for each flow using PMOO and TFA++. It is worth noting that the TFA++ still has the lowest time analysis for this avionics case study. This fact shows that our proposal still is a good value for this realistic case study, in comparison to the existing algorithms in terms of accuracy and complexity.

Figure 12 :

 12 Figure 12: Delay bounds for the avionics case study under the different analysis algorithms

A backlogged period ]s, t ] is an interval of time during which the backlog is non null, i.e., R(s) = R * (s) and ∀u ∈]s, t ], R(u) -R * (u) > 0

∀t, R * (t ) ≤ R ⊗ γ (t )

f ⊘ д(t ) = sup s ≥0 {f (t + s) -д(s)}

v(f , д) is the maximum vertical distance between f and д

h(f , д) is the maximum horizontal distance between f and д


Emmanuel Hermellin
email: hermellin@lirmm.fr

Fabien Michel
email: fmichel@lirmm.fr

GPU Delegation: Toward a Generic Approach for Developping MABS using GPU Programming

Keywords: MABS, GPGPU, CUDA, Hybrid architecture

 in 2008 was twofold: (1) data parallel execution capabilities of GPU can be used effectively in MABS and afford excellent speedup on models and (2) effective use of data parallel execution requires resolution of modeling and execution challenges at the cost of a decrease in modularity, ease of programmability and reusability. In this paper, we propose to study through experiments if the conclusions and issues outlined by Perumalla and Aaby are still true despite the evolution of GPGPU and MABS. To this end, we use the GPU environmental delegation principle on four models in order to compare CPU and GPU implementations. Then, we discuss and analyze the results from both a conceptual and a performance point of view.

INTRODUCTION

Multi-Agent Based Simulations (MABS) are used to study complex systems in many research domains (collective robotics, biology, economy, urban management, etc.) [START_REF] Michel | Multi-Agent Systems and Simulation: a Survey From the Agents Community's Perspective[END_REF]. Experimenting with MABS, performances often represent a major issue limiting the extent to which an Agent-Based Model (ABM) could be studied [START_REF] Parry | Large scale agent-based modelling: A review and guidelines for model scaling[END_REF]. Especially, it is interesting to not be too limited by the size of the environment or the number of agents, which may require a lot of computing resources.

That is why High Performance Computing (HPC) is more and more considered as a relevant way of speeding up MABS thanks to functional solutions that use Central Processing Appears in: Proceedings of the 15th International Conference Unit (CPU) clusters, computational grids, etc. A flagship example of this trend is the Repast HPC platform1 [START_REF] Collier | Repast HPC: A Platform for Large-Scale Agent-Based Modeling[END_REF], which allows to simulate millions of agents using CPU clusters.

Among HPC technologies, Graphics Processing Units (GPU) have been proved to be excellent computational platforms able to perform general-purpose computations [START_REF] Owens | A Survey of General-Purpose Computation on Graphics Hardware[END_REF]. General-Purpose computing on Graphics Processing Units (GPGPU) is thus a technology which relies on using the massively parallel architecture of usual PC graphics cards for accelerating very significantly the performance of programs (applications can be up to 100 times faster)2 [START_REF] Che | A performance study of general-purpose applications on graphics processors using CUDA[END_REF]. So, GPGPU is a cheap HPC solution which could be now used on almost every PC.

In 2008, Perumalla and Aaby's work [START_REF] Perumalla | Data parallel execution challenges and runtime performance of agent simulations on GPUs[END_REF] (P&A for short) was one of the first that studied the use of GPGPU in the context of MABS. The purpose of this work was to evaluate the advantages and drawbacks of this technology by implementing several representative ABM. So, P&A compared the runtime speed of GPU-based models against CPU-based implementations and showed that impressive speedup could be achieved (up to x40 in the paper). However, one major conclusion of this work was that this increase of performance could only be done at the expense of programmability, modularity and reusability. Indeed, it turns out that implementing ABM using GPGPU is very challenging because GPU programming relies on a highly specialized hardware architecture. Especially, usual object oriented features, which are common in ABM, are no longer available using GPGPU.

Since this study, graphic cards and their programming interfaces have dramatically evolved so that GPGPU is now recognized as a very efficient HPC solution in many domains. Still, considering GPGPU for MABS, it is worth noting that, despite some important milestones, very few works have been done: Implementing ABM using GPGPU remains very challenging and P&A's conclusions still hold [START_REF] Lysenko | A Framework for Megascale Agent Based Model Simulations on Graphics Processing Units[END_REF][START_REF] Erra | An Efficient GPU Implementation for Large Scale Individual-Based Simulation of Collective Behavior[END_REF][START_REF] Richmond | A High Performance Framework For Agent Based Pedestrian Dynamics On GPU Hardware[END_REF][START_REF] Aaby | Efficient Simulation of Agent-based Models on multi-GPU and Multi-core Clusters[END_REF]. Moreover, research works that only focus on performances are too hard to reuse, which does not ease the adoption of GPGPU in the MABS community.

Among the few works dealing with accessibility issues and that do not only focus on performance gains, almost all seek to hide GPGPU through a transparent use of this technology [START_REF] Richmond | Agent based GPU, a real-time 3d simulation and interactive visualisation framework for massive agent based modelling on the GPU[END_REF][START_REF] Richmond | High performance cellular level agent-based simulation with FLAME for the GPU[END_REF][START_REF] Laville | MCMAS: A Toolkit to Benefit from Many-Core Architecure in Agent-Based Simulation[END_REF]. However, this approach cannot take into account all the cases and needs that the implementations of MABS with GPGPU could require. Indeed, it is hard to create a framework or universal solution because of the wide variety of models.

In this paper, we present a completely different approach to these works: Instead of focusing on accessibility by means of transparency, we show that it is more interesting to use a design principle which transforms or adapts a model because it is thus possible to take into account a wider variety of models. To this end, we use GPU environmental delegation [START_REF] Michel | Translating Agent Perception Computations into Environmental Processes in Multi-Agent-Based Simulations: A means for Integrating Graphics Processing Unit Programming within Usual Agent-Based Simulation Platforms[END_REF] as a way to implement MABS with GPGPU. Based on an hybrid approach (the execution is shared between the CPU and the GPU), this principle put forward reusability rather than transparency [START_REF] Hermellin | GPU Environmental Delegation of Agent Perceptions: Application to Reynolds's Boids[END_REF]. The idea is to identify in the model some computations which can be transformed into environmental dynamics and then translated into GPU modules. The purpose of this article is to study the interest of this principle and show that it helps to address some of the problems outlined by P&A. This paper is organized as follows: Section 2 introduces the evolution of the GPGPU technology in the MABS domain and presents the GPU environmental delegation. Section 3 presents the models used for the experimentation and describes the adaptation process and implementations of these models using the GPU environmental delegation principle. Section 4 proposes a performance study and analyses the results of the experimentations while answering the questions expressed in this paper. Section 5 summarizes the work we have done and outlines our future researches.

GPGPU FOR MABS

Related Works

In 2008, the implementation made by P&A used graphical features of the card to perform generic computations. Indeed, at this time, using a specialized programming interface dedicated to GPU programming was not widespread, due to the youth of these tools. So, the few people interested by this technology had to divert the primary use of GPU by tinkering with graphical functions [START_REF] D'souza | SugarScape on steroids: simulating over a million agents at interactive rates[END_REF][START_REF] Richmond | Agent based GPU, a real-time 3d simulation and interactive visualisation framework for massive agent based modelling on the GPU[END_REF].

However, the release of two programming interface dedicated to GPGPU, CUDA3 (Compute Unified Device Architecture, 2007) and OpenCL 4 (Open Computing Language, 2008), have greatly simplified this technology while offering better support and improving its accessibility. GPGPU thus became a very appealing alternative for many domains in which computation time is critical. Consequently, the number of MABS using GPGPU has increased and new tools and frameworks have emerged.

Flame GPU [START_REF] Richmond | High performance cellular level agent-based simulation with FLAME for the GPU[END_REF] is a flagship example of the possibilities offered by the rise of specialized GPU programming tools. However, the majority of the existing works still start from scratch and put all their attention on acquiring the best computational gains without considering the accessibility, reusability and modularity. This is especially true in the scope of works that address the study of flocking models [START_REF] Erra | An Efficient GPU Implementation for Large Scale Individual-Based Simulation of Collective Behavior[END_REF][START_REF] Husselmann | Simulating Species Interactions and Complex Emergence in Multiple Flocks of Boids with GPUs[END_REF], crowd [START_REF] Richmond | A High Performance Framework For Agent Based Pedestrian Dynamics On GPU Hardware[END_REF], traffic simulations [START_REF] Strippgen | Multi-agent traffic simulation with CUDA[END_REF] or autonomous navigation and path planning algorithms [START_REF] Bleiweiss | Multi agent navigation on the GPU[END_REF][START_REF] Demeulemeester | Hybrid Path Planning for Massive Crowd Simulation on the GPU[END_REF].

However, as pointed out in [START_REF] Aaby | Efficient Simulation of Agent-based Models on multi-GPU and Multi-core Clusters[END_REF], implementing a model on GPU does not necessarily imply an increase of performances, notably in the field of MABS where many different and heterogeneous architectures can be conceived. Considering this issue, [START_REF] Aaby | Efficient Simulation of Agent-based Models on multi-GPU and Multi-core Clusters[END_REF] clearly shows the influence of implementation choices on the results and performances. So, despite the release of frameworks such as CUDA, achieving an effi-cient implementation still requires to take into account the specific programming model that comes with GPU. That is why, most of the previous cited works were realized in an ad hoc way and only offer one-off solutions.

Before 2011, most of works used an approach that consist in executing completely the model on the GPU (all-in-GPU ). This approach is useful when the main objective is to accelerate the simulation's execution. But from a software engineering point of view, this approach is not adapted because all development efforts are lost. In such a case, the implementation is designed for a specific model and therefore cannot be reused in other contexts. There is no possibility to generalize and reuse the work.

Considering these issues, hybrid systems begin to appear in 2011 and represent a very attractive alternative because they consist in sharing the execution of the MABS between the CPU and the GPU according to the nature of the computations. Despite the fact that an all-in-GPU implementation is more efficient than an hybrid implementation, the hybrid approach has many advantages.

Firstly, hybrid approaches enable a step further toward the achievement of more complex MAS models because one can choose what is executed on the GPU. It is thus possible to use agents with complex and heterogeneous architectures (e.g. [START_REF] Laville | Using GPU for Multi-agent Multi-scale Simulations[END_REF]). Moreover, they are inclined toward modularity, and for instance ease to concretely implement independently the agents and the environment (e.g. [18] [23]).

So, even if genericness is not necessarily an explicit objective of hybrid systems, such approaches actualy exhibit a software architecture promoting modularity and thus reusability. The MCMAS library [START_REF] Laville | MCMAS: A Toolkit to Benefit from Many-Core Architecure in Agent-Based Simulation[END_REF] is a recent and good example of this trend. Moreover, by removing the programming constraints related to all-in-GPU systems, hybrid approaches are by definition more flexible and open to other technologies so that accessibility is greatly increased [START_REF] Laville | Using GPU for Multi-agent Multi-scale Simulations[END_REF][START_REF] Michel | Translating Agent Perception Computations into Environmental Processes in Multi-Agent-Based Simulations: A means for Integrating Graphics Processing Unit Programming within Usual Agent-Based Simulation Platforms[END_REF].

Analyzing the Evolution of MABS and GPGPU

In their work, P&A concluded that the data parallel execution capabilities of GPU are very useful to speed up simulation of agent-based models. The empirical results from their experiments show that high simulation speeds can be achieved especially with very large agent populations (million of agents). The runtime speeds on the GPU are roughly two to three orders of magnitude higher than those of existing ABM toolkits. However, the excellent speedup obtained on simple models comes at the expense of modularity, ease of programmability and reusability.

From this overview, it is clear that the conclusions set out by P&A have been taken into account to some extent. But, research works are still divided into two distinct categories:

(1) Works that are only interested in performance gains and (2) works that try to answer the issues raised by P&A about modularity, ease of programmability and reusability. The works of this second category are mostly based on hybrid systems, which is clearly the most promising approach to address such issues. However, these works propose solutions that only focus on a transparent use of GPGPU in MABS. But, because of the wide variety of MAS models, making the use of GPGPU transparent cannot be generic enough and take into account all the cases and needs that MABS implementations with GPGPU could require.

So, the stated problems were not clearly resolved and the lack of reusability, accessibility and modularity are still present (despite work such as Flame GPU [START_REF] Richmond | High performance cellular level agent-based simulation with FLAME for the GPU[END_REF] or [START_REF] Laville | Using GPU for Multi-agent Multi-scale Simulations[END_REF][START_REF] Laville | MCMAS: A Toolkit to Benefit from Many-Core Architecure in Agent-Based Simulation[END_REF]). This is especially true considering that the proportion of works dealing with the use of GPGPU in MABS is still very small compared to the number of MABS-based publications and while the need in computing resources is strong and always increasing.

GPU Environmental Delegation

The purpose of the work presented in this article is to determine whether the conclusions stated by P&A in 2008 are still relevant today despite the evolution of GPGPU and its use in MABS. To this end, the GPU Environmental Delegation principle (GPU delegation for short) is used. Based on an hybrid approach, this design principle aims at easing the integration of GPGPU in MABS models rather than making its use transparent, in contrary to the works previously mentioned.

Overview of the Principle

GPU delegation is inspired by an Agent-Oriented Software Engineering (AOSE) trend which consists in using the environment as a first class abstraction in MAS [START_REF]4th International Workshop, E4MAS 2014 -10 Years Later[END_REF]. This idea is today well accepted and has proved to be a relevant approach for modeling and developing MAS [START_REF] Weyns | Environments for Multiagent Systems State-of-the-Art and Research Challenges[END_REF]. Especially, it could help to enhance the efficiency of agent interactions or simplify the behavioral process of the agents thanks to the environment that produce high level percepts. Therefore, at a high level, such an approach allows to design MAS with a clear separation of concerns.

So, at the implementation level, GPU delegation has to be related to other research works that consider the environment as a core concept of MAS and which reify parts of the agents' computations in external structures. Related examples are EASS (Environment As Active Support for Simulation) [START_REF] Badeig | A contextual environment approach for multi-agent-based simulation[END_REF], IODA (Interaction Oriented Design of Agent simulations) [START_REF] Kubera | IODA: an interaction-oriented approach for multi-agent based simulations[END_REF], environment-centered approach for MABS [START_REF] Payet | Environment as support for simplification, reuse and integration of processes in spatial MAS[END_REF] and the artifact approach [START_REF] Ricci | Environment programming in multi-agent systems: an artifact-based perspective[END_REF]. In the scope of our work, considering the environment as a first class entity enables us to keep the programming accessibility of the agent model in a GPU context and be able to scale up both the number of agents and the size of the environment.

Proposed in [START_REF] Michel | Translating Agent Perception Computations into Environmental Processes in Multi-Agent-Based Simulations: A means for Integrating Graphics Processing Unit Programming within Usual Agent-Based Simulation Platforms[END_REF], the GPU delegation principle is based on the fact that it is very difficult to deport all the behavior of agents on graphics cards (that is why it is especially reactive agents who evolve in simulations which use an all-in-GPU approach [START_REF] Richmond | Agent based GPU, a real-time 3d simulation and interactive visualisation framework for massive agent based modelling on the GPU[END_REF]). So, this principle consists in making a clear separation between the agent behaviors, managed by the CPU, and environmental dynamics, handled by the GPU. Especially, one major idea underlying this principle is to identify some computations (such as agent-level perceptions) which can be transformed into environmental dynamics.

GPU Delegation in Practice

GPU delegation has been used for the first time on a model of Multi-Level Emergence (MLE) [START_REF] Beurier | Model and Simulation of Multi-Level Emergence[END_REF] of complex structures [START_REF] Michel | Translating Agent Perception Computations into Environmental Processes in Multi-Agent-Based Simulations: A means for Integrating Graphics Processing Unit Programming within Usual Agent-Based Simulation Platforms[END_REF]. This very simple model relies on a unique behavior which allows to generate complex structures which repeat in a fractal way. The agent behavior is extremely simple and based on the perception, spread and reaction to pheromones. So, in this work, GPU modules dedicated to the perception and the spread of pheromones were proposed.

A second experimentation was proposed in [START_REF] Hermellin | GPU Environmental Delegation of Agent Perceptions: Application to Reynolds's Boids[END_REF] which further trials GPU Delegation by testing its feasibility and genericness on a classic ABM, namely Reynolds's boids5 . [START_REF] Hermellin | GPU Environmental Delegation of Agent Perceptions: Application to Reynolds's Boids[END_REF] especially shows that applying GPU delegation not only speeds up boids simulations but also produces an ABM which is easy to understand, thanks to a clear separation of concerns.

For both of these works, the integration of GPU computations was performed in the TurtleKit platform6 [START_REF] Michel | The TurtleKit Simulation Platform: Application to Complex Systems[END_REF]. TurtleKit is a generic spatial ABM, implemented with Java, wherein agents evolve in a 2D environment discretized in cells. The proposed hybrid approach integrated in TurtleKit focuses on modularity. In this context, this allows to achieve three objectives: (1) maintaining accessibility in the agent model while using GPGPU, (2) being able to scale and work with large number of agents on large environment sizes and (3) promoting reusability in the particular context of GPU programming.

Evolution of GPU Delegation

With respect to the underlying hybrid approach, GPU Delegation is about identifying specific behaviors which can be turned into environmental dynamics. Especially, GPU delegation states that agent perception computations that do not involve the agent's states could be translated into environmental dynamics and thus into a GPU module performing the computation.

However, the flocking model experimentation [START_REF] Hermellin | GPU Environmental Delegation of Agent Perceptions: Application to Reynolds's Boids[END_REF] shows that, according to the model used, GPU delegation could be extended with the aim at finding more computations to be delegated. So, the GPU delegation principle has evolved in order to be applied on a larger number of models and could be now stated as follows: If computations made within the behavior of the agent do not involve or do not modify the agent's states, they could be translated into environmental dynamics.

MODELS AND IMPLEMENTATIONS

To provide a comparison between the study conducted by P&A in 2008 and solutions available in 2015, especially hybrid systems, we use GPU Delegation to adapt and implement four models: Two already done by P&A (Conway's Game of Life and Schelling's segregation) and two taken from the Netlogo models library [START_REF] Sklar | NetLogo, a Multi-agent Simulation Environment[END_REF] (Fire model and DLA model). In this section, we present the associated technology and describe the implementation's process.

GPGPU Implementation with CUDA

To program on the graphics card and exploit its GPGPU capabilities, we use CUDA which is the GPGPU programming interface provided by Nvidia. The associated programming model relies on the following philosophy 7 : The CPU is called the host and plays the role of scheduler. The host manages data and triggers kernels, which are functions specifically designed to be executed by the GPU, which is called the device. The GPU part of the code really differs from sequential code and has to fit the underlying hardware architecture. More precisely, the GPU device is programmed to proceed the parallel execution of the same procedure, the kernel, by means of numerous threads. These threads are organized in blocks (the parameters blockDim.x, blockDim.y characterize the size of these blocks), which are themselves structured in a global grid of blocks. Each thread has unique 3D coordinates (threadIdx.x, threadIdx.y, threadIdx.z ) that specifies its location within a block. Similarly, each block also has three spatial coordinates (respectively blockIdx.x, block-Idx.y, blockIdx.z ) that localize it in the global grid. Figure 1 illustrates this organization for the 2D case. So each thread works with the same kernel but uses different data according to its spatial location within the grid 8 . Moreover, each block has a limited thread capacity according to the hardware in use.

Models and Application of GPU Delegation

In this section, the four selected models and their implementation are introduced. Then, the application of the GPU delegation on each of these models is described.

Game of Life

The first model, which was also in P&A study, is the Conway's Game of Life [START_REF] Gardner | Mathematical games: The fantastic combinations of John Conwayś new solitaire game life[END_REF]. This famous model shows that complex patterns can emerge from the implementation of simple rules. While this model does not contain agents (it is a cellular automaton), it defines an environmental dynamics which is representative of those encountered in MAS.

The environment of the Game of Life is a cellular automaton made of a two dimensional grid of square cells, each of which being in one of two possible states: Dead or alive. Each cell interacts with its eight neighbors according to the following rules: (1) Any living cell with less than two or more than three alive neighbors dies, (2) any living cell with two or three living neighbors stays alive and (3) any dead cell with exactly three living neighbors becomes alive.

In our experiments, the grid is initialized randomly 10 . At each time step, all the cells are updated according to the previous rules. Figure 2 illustrates the simulation steps.

Applying GPU Delegation.

The main computational part of the model is located in step (1) which consists in computing a sequential loop: It calculates the new state of each cell for the next step of the 8 Thread is similar to the concept of task: A thread may be considered as an instance of the kernel which is performed on a restricted portion of the data depending on its location in the global grid (its identifier). 9 e.g. http://www.jcuda.org 10 The probability for the cell to be dead or alive is the same. To create the GPU module correctly, it is necessary to focus on the representation of data. Especially, to avoid expensive transfers between CPU and GPU, the data need to be sent only once at each step. To this end, each cell write its state value (1 for alive and 0 for dead) in a 2D array (matching the size of the environment) according to its position. Then, this array is sent to the GPU. This module does the sum of the states (more precisely the sum of the number of alive cells) of all Moore neighborhood cells for each cell of the grid and stores the result in the 2D result array. Thus, each cell of the result array contains a value between 0 (no cell alive) and 8 (all neighbors are alive). The result array is then used to update the grid and compute the next simulation step according to the transition rules. Algorithm 1 presents an implementation of the corresponding GPU kernel.

Schelling's Segregation

The second model is a variant of the Schelling's Segregation model [START_REF] Schelling | Micromotives and Macrobehavior[END_REF]. This project models the behavior of two types of agents in a neighborhood: Red agents and green agents. These agents are scattered across a two dimensional grid. Their purpose is to be happy by staying near likecolored agents (each red agent wants to live near at least some red agents, and the same holds for the green agents). If they are dissatisfied at their position, the agents attempt to move to a new random vacant location for finding a better place according to their objective of happiness.

The simulation shows how these individual preferences ripple through the neighborhood, leading to large-scale patterns. In our experiments, green and red agents are scattered with the same proportion, and distributed randomly over the environment. Figure 3 illustrates the simulation steps.

Applying GPU Delegation.

The most intensive computations are in step (2) and (3) which consist, for each agent, in recovering a neighbor list and then computing its happiness according to the list. So, many sequential loops have to be done at each time step, and the required computation time clearly increases depending on the number of agents. According to the GPU delegation, the happiness computation can be deported into an environmental dynamics because these computations do not modify agents' states.

In the CPU model, the computation of happiness is done by testing the color of the agents present in the list of neighbors and counting the number of agents in each community according to their state: 1 for green agents and -1 for red agents. However, if we want to deport this computation, the data structure sent to the GPU needs to be adapted. So, at each time step, agents write in a 2D array (matching the size of the environment) their states depending on their position. This table is sent to the GPU that computes the sum of the neighboring states and returns in a result array the values for each cell of the environment. This result array thus contains values between -8 (all agents around are red) and 8 (all agents around are green). The agents then recover the value in the result array with respect to their position and act accordingly. Algorithm 2 presents an implementation of the corresponding GPU kernel.

Fire

Inspired by a model from the Netlogo library, this third model is called Fire and simulates the spread of a fire through a forest. It shows that the fire's chance of affecting the most possible number of trees in the forest depends critically on the density of trees.

In our model, all the trees are agents placed randomly in the environment (a two dimensional grid). Trees can be alive, burned or dead. When it burns, a tree releases heat which spreads in the environment. This heat can ignite other trees around: A tree ignites when the temperature is above a defined threshold. The threshold of each agent is randomly set within a range of values. A tree dies when its life reaches zero: The life of the tree decreases when it burns. Figure 4 illustrates the simulation steps. The most greedy computation loop of this model is in step [START_REF]4th International Workshop, E4MAS 2014 -10 Years Later[END_REF]. Indeed, the environment must compute the heat diffusion which requires making a global sequential loop over all the cells. Thus, the computation time increases very quickly according to the size of the environment. Because the heat diffusion is already an environmental dynamics, the delegation of this computation is very easy because it can be directly translated into a GPU module.

So, in the GPU model, the heat diffusion is done as follows: At each time step, agents add in a 2D array (matching the size of the environment) the heat that they release according to their state (alive, burn or dead). Then, this array is sent to the GPU module that compute the sum of heat values from neighboring cells (Moore neighborhood here) for each cell of the environment. More preciselly, the sum consists in adding the heat values already present in the environment (from the previous steps) and the heat generated by the agents, all modulated by a diffusion variable. Once the computation performed, the agents recover the heat value in the array with respect to their position and act accordingly. Algorithm 3 presents an implementation of the corresponding GPU kernel.

DLA

The last model is also inspired from a Netlogo model and is called DLA. DLA demonstrates diffusion-limited aggregation, in which randomly moving particles stick together to form beautiful treelike branching fractal structures. There are many patterns found in nature that resemble the patterns produced by this model: Crystals, coral, fungi, lightning, and so on.

In this model, particles are agents (initialy red), that move in a random way over the all environment (a two dimensional grid). Randomly, one of the agents stops and stays at the same position and in turn changes its color to green. Then, when a moving red agent encounters a motionless green agent, it stops and changes its color to green while the other agents continue to move in a random way. Figure 5 illustrates the simulation steps. The most greedy computations are in step ( 2) and ( 3) which consist for each agent in recovering a neighbor list and then searching within this list the nearest neighbors. So, the computation time of this model greatly increases according to the number of agents. According to GPU delegation, the agent action which consists in checking if one of the neighbor agents is green (motionless) can be converted into an environmental dynamics and then performed by a GPU kernel because this computation do not modify agents' states.

However, it is necessary to adapt the data structure to enable this transformation. Thus, all agents report, according to their position, their presence in a 2D array (matching the size of the environment): 1 if an agent occupies the cell, 0 otherwise. This array is sent to the GPU that calculates for each cell if there are agents around it. All the cells of the result array thus contain a value representing the number of neighbor agents (0 for an empty cell, 1 or more means there are neighbors around). So, the agents only have to recover the result value in the array and adjust their behavior accordingly. Algorithm 4 presents an implementation of the corresponding GPU kernel.

EXPERIMENTATIONS AND ANALYSIS

In this section, the main idea is not only to evaluate the performance of the experiments but to reflect on the accessibility, modularity and ease of use of the GPU delegation principle (the conceptual part of these implementations).

Experimentations

The GPU delegation principle produces implementations which are very different from those made by P&A. Moreover, the hardware and associated tools have also greatly evolved since P&A's work. So, we cannot directly compare the performance results of P&A with the experiences that we conduct here.

To test the performance of our implementations, each model is simulated in its CPU and hybrid (CPU + GPU) versions. Moreover, each model is simulated for different environment sizes and various densities of agents. Each simulation is executed several times over a period of 10,000 time steps. We average then the execution time for an iteration (lower execution time is better) allowing to compare the execution performance of each implementation.

For those tests, the configuration is composed of an Intel i7-4770 processor (Haswell generation, 3.40 GHz) and an Nvidia K4000 graphics card (Kepler architecture, 768 CUDA cores). Here is the list of experimentations (with setup) conducted and their results:

• Game Of Life (Figure 6): Environment size 256, 512, 1 024, 2 048; Fixed density of agents: 50%.

• Segregation (Figure 7): Environment size 256, 512, 1 024, 2 048; Fixed density of agents: 90%.

• Fire (Figure 8): Environment size 256, 512, 1 024, 2 048; Density of agents ranging from 10% to 100%.

• DLA (Figure 9): Environment size 256, 512, 1 024, 2 048; Density of agents ranging from 10% to 90%.

Analysis

In this section, we analyze the results from both a performance and a conceptual perspective by taking into account accessibility, genericness and ease of use of GPU delegation.

Performances

From Figure 10, we notice that performance gains vary significantly depending on the simulated model and the size of the environment. With GPU delegation, the gain can reach x14 but is more likely between x2 and x5 which is still a good result. Concerning the DLA model, the performance decreases when increasing the size of the environment (which is the opposite of what can be seen with the other models). This decrease is explained by the consumption of resource by the agents that slows the simulation.

P&A's implementations are much more efficient and allow to have greater performance gains that can reach x40 (particularly for the Game of Life model) but generally the gain is around x15 [START_REF] Perumalla | Data parallel execution challenges and runtime performance of agent simulations on GPUs[END_REF]. The differences are huge but must be put in perspective because we have no informations about how the models and CPU implementations were created. Nonetheless, the order of magnitude of performance gains cannot be ignored. This difference is explained by the implementation technique used by P&A. Indeed, they used an all-in-GPU approach that consists in running entirely the model on the GPU with the objective of obtaining maximum performances. GPU delegation offers less impressive performances but brings other significant benefits.

Genericness, Modularity and Reusability

GPU delegation brings more genericness thanks to its hybrid conception. Especially, the created GPU modules can be reused in other simulations. That is the case here: The module created for the Game of Life model was reused in the DLA model. Only the data sent have been adapted. This is a major difference compared to the work and implementations of P&A, and more generally with works using an allin-GPU approach. Indeed, with such works, the model used is completely redesigned and reprogrammed to run entirely on the GPU and thus cannot be easily reused. Therefore, the solution created is only punctual, which is not satisfying from a software engineering perspective.

On the contrary, more genericness could be obtained with an hybrid approach. This is especially the case with GPU delegation because it allows to reuse the created GPU modules, thus saving a lot of development time.

From an accessibility perspective, GPU delegation increases GPGPU accessibility because it promotes a modular design of the model that produces very simple kernels requiring only very little GPGPU knowledge. Moreover, by deporting only a specific part of the agent computations, it is thus possible to take advantage of the computing power of GPGPU without changing the agent model.

Finally, this modular design allows to choose what is run on the CPU and executed by the GPU. The strategy put forward by GPU delegation is to identify in the model the computations that can be transformed into environmental dynamics and thus translated into GPU modules. Considering the difficulties faced when implementing the behavior of agent on GPU or providing generic tools to do it, this approach greatly eases the use of GPGPU in MABS context.

Toward a Dedicated Methodology

The application of GPU delegation as generic approach for using GPU programming in MABS has shown many advantages whether in [START_REF] Michel | Translating Agent Perception Computations into Environmental Processes in Multi-Agent-Based Simulations: A means for Integrating Graphics Processing Unit Programming within Usual Agent-Based Simulation Platforms[END_REF] or [START_REF] Hermellin | GPU Environmental Delegation of Agent Perceptions: Application to Reynolds's Boids[END_REF]. The experiment conducted in this article confirms it. However, in view of these studies, there are still some open questions such as the difficulty of implementation for a new GPGPU user and the types of models on which GPU delegation could be applied.

Indeed, for now, we have only tested GPU delegation on models that we know compatible. This type of model has the following features: The agents are reactive and operate in a discretized environment. These agents communicate and exchange informations through the environment. Next step will be to consider models with continuous environments or agents having cognitive abilities.

About the implementation difficulties, GPU delegation aims at easing the integration of GPGPU in MABS models by promoting a real delegation strategy that ends with very simple kernels wich are easy to implement and require only very little GPGPU knowledge. However, it will be necessary to improve the support around this approach.

But more generally, our objective is to offer a generic approach which can be applied on a wide variety of models and that eases the use of GPGPU. An approach that is also able to define if the model can benefit from GPGPU power before starting the adaptation process and could determine whether there is an interest in applying the GPGPU delegation on the selected model.

From the previous experiment, a recurring process has appeared during the application of GPU delegation. We have indeed followed an iterative process that allows us to apply our approach on each model. Once formalized, we believe that this process can become a methodology dedicated to GPGPU for MABS. A methodology that would (1) help potential users to decide if they could benefit from GPGPU considering their models and (2) support the translation of MABS model into GPU programming without hiding the underlying technology. And finally such a methodology would help to spread the GPGPU in the MABS community.

SUMMARY AND FUTUR WORKS

In this paper, we proposed to study through experiments if the conclusions and issues outlined by P&A in 2008 are still true despite the evolution of GPGPU and MABS. This paper was motivated by the fact that there are actually only very few publications dealing with the use of GPGPU in a MABS context, despite the computing power that these simulations require.

To address some of the challenges that remains when using GPGPU for MABS, we apply GPU delegation on four models and show that it helps to solve some of the problems outlined by P&A. Based on a hybrid approach, GPU delegation proposes to identify in the model what are the computations that can be transformed into environmental dynamics and then translated into GPU modules.

Instead of focusing on accessibility by hiding the use of GPGPU, GPU delegation aims at easing the direct integration of GPU programming in MABS models. So, GPU delegation helps he user during the adaptation process contrary to some solutions which seek to make transparent the use of GPGPU.

This principle was applied on four models and an experimentation which compares the GPU-based runtime speed with the speed of equivalent CPU-based models was conducted. The results of experiments and analysis were then proposed.

From a performance perspective, GPU delegation offers less significant performance gains compared to all-in-GPU solutions. As highlighted by P&A, this performance gap can only be obtained at the expense of the accessibility, modularity and genericity and reflects of a totally different use of the GPU architecture.

From a conceptual point of view, we have seen that GPU delegation improves genericity (reusability of the created GPU modules), provides a more modular design (with the hybrid approach) and enhances accessibility thanks to the fact that we are no longer using an all-in-GPU approach. More precisely, its reusability and design ease the direct use of GPGPU in a MABS context. Given that GPU does not run the entire simulation but only some parts, the adaptation of the model is more simple and the knowledge required is less important. Indeed, the reusability allows to save the time invested and the efforts done to use this technology because it concretely relies on an AOSE perspective thate promotes a clear separation of concerns, and thus modularity and reusability.

Thanks to the evolution of both the GPGPU field and associated tools, GPU delegation shows that it can be easy to implement and efficient while providing significant conceptual advantages. The next step is to develop the user support so that GPU delegation can be simple to use. To this end, it is now necessary to formalize its use and all conception steps.

Especially, during implementation of the four models, we identified some repetitions in the application of the principle. Moreover, some patterns and similarities in the source code can be observed. So, an iterative process begin to appear: A list of the steps necessary to transform the model in order to use GPGPU through GPU delegation.

Based on this work, our next objective is to create a comprehensive design methodology for GPGPU in a MABS context based on the GPU delegation principle. This methodology will allow to take a model and then will help in every step necessary in the translation process in order to let the model takes advantage of GPGPU power. Such a methodology should allow to use GPU programming on a wide variety of models and spread this technology in the MABS community.

Figure 1 :

 1 Figure 1: Thread, blocks, grid organization

Figure 2 :

 2 Figure 2: Game of Life simulation

Algorithm 1 :

 1 Alive cells Kernel input : width, height, statesArray output: resultArray (the number of alive neighbors) 1 i = blockIdx.x * blockDim.x + threadIdx.x ; 2 j = blockIdx.y * blockDim.y + threadIdx.y ; 3 sumOf State = 0 ; 4 if i < width and j < height then 5 sumOf State = getNeighborsState(statesArray[i, j]); 6 end 7 resultArray[i, j] = sumOf State ;

Figure 3 :

 3 Figure 3: Segregation simulation

Figure 4 :

 4 Figure 4: Fire simulation

Algorithm 3 :

 3 Heat diffusion Kernel input : width, height, heatArray, radius output: resultArray (the quantity of heat) 1 i = blockIdx.x * blockDim.x + threadIdx.x ; 2 j = blockIdx.y * blockDim.y + threadIdx.y ; 3 sumOf Heat = 0 ; 4 if i < width and j < height then 5 sumOf Heat = getNeighborsHeat(heatArray[i, j], radius); 6 end 7 resultArray[i, j] = sumOf Heat * heatAdjustment ;

Figure 5 :

 5 Figure 5: DLA simulation

Algorithm 4 :

 4 Presence detection Kernel input : width, height, presenceArray output: resultArray (the number of neighbors around) 1 i = blockIdx.x * blockDim.x + threadIdx.x ; 2 j = blockIdx.y * blockDim.y + threadIdx.y ; 3 sumOf Agents = 0 ; 4 if i < width and j < height then 5 sumOf Agents = getNeighborsPresence(presenceArray[i, j]); 6 end 7 resultArray[i, j] = sumOf Agents ;

Figure 6 :Figure 7 :

 67 Figure 6: Results for the Game of Life model

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Results for the Fire model

e.g. http://repast.sourceforge.net/repast hpc.php

e.g. https://developer.nvidia.com/about-cuda

e.g. https://developer.nvidia.com/what-cuda

e.g. http://www.khronos.org/opencl

http://www.lirmm.fr/~hermellin/Website/Reynolds_ Boids_With_TurtleKit.html

e.g. http://www.turtlekit.org

e.g. http://docs.nvidia.com/cuda/