Abdallah Arioua

Patrice Buche

Madalina Croitoru

A Datalog+/-Domain-Specific Durum Wheat Knowledge Base

We consider the application setting where a domain-specific knowledge base about Durum Wheat has been constructed by knowledge engineers who are not experts in the domain. This knowledge base is prone to inconsistencies and incompleteness. The goal of this work is to show how the state of the art knowledge representation formalism called Datalog± can be used to cope with such problems by (1) providing inconsistency-tolerant techniques to cope with inconsistency, and (2) providing an expressive logical language that allows representing incomplete knowledge.

Introduction

The Dur-Dur research project 1 aims at restructuring the Durum Wheat agrifood chain in France by reducing pesticide and fertilizer usage while providing a protein-rich Durum Wheat. The project relies on constructing a multidisciplinary knowledge base (involving all actors in the agrofood chain) which will be used as a reference for decision making. This knowledge base is collectively built by several knowledge engineers from different sites of the project. Due to various causes (errors in the factual information due to typos, erroneous databases / Excel files, incomplete facts, unspoken obvious information "everybody knows" etc.) the collectively built knowledge base (KB) is prone to incompleteness and inconsistencies. Incompleteness has many forms, in our case it reflects itself as a lack of precision and explicitness. For instance, an expert may say that the Durum Wheat is contaminated by a mycotoxin but he/she may, for some reasons, do not specify which mycotoxin. Inconsistency appears as logical contradictions due to the causes stated above. The problem is that in presence of inconsistencies the knowledge base becomes unreliable and not trustworthy, let alone the fact that reasoning under inconsistency is challenging for logical formalisms.

To solve the above mentioned problems, we propose in this paper a methodology of representing Durum Wheat knowledge in the logical framework of Datalog± [START_REF] Calì | A general datalog-based framework for tractable query answering over ontologies[END_REF][START_REF] Mugnier | An introduction to ontology-based query answering with existential rules[END_REF]. Datalog± is expressive enough to allow the representation of unknown individual in the knowledge base and cope with heterogeneous data as it allows for n-ary predicates. Moreover, Datalog± has an interesting equivalent relation with conceptual graphs [START_REF] Sowa | Conceptual graphs for a data base interface[END_REF], in fact any logical formula in Datalog± can be translated to a graphical representation, which significantly helps experts in other domains in the process of knowledge acquisition. We present with detailed examples how this methodology is used to construct the Durum Wheat knowledge base in the French project Dur-Dur. The knowledge base is available online at http://www.lirmm.fr/ ˜arioua/dkb/ where the reader can find downloadable materials.

The Logical Language Datalog±

There are two major approaches in the knowledge representation community: Description Logics (DL) (such as EL [START_REF] Baader | Pushing the el envelope[END_REF] and DL-Lite [START_REF] Calvanese | Tractable reasoning and efficient query answering in description logics: The dl-lite family[END_REF] families) and rule-based languages (such as Datalog± language [START_REF] Mugnier | An introduction to ontology-based query answering with existential rules[END_REF][START_REF] Calì | A general datalog-based framework for tractable query answering over ontologies[END_REF], a generalization of Datalog that allows for existentially quantified variables in rule's head). Despite its undecidability when answering conjunctive queries, different decidable fragments of Datalog± are studied in the literature [START_REF] Baget | Walking the complexity lines for generalized guarded existential rules[END_REF]. These fragments generalize the above mentioned DL families and overcome their limitations by allowing any predicate arity as well as cyclic structures.

The Datalog± corresponds to the positive existential conjunctive fragment of firstorder logic, which is composed of formulas built with the connectors (∧, →) and the quantifiers (∃, ∀), with constants but no function symbol.

An atom is of the form p(t 1 , . . . , t k ) where p is a predicate of arity k and the t i are terms, i.e., variables or constants (we use vectors , e.g. -→ x , to denote a sequence of variables). A finite set of atoms F is called an atomset (a fact), we denote by terms(F ) (resp. vars(F )) the set of terms (resp. variables) that occur in

F . A homomorphism π from two atomsets A 1 to A 2 is a substitution of vars(A 1 ) by terms(A 2 ) such that π(A 1 ) ⊆ A 2 . An existential rule (or a rule) is of the form R = ∀ - → x ∀ - → y (B → ∃ - → z H),
where B and H are conjunctions of atoms, with vars(B) = -→ x ∪ -→ y , and vars(H) = -→ x ∪ -→ z . B and H are respectively called the body and the head of R. Chase is the mechanism by which one deduce new facts by rule application on the initial set of facts F. We denote by Cl R (F) the set of all facts that can be deduced from F by a set of rules R. A knowledge base K = (F, R, N ) is composed of a finite set of facts F, rules R and negative constraints N (i.e. a rule whose head is set to ⊥). A Boolean conjunctive query (BCQ or query in the following) has the form of a fact. We say a query Q is entailed from

K iff Cl R (F) |= Q. We say K is inconsistent iff Cl R (F) |= ⊥.

The Durum Wheat Knowledge Base

The Durum Wheat knowledge base has been constructed within the French National Project DUR-DUR. The goal of this knowledge base is to integrate scientific knowledge acquired from different tasks during the project to redesign the durum wheat chain. The Dur-Dur project suggests developing a systematic approach to investigate issues related to the management of the nitrogen, energy and contaminants, to guarantee a global quality of products throughout the production and the processing chain. Started in 2014 and planned over 4 years, the project aims at integrating the 3 dimensions of the sustainability (environmental, economic, and social), at 4 levels of investigation (4 tasks) with a complementary task (task 5). Figure 1 depicts the different tasks of the project where the fifth task's central role is to integrate knowledge from different tasks. The Durum Wheat knowledge base is the product of the fifth Task. It will be used in many computational tasks, notably analyzing and comparing the alternative innovative technical itineraries proposed in the project to reduce the use of chemical inputs (nitrogen fertilizers and pesticides). The knowledge base represents domain-specific knowledge about Agronomy. It is composed of four main parts:

-Vocabulary: it contains knowledge about concepts and relations.

-Rules: they represent rules that encode generic knowledge.

-Negative constraints: this part contains constraints about crops and Agronomyrelated constraints. -Facts: this part contains factual knowledge about Agronomy-related subjects (fertilizers, pesticides, diseases,etc.). In the next section we start by highlighting the guidelines which were followed to author the knowledge base (Subsection 3.1) then we turn to the the internal structure or the architecture of the knowledge base including the vocabulary, the rule-base alongside with the constraints and the factual knowledge (Subsection 3.2).

The Authoring

A multidisciplinary process of knowledge acquisition and representation was deployed to author the knowledge base. We used technical reports to define the scope of the knowledge base and the relevant concepts of our vocabulary. Taking into account the recommendation of [START_REF] Thunkijjanukij | Lesson learned for ontology construction with thai rice case study[END_REF], we followed three steps specification, conceptualization and formalization to build the knowledge base. pesticides and optimize energy consumption using a systematic approach that makes use of innovative technical itineraries. The contribution of the knowledge base lays in offering an expressive way of representing domain-knowledge.

Conceptualization. The concepts and the relations among them alongside to rules, facts and constraints have been defined and collected from technical reports (see Figure 3) and online materials (see [START_REF] Arvalis | Les fiches arvalis: Varits produits accidents en grandes cultures[END_REF]). It is worth mentioning that in the vocabulary part we have built on the vocabulary of Agropedia indica ( [START_REF] Sini | Building Knowledge Models for Agropedia Indica v 1.0 Requirements, Guidelines, Suggestions[END_REF]) with an increase (and modification) in content that approximates 60%. 2Formalization. Since understanding logical formulas is quite difficult for experts who are not familiar with KRR formalism we have chosen a graphical framework (conceptual graphs; [START_REF] Sowa | Conceptual graphs for a data base interface[END_REF][START_REF] Sowa | Conceptual structures: information processing in mind and machine[END_REF]) to author the knowledge base. Moreover, the conceptual graphs (CGs) made it easy for the Agronomy experts to understand the content of the knowledge base. Furthermore, CGs enjoy the same expressive power as Datalog±. In fact, it is an equivalent formalism of Datalog± as shown in [START_REF] Chein | Graph-based Knowledge Representation -Computational Foundations of Conceptual Graphs[END_REF]. Therefore, our choice was to choose CGs for knowledge acquisition and Datalog± as a framework for reasoning. For CGs, we used CoGui 1.6b which is an IDE for representing and reasoning with CGs. 3 We shall explain in-depth in Section 3.2 the graphical and logical representation for each part of the knowledge base. The facts within the knowledge base are exported to an RDF/XML format whereas the vocabulary, rules and constraints are exported as DLGP format (DataLoG Plus; [START_REF] Leclère | Kiabora: an analyzer of existential rule bases[END_REF]). The vocabulary of the knowledge base contains 279 concepts and 116 relations, the rule-base contains 23 rules and the constraints part contains 25 constraints. The factual part has around 900 atoms. The knowledge base is available online at http://www.lirmm.fr/ ˜arioua/dkb/. Fig. 3: Some snapshots of the technical reports.

The Structure

As depicted in Figure 2 the knowledge base is composed of four parts. It is worth mentioning that on the logical level the vocabulary and the rule-base are the same. However, we adapt here the the Semantic web notation and we differentiate between them. Therefore, we distinguish between those rules that express logical consequences (in the rule-base) and those that encode generalizations and classes inclusions (in the vocabulary).

The vocabulary The vocabulary represents an explicit specification of the terms and concepts used in Agronomy. The vocabulary is composed of two parts: (1) concept types hierarchy and (2) relation types hierarchy.

1. Concept types hierarchy: concepts are organized within a hierarchy as superconcepts and sub-concepts. For instance, the concept disease and its sub-concepts (e.g. viral disease, fungal disease, etc.), types of pesticides (e.g. herbicide, insecticide, fungicide) are all of organized in a hierarchy. 2. Relation types hierarchy: in CGs the concepts are related by relationships. Since concepts are divided into super-concepts and sub-concepts, relationships are divided in the same way. In the relation types hierarchy we find super-relations and sub-relations. For instance, the relation "useSowingProcess" which relates the seeding and sowing production step with the process of sowing (which is a superconcept of broadcasting, behind plough and a sub-concept of process). This relation is a sub-relation of the super-relation "useProcess" that relates any production step with any process.

In CGs the hierarchy of concept types is represented as in the upper graph of Figure 4. Rectangles represent concepts and the arrow represents the generalization between them where the source of the arrow is the sub-concept and the target of the arrow is super-concept. In the relation types hierarchy (the lower graph), the circles are the relations and the arrows are generalizations.

To better illustrate the relation between existential rules and CGs let us take an example that shows the transformation of some part of the graphs of Figure 4 to their logical form.

Example 1. The left-most part of the concept types hierarchy that indicates that "Viral disease is a disease" is represented logically by a rule as follows:

-∀x(V iral disease(x) → Disease(x)).
The part of the relation types hierarchy that indicates that "Using Herbicide is using Pesticide" is represented logically by a rule as follows:

-∀x, y(useHerbicide(x, y) → useP esticide(x, y)). The mycotoxin z is unknown (it could be Aflatoxins, Deoxynivalenol, etc.) but still the information that "there is necessarily a mycotoxin" is present, which is an important information when it comes to risk management where a possible contamination by any mycotoxin is taken to be critical. Moreover, the importance of such representation manifests also in helping knowledge elicitation where the knowledge base can make use of incomplete information and then be updated incrementally by identifying the existential variables.

In conceptual graphs the rule (b) is depicted in Figure 5. In a rule, the rectangles are called concept nodes and the circles are called relation nodes. A concept node has a concept type and a marker which can be either an individual marker (constant) or a generic marker (a variable denoted as *). For instance, the concept richOf has a generic marker (*) which represents a variable. If the marker were an individual marker we should have found a constant name like Nitrogen. The relation nodes are predicates that relate different concepts. A rule in conceptual graphs is composed of two parts, a hypothesis (left) and a conclusion (right). The dashed lines link those concepts that share the same variables (called frontier variables). That means, variables that appear in the hypothesis and in the conclusion. In the rule (b), the concept weed in the hypothesis part shares the same variable with the concept weed in the conclusion. As we said earlier, certain classes of Datalog± render the inference undecidable. However, there are some classes that ensure decidability. Most notably, FUS (Finite Unification Set) and FES (Finite Expansion Set) classes. The online tool Kiabora deploys syntactical and semantic analysis on any set of rules written in the DLGP format. 4The tool, for a given rule-base, classifies all the rules with respect to the known classes. From the analysis we found that our rule-base lays within the decidable classes. Specifically, FUS and FES.

The negative constraints

Representing what cannot be allowed within certain domain of interest is called negative constraint (or constraint). Consider the following negative constraint:

(c) ∀x, y, z(Soil(x) ∧ M aize(y) ∧ Durum wheat(z) ∧ hasP recedent(x, y) ∧ isCultivatedOn(z, x) → ⊥).
This negative constraint forbids using Maize as a precedent on a soil if we want to cultivate Durum Wheat on this soil. Figure 6 represents the CGs representation of this negative constraint. Besides this type of constraints we have the banned types constraints. These are particular forms of constraints that express concept disjointness. For Fig. 6: The negative constraint (c) in the CGs framework. instance, a soil x cannot be a disease, ∀x(Soil(x) ∧ disease(x) → ⊥). In the Durum Wheat knowledge base all concepts are disjoint except those concepts which have a generalization/specialization relations among them.

The factual knowledge In the Durum Wheat knowledge base the factual part represents domain-specific knowledge. This knowledge is divided into two parts: (1) general factual knowledge and (2) knowledge about different technical itineraries. According to [START_REF] Sebillotte | Itineraires techniques et evolution de la pensee agronomique[END_REF] a technical itinerary is a "logical organized course of technical actions applied to a cropped species".

General factual knowledge is the part of the knowledge base that represents general facts about the domain, for instance, Miradoux is a variety of Durum Wheat or the fungal disease Fusarium Flag smut is cause by, among other causes, the fungi Urocyctis agropyri of the family Fusarium. The following is an example of a set of facts. Recall that commas are interpreted as conjunctions.

(d) {F ungal disease(F lag smut), isCausedBy(F usarium ear blight, U rocyctis agropyri), f ungi(U rocyctis agropyri)}.

Here we have the relation isCausedBy instantiated on the individuals F lag smut and U rocyctis agropyri. The former is a fungal disease as stated by the concept F ungal disease and the latter is a fungi. Figure7 depicts the set of facts in the CGs framework. The second part of the factual knowledge part are those facts about the technical itineraries. In what follows we give a real-world example of a well-known technical itinerary in France. Example 2. This example represents the reference technical itinerary in France which is followed by farmers to cultivate their fields.

"The variety to be seeded in the soil is Miradoux, the culture precedent is sunflower. The soil is prepared means of harrowing. The seeding is done with density of 280 grains/m 2 . Fertilization is to be performed at the growing stage when the tiller begins with dose 40u and 50u at the end of the tiller."

This technical itinerary is a set of facts, e.g. "variety is Miradoux", "Fertilization is to be performed at the growing stage", etc. However, not any any set of facts. Particularly, it is a precise set of describing facts. Actually, any ITK (according to the studied reports) should precisely account for the following steps:

1. Variety to bee seeded. 2. Date of seeding alongside the density. 3. Cultural precedent. 4. Inter-cropping techniques. 5. Soil preparation method. 6. Disease management method. 7. Weed management method. 8. Insect control method.

Thus a technical itinerary should be mainly composed of these describing facts. The fowling is a snippet of the technical itinerary described in Example 2.

FIT K =                        Soil(Soil1)
Durum wheat(D1) isOf V ariety(D1, M iradoux)

V ariety(M iradoux) isCultivatedOn(D1, Soil1)

Seeding and sowing(Seeding1) Seed(Seed1) useSeed(Seeding1, Seed1) seedOf (Seed1, Durumw1) isAppliedOn(Seeding1, Soil1) withDensity(Seeding1, Density1) Density(Density1) U nit(grain mm) hasV alue(280) V alue(280)

What has been presented so far is the Durum Wheat knowledge base we have constructed manually within the project, which can be seen as a contribution on itself. since what is mainly proposed by researchers in this field are ontologies. Besides that, our knowledge base provides querying facilities not only on the ontological layer but also on the factual layer where real knowledge about the domain is represented in form of facts.

In the Durum Wheat knowledge base each technical itinerary is stored separately from the other technical itineraries (we have three in total). On can query them all together or separately.

Reasoning

The first and foremost reason to acquire knowledge and store it in knowledge bases is to provide querying facilities for the end-user. Like in classical database systems, in Datalog± the main reasoning task is query answering. The main and important difference is that in our case the querying is enriched by a rule-base layer. Thus the reasoner takes into account the domain-knowledge represented within the rules while querying.

Formally, a conjunctive query has the form of a fact but with possibly free variables. For instance Q(x) = F ungal disease(x) ∧ isCausedBy(x, culmorum) is a conjunctive query that looks for "the fungal disease that is caused by culmorum".

In order to perform reasoning in forward chaining in presence of rules, the reasoner applies all the rules in the rule-base on the set of facts in the factual part then query the knowledge in a classical manner.Given a set of facts F and a set of rules R This means that the chase computes all deducible knowledge of F by the application of all the rules of R on all the facts on F until no rule will be applicable. This process is also called saturation. Note that if the closure of a set of facts F is the same as F, i.e. Cl R (F) = F, then we say that F is closed under the application of rules (or deductively closed). A query Q has an answer within a knowledge base K iff Cl R (F) |= Q where |= refers to the usual first-order entailment. It may happens that the set facts F contains contradictory knowledge (i.e. inconsistencies). We say that a set of facts is inconsistent iff Cl R (F) triggers a negative constraint. The solution [START_REF] Lembo | Inconsistency-tolerant semantics for description logics[END_REF] is to construct maximal (with respect to set inclusion) consistent subsets of F. Such subsets are called repairs and denoted by Repair(K). They actually represent possible distribution of facts to restore consistency. Once the repairs are computed, different semantics can be used for query answering over the knowledge base. After repairing the knowledge base we can query it using different semantics. The most common semantics is to query the intersection of all repairs. This is a cautious strategy because the intersection is practically those facts which are not involved in any inconsistency.

Conclusion

In this paper we have presented a general methodology to build Durum Wheat knowledge bases within the logical language Datalog±. We presented detailed examples and a real-world Durum Wheat knowledge base which has been built within the French national project Dur-Dur. The expressiveness of Datalog± lays in its ability to deal with incompleteness and inconsistency. Moreover, it has an interesting relation with Conceptual Graphs which makes it easy to non-experts to manipulate and understand logical formulae. In addition, DLGP format (DataLoG Plus; [START_REF] Leclère | Kiabora: an analyzer of existential rule bases[END_REF]) can be translated to semantic web languages as OWL/RDFS using COGui or GRAAL framework [START_REF] Baget | Graal: A toolkit for query answering with existential rules[END_REF].

Fig. 1 :

 1 Fig. 1: The different tasks of the Dur-Dur project. The knowledge tasks aims at integrating multidisciplinary knowledge from other tasks.

Specification.

  Fig. 2: An overview of the Durum Wheat knowledge base. The circles contains knowledge examples represented in the conceptual graph framework.

Fig. 4 :

 4 Fig. 4: Concept and relation types hierarchy.

Fig. 5 :

 5 Fig. 5: The rule (b) in the CGs framework.

Fig. 7 :

 7 Fig. 7: The set of facts (d) about fungal diseases and fungus.

Example 3 .

 3 Consider the following knowledge base K:F = {D(a), S(b)}, R = {∀x(D(x) → C(x)), ∀x, y(S(x) ∧ C(y) → M (x, y))}, N = ∅. The closure is Cl R (F) = {D(a), S(b), C(a), M (a, b)}.

Example 4 .

 4 Consider: F = {D(a), S(b), P (c)}, R = {∀x(D(x) → C(x))}, N = {∀x, y(S(x) ∧ C(x)) → ⊥}. Then the negative constraint will be triggered after the application of the rule which infers C(a). Therefore our repairs would be A 1 = {D(a), P (c)} and A 2 = {S(b), P (c)} and Repair(K) = {A 1 , A 2 }. While Cl R (A 1 ) = {D(a), C(a), P (c)} and Cl R (A 2 ) = A 2 .

http://www.agropedia.net.

http://www.lirmm.fr/cogui/.

Kiabora 0.1 website: http://www.lirmm.fr/ ˜mugnier/graphik/kiabora/, see[START_REF] Leclère | Kiabora: an analyzer of existential rule bases[END_REF] for a detailed explanation.


Hugo Daigmorte
email: hugo.daigmorte@onera.fr

Marc Boyer
email: marc.boyer@onera.fr

Traversal time for weakly synchronized CAN bus

Keywords: reference → Performance, •Mathematics of computing → Numerical analysis, •Computer systems organization → Real-time systems, worst-case traversal times, network calculus, CAN bus

Scheduling frames with offsets has been shown in the literature to be very beneficial for reducing response times in realtime networks because it allows the workload to be better spread over time and thus to reduce peaks of load. Maintaining a global synchronization amongst the stations induces substantial overhead and complexity in networks not providing a global time service such as CAN. Indeed, on CAN, no global clock is implemented in practice and each station possesses its own local clock. Without a global clock, the de-synchronization between the streams of frames created by offsets remains local to each station. The first contribution of this work is to show that important gains with respect to the communication latencies, around 40% in our experiments, can be achieved if we implement bounded clock desynchronization, also refered to as bounded phases, between the stations. The second contribution of this work is to provide a set of network-calculus based timing analyses to handle systems with bounded phases and compare their performances.

INTRODUCTION

Controller Area Network (CAN), a serial communication bus network, was initially developed for automotive applications. Due to the many advantages of CAN, including its high reliability and cost effectiveness, it has found application in other industries. In particular the standard AR-INC 825, was developed to standardize the use of CAN in aerospace domain [START_REF] Committee | Arinc specification 825-2: General standardization of CAN (Controller Area Network) bus protocol for airborne use[END_REF].

Figure 1: Example of CAN topology

It is well known that the use of offsets reduces response time [START_REF] Tindell | Adding time-offsets to schedulability analysis[END_REF] and increases the bus usage [START_REF] Grenier | Pushing the limits of CAN-scheduling frames with offsets provides a major performance boost[END_REF][START_REF] Yomsi | Controller Area Network (CAN): Response time analysis with offsets[END_REF], because offsets allow the workload to be spread over time and thus to reduce peak load and avoid contention and so to reduce the worst-case response times (WCRT) and to permit a better bandwidth utilization [START_REF] Grenier | Pushing the limits of CAN-scheduling frames with offsets provides a major performance boost[END_REF]. However using offsets requires a clock. In distributed systems there are two main solutions: each node having a local clock or all nodes sharing a global clock. In both, each message is sent at an offset with regard to a clock. In case of global clock each node has (up to some precision) the same clock value, and no contention occurs, neither between flows from the same node, nor from different nodes, for instance in TTCAN [START_REF] Leen | TTCAN: a new time-triggered controller area network[END_REF] or TTP [START_REF] Kopetz | TTP-a protocol for fault-tolerant real-time systems[END_REF]. However synchronized local clocks require synchronization mechanisms, and the clock precision must be smaller than the sending time of one frame. In case of local clocks, the scheduling remains local, synchronization mechanisms between stations are no longer a concern. Using local clocks avoids the contentions between flows from the same network, and reduces the contentions between flows from different nodes, as will be illustrated in Section 3.

Inspired by the "bandwidth management" recommendation of [START_REF] Committee | Arinc specification 825-2: General standardization of CAN (Controller Area Network) bus protocol for airborne use[END_REF], we wondered what could be the gain of some intermediate solution: a system with a global clock but a weak precision, that can also be seen as a system with local clocks, where the phases between the clocks are bounded (see Figure 1).

Nevertheless, there was no method, up to our knowledge, to compute response time of such a system. Then, we have adapted network calculus to model an exact synchronization between flows from the same node, and bounded synchronization between flows from different nodes.

Section 2 will present an overview of the related work, in the area of offsets, synchronized networks and the associated analyze methods. Section 3 will give more details on the approach: the hypotheses, the expected benefits and the kind of systems it can model. Section 4 will recall some network calculus background and present some new properties, required in this context. Finally, numerical results on an example will be given in Section 5.

RELATED WORK

The timing analysis of CAN has been the object for various studies in the past. The worst case response times were first provided in [START_REF] Tindell | Analysing Real-Time Communications: Controller Area Network (CAN)[END_REF] and then revisited in [START_REF] Davis | Controller Area Network (CAN) schedulability analysis: Refuted, revisited and revised[END_REF] but without considering offsets. They gave the exact response time for sporadic messages. In [START_REF] Klehmet | Delay bounds for CAN communication in automotive applications[END_REF][START_REF] Sofack | Non preemptive static priority with network calculus: Enhancement[END_REF] it has been shown that application of Network Calculus can bound the worst case response times, this bound is the exact response time on a large set of tests but there is no guaranty, these results can be applied to non-periodic flows.

Response times on CAN with offsets and local clocks have also been studied first with approximate but lower-complexity forms of analyse in [START_REF] Tindell | Adding time-offsets to schedulability analysis[END_REF] and then an effective worst-case response time analysis in the non-preemptive case with offsets has been given in [START_REF] Yomsi | Controller Area Network (CAN): Response time analysis with offsets[END_REF]. Multi-hop systems with local clocks and offsets have also been studied using network calculus in the case of AFDX [START_REF] Li | Improving end-to-end delay upper bounds on an AFDX network by integrating offsets in worst-case analysis[END_REF].

In case of global clock, each node has the same clock value, and no contention occurs, for instance: TTCAN [START_REF] Leen | TTCAN: a new time-triggered controller area network[END_REF], TTP [START_REF] Kopetz | TTP-a protocol for fault-tolerant real-time systems[END_REF] or TTEthernet [START_REF] Steiner | TTEthernet dataflow concept[END_REF]. These systems require the ability to build a global communication schedule and to synchronize local clocks.

Generating a communication schedule is NP-complete and involves non-trivial optimization algorithms, see for instance [START_REF] Tamas-Selicean | Synthesis of communication schedules for TTEthernet-based mixed-criticality systems[END_REF], but this schedule gives the response time of frames by avoiding contentions.

There are several solutions to synchronize clocks throughout a computer network. Some standard protocols exists for example the Precision Time Protocol (PTP), defined in [START_REF] Lee | IEEE 1588-standard for a precision clock synchronization protocol for networked measurement and control systems[END_REF]. Achieving clock synchronization can also be done with dedicated protocols, implemented in software (TTCAN level-1 [START_REF] Leen | TTCAN: a new time-triggered controller area network[END_REF][START_REF] Hartwich | Timing in the TTCAN network[END_REF]) or requiring specific hardware support (TTEthernet [START_REF] Steiner | TTEthernet dataflow concept[END_REF], TTP [START_REF] Kopetz | TTP-a protocol for fault-tolerant real-time systems[END_REF], TTCAN level-2 [START_REF] Leen | TTCAN: a new time-triggered controller area network[END_REF][START_REF] Hartwich | Timing in the TTCAN network[END_REF]...).

BOUNDED PHASES BETWEEN STATIONS

Model assumptions

In this paper we propose to study a network with local clocks coarsely synchronized.

We assume a CAN bus, with a set of nodes N and a set of periodic flow F. Each node n ∈ N has a local clock: if t is the absolute time, cn(t) is the node clock value at this instant. For each flow f ∈ F is characterized by a sender and a tuple (O f ,T f ,lmax f ) where O f is its offset, T f its period and lmax f its maximal frame size. The k-th instance of a flow f is generated each time the local clock cn(t) is equal to

kT f + O f .
Moreover, some bound is known on the difference between local clocks from different nodes, cn(t) -c n (t), called phase.

If all the phases are null, then we are in the case of global clock, if they are completely unknown, then we are in the case of local clocks.

For sake of simplicity, we assumed that there is no jitter at frame sending, and we considered only periodic flow, no sporadic. These are not limitations of the analysis method, as will be explained in the conclusion.

Expected benefits

Bounded phases are a trade-off between global clock and local clock. In this section, a small example (two nodes and three flows) is presented, illustrating the benefits of the three approaches. We consider a system of two nodes. The flows Using a global clock allows to avoid contentions and then provides small network delays. It requires building a global frame schedule such that there is no contention: such a schedule is given in Figure 2 at the top. A time slot is dedicated to each message, and no contention occurs between the flows X, Y, Z.

A,1 B,1 A,2 C,3 C,1 C,2 A,1 B,1 A,2 C,1 C,2 C,3 BUS N1 N2 BUS N1 N2 BUS BUS A,1 B,1 A,2 C,3 C,1 C,2 A,1 B,1 A,2 C,1 C,2 C,3
In case of local clocks there is no control on phases between the stations. Nevertheless offsets with respect to a local clock create some shaping and reduces contentions. These effects can be illustrated in Figure 2 at the bottom. Contentions between flows A and B from node 1 cannot happen. However contentions between flows from different node can happen: both A and B can be delayed by C. Nevertheless, local clocks create some shaping and reduces contentions between nodes: C can be delayed by at most A or B but never both of them.

Considering a network with local clocks coarsely synchronized, i.e. where the phases between nodes is not perfectly known but bounded, some contentions can be avoided. Like in the case of local clocks, no contention will occur between the flows A and B. But if the phase between N 1 and N 2 is small enough, no contention can occur between flow B and C. This shows it is possible to benefit from some of the advantages of a global clock with fewer constraints on the synchronization between nodes.

Targeted systems

Several synchronization protocols have been defined, but their implementation on a COTS can be costly, and the use of dedicated hardware is not always possible.

But even having a global reference clock can be insufficient to synchronize the bus usage: the payload of the frames is produced by some tasks, and sending a frame wrt to a clock requires that the data is produced before the frame sending time, and then gives some constraints on the task scheduling itself.

The bandwidth management requirement on ARINC 825 [START_REF] Committee | Arinc specification 825-2: General standardization of CAN (Controller Area Network) bus protocol for airborne use[END_REF] is illustrated by the use of a cyclic scheduling where a major frame (MAF) is divided as a sequence of minor frames (MIF) and possibly some idle time at end of the MAF.

Due to hardware and software latency, it may be hard to start all MAF at exactly the same time, but it could be possible to have a bounded phase between them. In the standard examples, the MAF is 50ms long, and a MIF is 10ms long.

The network calculus theory will be used to compute upper bound on the CAN bus with bounded phases. We assume that the system respects the nominal CAN behavior, and behaves as a non-preemptive static priority scheduler, where each node is always inserting its frames according to their priorities. Then, from network calculus point of view, the response time of such a system can be computed adapting existing results on non-preemptive static priority scheduler.

This section presents how it has been adapted to handle a system with offsets and bounded phases. Section 4.1 presents some reminds on network calculus. Section 4.2 gives a new expression to upper bound delays: as it will be shown in the experiments, this bound gives good results as long as the phases are small. The core of the contribution is the computation of arrival curves for aggregate flows, taking into account the offsets in a node, and the phases between nodes: this is presented in Section 4.3.

Network Calculus reminds

Mathematical background: (min;+) dioid

Network calculus is a theory to get deterministic upper bounds in networks. Network calculus mainly handles non decreasing functions, null before 0 : F. It is mathematically based on the (min,+) dioid and beyond classical operations like addition or minimum, network calculus relies on two basics operators the convolution and the deconvolution.

Notation 1. The ceiling function is denoted • , and we also introduce [x] + = max(x, 0).

Definition 1. (Min-plus convolution and deconvolution)

The min-plus convolution * and deconvolution of two functions f and g are defined by

(f * g)(t) = inf 0≤s≤t {f (t -s) + g(s)} (f g)(t) = sup 0≤u {f (t + u) -g(u)} Property 1. ([9, Th. 3.1.12]) For any f, f , g, g ∈ F : g ≥ g ⇒ f g ≤ f g (1) f ≥ f ⇒ f g ≥ f g (2) (f g) g = f (g * g ) (3) Property 2. Let f, g, h ∈ F three non decreasing func- tions, then (f + g) * h ≤ (f * h) + g. Proof. ((f + g) * h)(t) = inf 0≤s≤t {f (s) + g(s) + h(t -s)} ≤ inf 0≤s≤t {f (s) + g(t) + h(t -s)} ≤ (f * h)(t) + g(t)
Definition 2. Let d, P, h ∈ R be some parameters. Then the functions latency δ d and stair ν P,h are defined by: The servers are just relations between some input and output flows (S ∈ F × F ). Then (A, D) ∈ S, means that a server S receives an arrival/input flow, A(t), and delivers the data after some delay, it is the departure/output flow, D(t). We have relation D ≤ A, meaning that data goes out after being entered. Figure 5 shows input and output functions for a single server queue.

δ d (t) = 0 if t ≤ d ∞ otherwise. ν P,h (t) = h t P + Note that (f * δ d )(t) = f ([t -d] + ), (f δ d )(t) = f (t + d).
The main network calculus performance measures are backlog and delay (see Figure 5). Definition 3. (Backlog period) An interval I is a backlogged period iff the backlog (number of bits that are inside the system) is not null in this interval, i.e. ∀t ∈ I, D(t) -A(t) < 0 Definition 4. Delay at a time t is the delay that a bit entered at time t will wait until going out. It is the horizontal deviation (hDev) (see Figure 5).

hDev(A, D, t) = inf{d ∈ R + |A(t) ≤ D(t + d)} hDev(A, D) = sup t∈R + {hDev(A, D, t)}
However the exact input/output data flows are in general unknown at design time, or too complex, and the calculus of these bounds cannot be obtained.

Network calculus: contract modeling

The evolution of input/output data flows can be determined considering traffic contract on the traffics and the services in the network. For this purpose, network calculus provides the concepts of arrival curve and service curve.

∀(t, d) ∈ R 2 + , A(t + d) -A(t) ≤ α(d) (4)
The eq. ( 4) is equivalent to

A ≤ A * α [9, Lem. 1.2.3].
The contract on flow can be seen either from an engineering point of view (bound on any interval, see Figures 6) or from a mathematician point of view (convolution-based relation).

Property 3. Let A be a flow and α an arrival curve of A.

• If α ≥ α then α is an arrival curve of A.

• A A is the "best" arrival curve for A, i.e. A A is an arrival curve and A A ≤ α

It is important to notice that for a flow A, an infinite number of arrival curves exists and depending on the one which will be used to compute bounds on worst case traversal time the result may suffer a lack of accuracy. And so it is essential to have, if possible, the "best" arrival curve possible. Property 4. Let A be a flow and α an arrival curve of A.

If ∀d ∈ R+, ∃ t ∈ R+ such as α(d) = A(t + d) -A(t)
then α is the best arrival curve for A.

Proof. Suppose d ∈ R+ such a α(d) > A A(d) then there is t ∈ R * such as α(d) = A(t + d) -A(t) > A A(d) which is contradictory. So α = A A.
For the service definition, the two points of view exist, but lead to slightly different definitions: the strict minimal service curve and the min-plus minimal service curve. Definition 6. (Strict minimal service) A server S offers a strict service β iff for all input/output A,D and for all backlogged period (s,t] Let us now present the main network calculus results which allows, considering contracts, to compute bounds on delay.

D(t) -D(s) ≥ β(t -s)
Theorem 1. (Delay bound) Let S be a server offering a min-plus minimal service curve β. If the input flow A has an arrival curve α, then, the delay can be bounded by

hDev(A, D) ≤ hDev(α, β) hDev(A, D) ≤ (β -α) -1 (0+)
where f -1 is the pseudo-inverse of the function f . The bound hDev(α, β) is the most commonly used bound, and (β -α) -1 (0+) an upper bound on the busy period, in general greater than the other bound: cf. [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF] at Section 1.6.2, 3.1 and 6.

Theorem 2. (Non-preemptive static priority, [START_REF] Bouillard | Service curves in Network Calculus: dos and don'ts[END_REF]) Let S be a server offering a strict minimal service β, shared by three flows, A, AH , AL, AH having a higher priority than A, and AL a lower. Then, if αH is an arrival curve for AH and L max L is an upper bound on the frame size of AL, the flow A receives a residual min-plus service of curve

[β -αH -L max L ] + .
Corollary 1. Let S be a server offering an aggregate minimal strict service of curve β to a set of flows A1, . . . , An, with a non preemptive static priority scheduling policy (taking the natural order as priority order, A1 having the highest priority, An the lowest). For any Ai, if α1..i-1 (resp. α1..i) is an arrival curve for the aggregate flow i-1 k=1 A k (resp. i k=1 A k ), then two upper bounds on the delay of the flow i are:

hDev(αi, β -α1..i-1 -l max N P ) (5a) (β -α1..i -l max N P ) -1 (0+) (5b)
where l max N P represents the maximal size of a frame of low priority.

The corollary is just an application of Theorem 1 and 2. One contribution of this article is the computation of an accurate α1..i function for systems with phases.

A new bound on delay

The common bounds used in network calculus involve the arrival curves of the flows. But when the data flows are well known, when the system is deterministic, or quite deterministic, considering the real flow can lead to better results.

The Theorem 3 gives a new expression of the delay, based on the real arrival curves, Ai. In the model considered, the exact behaviour is unknown, but it is possible to bound the difference between the ideal systems and the real system. Then, the corollary 2 generalises it by considering bounds on phases, and Property 5 generalises it to variable frame sizes.

Theorem 3. Let S be a non-preemptive SP n-server taking the natural order as priority order, offering an aggregate minimal strict service of curve β. Then an upper bound on the delay of the flow i is:

hDev   Ai, ( j≤i Aj) * [β -l max N P ] + - j<i Aj   (6)
Proof.

Di(t) = j≤i Dj - j<i Dj ≥ j≤i Dj - j<i Aj ≥ ( j≤i Aj) * [β -l max N P ] + - j<i Aj (7)
Moreover, given f, g, g , three flows. If g ≥ g then hDev(f, g) ≤ hDev(f, g ) And so we can bound hDev(Ai, Di).

Theorem 3, contrary to property 1, bounds the traversal time without using arrival curve and therefore avoids an approximation but it requires to know the exact behavior of the flows.

If phases are unknown but can be bounded, it can be adapted to:

Corollary 2. hDev   A i , ( j≤i A j δ φ(i,j) ) * [β -l max N P ] + - j<i A j * δ φ(i,j)   (8)
Where φ(i, j) is the phases between node which sends Ai and node which sends Aj.

Proof. We will use the clock of the node which sends Ai as reference. Considering a flow Aj we can bound it by: Aj δ φ(i,j) ≤ Aj ≤ Aj * δ φ(i,j) So we can deduce:

Di ≥ ( j≤i Aj) * [β -l max N P ] + - j<i Aj ≥ ( j≤i Aj δ φ(i,j) ) * [β -l max N P ] + - j<i Aj * δ φ(i,j)
And then bound hDev(Ai, Di) in the case of phases.

Each periodic flow is also characterized by its maximal size of frame, but the real size can vary. These variations can be due to different payload or to bit stuffing (bit stuffing is the insertion of non information bits into data and it can increase by 25% at most the packet length). Property 5. Let A ≥ A, Theorem 3 can be extended to:

hDev(Ai, Di) ≤ hDev(Ai, ( j≤i Aj) * [β -l max N P ] + - j<i Aj) (9)
Proof. First we will only consider that Ai is known but that Aj can only be bounded if j < i.

(Ai + j<i Aj) * [β -l max N P ] + - j<i Aj (10) = (Ai + j<i (Aj -Aj + Aj)) * [β -l max N P ] + - j<i Aj (11)
and using property 2,

≤ ( j≤i Aj) * [β -l max N P ] + + j<i (Aj -Aj) - j<i Aj (12) ≤ ( j≤i Aj) * [β -l max N P ] + - j<i Aj (13)

≤ Di

Then, Di is smaller than the expression of eq. ( 10), leading to

hDev(A i , D i ) ≤ hDev(A i , (A i + j<i A j ) * [β -l max N P ] + - j<i A j )
Let us now decompose the i-th flow: A i = A i + A i and D i = D i + D i , A i and A i have the same priority. For flows of same priority the server apply a FIFO policy. Using the result of [START_REF] Boudec | Network calculus: a theory of deterministic queuing systems for the internet[END_REF]:

hDev(A i , D i ) ≤ hDev(A i , D i )
We can deduce:

hDev(A i , D i ) ≤ hDev(A i , ( j≤i A j ) * [β -l max N P ] + - j<i A j )
4.3 Accurate arrival curve for flows with phases Notation 2. Let A1, .., An be a set of flows. In the following, α1..i will denote an arrival curve for the aggregate flow i j=1 Aj.

Notice that i j=1 αj is one particular arrival curve of the aggregate flow, but the core of this paper will be to find better arrival curve, modelling offsets and phases between flows.

To do so, we first simplify the problem by considering a common global hyper-period, i.e. decomposing each flow as a sum of sub-flow with the same hyper-period and specific offsets.

A = n A i=1 ν T,h A * δ O A +(i-1)T A = n A i=1
Ai. And so we can define our functions as the sum of TAB-periodic functions.

Obviously Lemma 1 can be extend to more than two flows. We have seen with Lemma 1 that we can consider an unique common period for different flows so we will consider that all the messages are sent periodically every T seconds. First let us suppose that the synchronization is perfect, for each message m its offset Om is known (we can consider that Om < T ), meaning that the message m will be sent on predefined instants: the first time at time Om and then every T . We will also consider that at this point we are working on the maximal frame size. So the message m can be associated with the flow M (see Figure 7 Property 6. Let A1, . . . , An a set of periodic flows, such as Ai = ν T,h i * δO i with ∀(i, j) ∈ {1, n}, i ≤ j ⇒ 0 ≤ Oi ≤ Oj < T . We define On+1 such as On+1 = O1 + T . Ai the aggregate flow, such as

Flow with known offset difference

∀j ≤ n, ∀t ∈ [0, T [, Aj(t) = Aj(O next(t) ) (14) Proof. Let t ∈ [0, T [, j < next(t) ⇒ Oj + T ≥ Onext ≥ t > Oj j ≥ next(t) ⇒ t ≤ O next(t) ≤ Oj ⇒ Aj(t) = Aj(Oj + T ) = A(O next(t) ) Aj(t) = Aj(Oj) = A(O next(t) ) Theorem
Ai = ν T,h i * δO i with ∀(i, j) ∈ {1, n}, i ≤ j ⇒ 0 ≤ Oi ≤ Oj < T .
We define On+1 such as On+1 = O1 + T . Then α is an arrival curve of A with:

α(d) = max 1≤i≤n (A(Oi + d) -A(Oi)) ( 15
)
And α is the best arrival curve of A.

Proof. Let t ∈ R * such as t = kT + t0 with 0 ≤ t0 < T and k ∈ N and d ∈ R * , then

A(t + d) -A(t) = A(t0 + d) -A(t0) = A(t0 + d) -A(O next(t 0 ) ) ≤ A(O next(t 0 ) + d) -A(O next(t 0 ) ) ≤ max 1≤i≤n (A(Oi + d) -A(Oi)) So max i (A(Oi + d) -A(Oi)) is an arrival curve of flow A.
Using property 4 we can deduce that it is the best arrival curve.

Property 7. Let A1, . . . , An a set of periodic flows, such as Ai = ν T,h i * δO i with ∀(i, j) ∈ {1, n}, i

≤ j ⇒ 0 ≤ Oi ≤ Oj < T . i (A(Oi + d) -A(Oi)) = max i n j=1 ν T,h j * δ O j -O i + O i -O j T T (d) (16)
Proof. Let i ≤ n, A(Oi + d) -A(Oi) = n j=1 Aj(Oi + d) - Aj(Oi) • j < i : Aj(Oi + d) -Aj(Oi) = hj d+(O i -O j ) T -hj = ν T,h j * δO j -O i (d)-hj = ν T,h j * δ O j -O i + O i -O j T T (d) Because 0 ≤ Oj ≤ Oi < T ⇒ O i -O j T T = T • j ≥ i : Aj(Oi + d) -Aj(Oi) = hj d+(O i -O j ) T = ν T,h j * δO j -O i (d) = ν T,h j * δ O j -O i + O i -O j T T (d) Because 0 ≤ Oi ≤ Oj < T ⇒ O i -O j T T = 0
We can deduce that:

n j=1 Aj(Oi+d)-Aj(Oi) = n j=1 ν T,h j * δ O j -O i + O i -O j T T (d)
This last property gives us the arrival curve for a set a flows considering that the synchronization is perfect.

Flows with bounded offset difference

We will now consider that the synchronization is not perfect, there is a phase between each calculator's clock and the global time. However this phase can be bounded, we can consider the offset difference between all the messages.

= ν T,h i * δO i with ∀(i, j) ∈ {1, n}, -T < Φij ≤ Oi -Oj ≤ Φ ij < T ( 17
)
Then α is an arrival curve of A with:

α = max i n j=1 ν T,h i * δ max(0,-Φ ij + Φ ij T T ) (18)
Proof. Theorem 4, have to be adapted since Oi and Oj are unknown. By definition of Φ ij , Φij, Oj-Oi+

O i -O j T T ≥ -Φ ij + Φ ij T
T . Moreover, the ceiling function ensures

x T ≥ x T , to -x + x T T ≥ 0, leading to Oj -Oi + Oi -Oj T T ≥ max(0, -Φ ij + Φij T T )
Then we deduce that:

δ O j -O i + O i -O j T T ≤ δ max(0,-Φ ij + Φ ij T T )
And so to conclude:

α(d) = max i n j=1 ν T,h j * δ O j -O i + O i -O j T T (d) ≤ max i n j=1 ν T,h i * δ max(0,-Φ ij + Φ ij T T ) (d)
We have found a function bigger than an arrival curve so it is also an arrival curve.

EXPERIMENTAL RESULTS

We evaluate in this section the actual gain that can be expected using frame offsets in a network where local clocks have bounded phases. We consider a CAN network where the phases between nodes successively amounts to 0ms, 1ms, 5ms and finally 10ms, and compare the bounds on the frame communication latencies obtained with different timing analyses.

Configuration under study

The CAN network under study runs at 250kbit/s and is made of 10 nodes. The periods of the frames are uniformly chosen from the set {20, 50, 100, 200, 500, 1000}ms and the data payloads sizes are drawn at random according to a uniform distribution between 1 and 8 bytes. Figure 8 shows for the set of frames the number of flows per period and payload, as indicated by the size of circle and the number inside. The total load on the network is 35%, and there is no difference between stations, each of them generating similar loads. The frame offset assignment that is used in this study is the SOPA algorithm available in the RTaW-Pegase software from the company RTaW, whose results have been used as a benchmark to validate the analysis developed in this and the gain obtained with bounded phases. The SOPA algorithm has been chosen because our experience is that it consistently outperforms the few other offset algorithms available [START_REF] Grenier | Pushing the limits of CAN-scheduling frames with offsets provides a major performance boost[END_REF][START_REF] Grenier | Near-Optimal Fixed Priority Preemptive Scheduling of Offset Free Systems[END_REF], and thus provides us with an estimate of the best possible gains that can be achieved in practice with offsets.

Methodology and timing analyses

Hereafter we evaluate the impact of the phases between the nodes by comparing bounded phases, ideal global clock, and no time synchronization at all between nodes, like usually on CAN. This will be performed by gradually increasing the upper bound of the phases between the nodes. To gain a better understanding, we consider that the bound on the phases is identical for all nodes and symmetrical. For example, a phase of 1ms means that the clock of a station can be in advance, or late, by at most 1ms compared to any other station in the network. In other words from the point of view of a given node, the other nodes can send their frames over a 2ms window that is centered on the expected date.

There are several ways to bound the communication latencies in systems with bounded phases, all of which making use of results that have been introduced before:

• Method 1 (M.1): Use of residual service (5a), and calculations of arrival curves with phases (Theorem 5).

• Method 2 (M.2): Bound by size of busy window (5b), and again calculations of arrival curves with phases (Theorem 5).

• Method 3 (M.3): New method to bound delay without using arrival curves, relying on Property 5.

We then include on the comparison results obtained with local and global clocks. For local clocks:

• Method 4 (M.4): Use of the algorithm published in [START_REF] Yomsi | Controller Area Network (CAN): Response time analysis with offsets[END_REF] that is available in RTaW-Pegase. This algorithm is exact (except for an overestimation of the blocking factor of CAN frame) and takes into account frame offsets but consider that the phases between nodes are unknown.

For global clocks:

• Method 5 (M.5): From the network point of view, a single node or several with perfectly synchronized clocks leads to the exact same workload submitted to the network. The same algorithm as in method 4 is thus used, but this time with a single node on the bus.

And in order to be able to evaluate the gain brought by offsets:

• Method 6 (M.6): We compute the exact response time, presented in [START_REF] Davis | Controller Area Network (CAN) schedulability analysis: Refuted, revisited and revised[END_REF], that does not consider offsets.

• Method 7 (M.7): Network calculus without offsets

(α1..i-1 = i-1 k=1
α k ) with (5a) is used.

Experimental results

For each experiment, the results obtained with the different methods that are compared are reported in a figure and a table. The table contains the average and maximal delay bounds obtained with the different methods with the flows grouped by priority ranges.

Global clock (null phases).

Let us start by the global clock configuration, i.e. null phases between stations, method 5 in Figure 9. This configuration allows to evaluate the extent to which bounded phases degrade the performances in terms of communication latencies. Table 1 summarizes the results with the different methods.

When phases are null, (M.3) is as good as (M.5).

The results with (M.1)(M.2) are slightly greater than with (M.5) but still close: the difference between the average delay bound with (M.2) and (M.5) is less than 0.5ms (cf. Table 1), which is the transmission time of a frame with a payload of 8 bytes on a 250kbit/s bus.

This experiment also shows the gain with global clock (methods M.1,2,3,5) with regard to local clocks (M.4). As expected from other results in the literature [START_REF] Grenier | Pushing the limits of CAN-scheduling frames with offsets provides a major performance boost[END_REF], the latter is also much better than without any synchronization at all, i.e. methods (M.6) and (M.7). Without offsets, the latency increases almost linearly. The discontinuity at 20ms can be simply explain: it is the shortest transmission period, and from this point on frames can be delayed by more than one frame instance by higher priority flow.

Lastly, we note that network calculus without offset (M.7) leads to the exact same response times as with the analysis in [START_REF] Davis | Controller Area Network (CAN) schedulability analysis: Refuted, revisited and revised[END_REF] (M.6). In the following we will refer to these methods as (M.6,7).

Bounded phases of maximum 1ms.

In the next experiment the system configuration stays the same, the only difference is that we suppose that now our phases are unknown but bounded by 1ms. The results for this configuration are shown in Figure 10. This choice of phases bound has been done because TTCAN level-1 have a precision in the order of 1 ms [START_REF] Rodrıguez-Navas | Clock synchronization in CAN distributed embedded systems[END_REF]. In these setup, the analysis for global clock (M.5) cannot be used anymore. First consider (M.3): for high priority messages (0-35), it is still the method which provides the minimal bounds, however for flows of lower priorities it returns a bound of variable accuracy. Sometimes this is the best we have, but sometimes it is even worst than when offsets are not used (M.6,7). In the tables, in order to avoid extreme values which have no interest, we consider the minimum between (M.3) and (M.6,7).

Then let us consider (M.1) and (M.2): this phase of 1ms increases the delay bound by an average of 0.5ms compared to a global clock. It should be noted that these bounds are still well below those of (M.4), which is the exact offset analysis without accounting for inter-node synchronization.

Bounded phases of maximum 5ms.

In the next experiment we increase significantly our bound on phases to 5ms (see Figure 11 for the results obtained with this configuration). First let us consider (M.3): the previous observation is confirmed, for very high priority frames (0-12), this is still the method which provides the most accurate bound, however for lower priority frames, the bounds returned are unusable.

Then consider (M.1) and (M.2): compared to the bounds given by (M.4), accouting for bounded phases permits to reduce the bound by almost 40%, which is very significant for many real-time distributed applications.

Bounded phases of maximum 10ms.

In the last experiment presented in this paper the scheduling still remains the same, but now phases are bounded by 10ms (see Figure 12 for the results obtained with this configuration). We remind that phases bounded by more or less 10ms implies a transmission window of length 20ms. The smallest period of emission of our network is 20ms, meaning that offsets have no more influence on flows with the smallest period.

It is confirmed that (M.3) should not be used for too large phases, even if it still gives the best result for 8 highest priority messages.

This experiment shows that even if (M.2) very often return better results than (M.1), it is not always the case. Moreover, if these methods most of the time give better results than considering fully unsynchronized nodes. However, especially for the lower priority flows, sometimes the bounds calculated using these functions is larger than the one of (M.4).

CONCLUSIONS AND DISCUSSION

It has been largely evidenced in the litterature that offsets, by reducing contentions, are an efficient way to improve the response times of tasks and messages in periodic realtime systems. In distributed systems one distinguish two main approaches: either the global clock paradigm, leading to a globally time-triggered system, and requiring synchronization mechanisms with a precision much smaller than the frame sending time, or resorting to local clocks only and thus without the requirement of inter-node synchronization.

In this paper, we investigate a trade-off: local clocks with bounded phases between nodes and apply it to CAN networks. Since, to the best of our knowledge, no results were readily available in the litterature for estimating the worstcase traversal times in such systems, the first contribution of this work is an extension of an existing network calculus analysis to handle such systems. Three methods to compute such bound have actually been proposed. On the case study, it appears that they are incomparable, i.e. there is no method that always outperforms the others, even if it appears that method (M.3) returns more accurate results when the phases are small (1ms in our example), and that method (M.2) give the best results, in average.

For the sake of understandability, some simplifications have been done: the traffic associated to the synchronization mechanism has not been considered, no event-triggered traffic is considered, and frame transmissions do not suffer jitters. These are not limitations of the method, but methodological choices. Adding other kind of data flow can be done by using the Theorem 2 on non-preemptive static priority. Considering a data flow of priority i, if it exists some sporadic flow of higher priority frames with arrival curve α s , we have to replace β by β -α s in the equations used to compute the delay on the flows with phases. And the same can be done for the flow dedicated to clock synchronization. Adding a jitter d to an arrival curve α is just a deconvolution by δ d . Then, modeling such a jitter in each node can be done by simply adding a δ d in eq. [START_REF] Tindell | Analysing Real-Time Communications: Controller Area Network (CAN)[END_REF].

The experiments also have produced some insights in the beneficial impact of bounded phases. Comparing bounded phases to the global clock solution, it appears that a phase of 1ms leads to limited additional latencies, smaller than 3ms for the 250kbit/s network used in the experiments. And comparing phases to the local clock, a phase of 5ms already provides an average gain of around 40%. These are very promising results.

In the future, we would like first to enhance the analysis method itself: currently, some information is lost when a flow is decomposed as a sum of sub-flows, in order to have a common period of all sub-flows. We would like also to compare the analytic results to simulation results. But the main work is about further evaluating the benefits of using bounded phases: how does the gain change with respect to the bus load or with respect to the period distribution? Could offset assignment algorithm tailored to systems with phases lead to substantial improvements?

Figure 2 :

 2 Figure 2: Impact of phasing between two nodes.

Figure 3 :Figure 4 :

 34 Figure 3: Latency and stair functions

Figure 5 :Figure 6 :

 56 Figure 5: Backlog and delay

Definition 7 .

 7 (Min-plus minimal service curve) A server S offers a min-plus minimal service curve β iff for all input/output A,D D ≥ A * β These two notions of service curve are not equivalent but are both of interest and are related: a server offering a minimal strict service β also offers the minimal min-plus service β.

Lemma 1 .

 1 (Common period) Let A and B two periodic flows, such as A = ν T A ,h A * δO A and B = ν T B ,h B * δO B with TA, TB ∈ N. Then we can define TAB and (nA, nB) ∈ N such as TAB = nATA = nBTB. Now A (resp. B) can be express as:

Figure 7 :

 7 Figure 7: Example of flow M.

  ): M = ν T,hm * δO m Definition 8. Let O1, . . . , On ∈ [0, T [ such as 0 ≤ O1 ≤ . . . ≤ On < T and On+1 = O1 + T then we define for all t ∈ [0, T [: next(t) = min 1≤i≤n+1 {i Oi ≥ t}. next(t) is the next flow that will send a message after time t.

4 .

 4 (Periodic flow with offset) Let A1, . . . , An a set of periodic flows, A = n i=1

Theorem 5 .

 5 (Periodic flow with bounded offset difference) Let A1, . . . , An a set of periodic flows, A = n i=1 Ai the aggregate flow, such as Ai

Figure 8 :

 8 Figure 8: Distribution of the size of the data payloads and periods for the tested configuration.

Figure 9 :

 9 Figure 9: Delay bounds with null phases.

Figure 10 :

 10 Figure 10: Delay with phases bounded by 1ms.

Figure 11 :

 11 Figure 11: Delay bounds with phases ≤ 5ms.

Figure 12 :

 12 Figure 12: Delay bounds with phases ≤ 10ms.

Table 1 :

 1 Delay bounds when phases are null.

	Priorities	2-16	17-31	32-46	47-61	2-61
	Average delay bounds (in ms)
	Phases: M.1	1.47	1.49	1.82	1.94	1.68
	Phases: M.2	1.04	1.05	1.40	1.54	1.25
	Phases: M.3	0.97	0.98	0.94	0.91	0.95
	Local clocks: M.4	3.00	5.84	8.04	10.88	6.94
	Global clock: M.5	0.97	0.98	0.94	0.91	0.95
	No offsets: M.6,7	4.80	11.23	17.96	27.11	15.27
	Maximum delay bounds (in ms)
	Phases: M.1	1.54	1.62	2.50	2.46	2.50
	Phases: M.2	1.04	1.08	1.54	1.54	1.54
	Phases: M.3	1.04	1.08	1.08	1.04	1.08
	Local clocks: M.4	4.26	6.64	10.04	12.28	12.28
	Global clock: M.5	1.04	1.08	1.08	1.04	1.08
	No offsets: M.6,7	7.84	14.46	23.66	30.30	30.30

Table 2 :

 2 Delay bounds with phases ≤ 1ms.

	Priorities	2-16	17-31	32-46	47-61	2-61
	Average delay bounds (in ms)
	Phases: M.1	1.47	1.85	2.38	3.09	2.20
	Phases: M.2	1.04	1.45	1.76	2.46	1.68
	Phases: min(M.3,M.7)	0.97	0.98	3.86	13.88	4.93
	Local clocks: M.4	3.00	5.84	8.04	10.88	6.94
	No offsets: M.6,7	4.80	11.23	17.96	27.11	15.28
	Maximum delay bounds (in ms)
	Phases: M.1	1.54	2.12	2.88	3.72	3.72
	Phases: M.2	1.04	1.58	1.96	2.88	2.88
	Phases: min(M.3,M.7)	1.04	1.08	17.3	30.3	30.3
	Local clocks: M.4	4.26	6.64	10.04	12.28	12.28
	No offsets: M.6,7	7.84	14.46	23.66	30.30	30.30

Table 3 :

 3 Delay bounds with phases ≤ 5ms.

	Priorities	2-16	17-31	32-46	47-61	2-61
	Average delay bounds (in ms)
	Phases: M.1	1,77	3,17	4,91	7,02	4,22
	Phases: M.2	1,39	2,82	4,57	6,48	3,82
	Phases: min(M.3,M.7)	2,70	10,58	17,96	27,12	14,59
	Local clocks: M.4	3,01	5,85	8,04	10,88	6,94
	No offsets: M.6,7	4,80	11,23	17,96	27,12	15,28
	Maximum delay bounds (in ms)
	Phases: M.1	2.38	3.80	6.40	8.00	8.00
	Phases: M.2	1.88	3.26	5.68	7.66	7.66
	Phases: min(M.3,M.7)	7.84	14.46	23.66	30.30	30.30
	Local clocks: M.4	4.26	6.64	10.04	12.28	12.28
	No offsets: M.6,7	7.84	14.46	23.66	30.30	30.30

Thanks

The authors would like to thank Nicolas Navet, for its great help on a preliminary version of this work.

This work has been partially funded by French agency DGA, under project IREHDO 2.


Jean-Baptiste Chaudron
email: jean-baptiste.chaudron@isae.fr

David Saussié
email: d.saussie@polymtl.ca

Pierre Siron
email: pierre.siron@isae.fr

Martin Adelantado
email: martin.adelantado@onera.fr

How to solve ODEs in real-time HLA distributed simulation

Keywords: Simulation, ODEs, HLA, RTI, real-time, scheduling

In the context of the Research Platform for Embedded Systems Engineering (PRISE) Project, we are developing and maintaining a complete aircraft flight simulation using the High Level Architecture (HLA), an IEEE standard for distributed simulation. This complex distributed simulation is composed of different distributed HLA simulators (e.g., Flight Dynamics, Sensors), whose dynamic behaviors are implemented as Ordinary Differential Equations (ODEs). The resolution of these equations is done, locally for each simulator, by numerical integration with methods like Euler or Adams-Bashforth. The global behavior of this distributed simulation, where each component runs its own local resolution, is a key challenge. The main problem is to ensure the global simulation consistency and, in particular, the specific data flows between components with the correct temporal real-time behavior. This paper specifically addresses the problem of solving ODEs over an HLA distributed architecture and offers a complete study (specifications, implementation and validation) where several theoretical concepts and methods are discussed.

Introduction

The main goal of the Research Platform for Embedded Systems Engineering (PRISE) project is to provide students and researchers with a platform for the study, evaluation and validation of new embedded system concepts, architectures and techniques through a dedicated hardware and software environment. Modern embedded systems become more and more complex with an increasing number of both components and interactions between them and the use of simulation is essential for relevant studies. In this context, we implemented from scratch and we are currently maintaining an aircraft component-based simulation composed of several distributed simulators [1], [2]. Each simulator represents a specific part of the aircraft (e.g., engines, actuators, control devices, flight control laws), its dynamic in the environment (flight dynamics) or the environment itself. An overview of this distributed simulation is depicted in Figure 1. This flight simulation, called SDSE (French acronym for Distributed Simulation of Embedded Systems) has been implemented by using the HLA IEEE standard [START_REF]IEEE Std 1516 TM -2010, IEEE Standard for Modeling and Simulation (M&S) High Level Architecture (HLA) -Framework and Rules[END_REF], [START_REF][END_REF], [START_REF][END_REF] which provides a well-known simulation framework for composability, maintenance and reusability for this complex simulation. Moreover, HLA provides time management mechanisms [START_REF] Fujimoto | Time Management in the High Level Architecture[END_REF], guaranteeing a consistent global logical time throughout the whole simulation, which is one of the main benefits of this simulation standard. Previous works showed that these time management mechanisms could be used to ensure the good behavior during run time, especially in real-time simulations [START_REF] Phillips | Time Management Issues and Approaches for Real Time HLA based simulations[END_REF], [START_REF] Fujimoto | Repeatability in real-time distributed simulation executions[END_REF], which is typically an aspect of the SDSE simulation. From the software point of view, we are using the open source CERTI [START_REF] Siron | CERTI : an open Source RTI, why and how[END_REF] Run Time Infrastructure (HLA middleware implementation) because we have a complete control on its implementation. We have shown that it is suitable for interconnection of simulators (i.e., HLA federates) with short computation and communication cycles [START_REF] Noulard | Running Real Time Distributed Simulations under Linux and CERTI[END_REF], [START_REF] Chaudron | HLA High-Performance and Real-Time Simulation Studies with CERTI[END_REF]. Additionally, we have also investigated the usage of CERTI time management techniques and implementations for real-time simulations [START_REF] Chaudron | Design and model-checking techniques applied to real-time RTI time management[END_REF] and especially for SDSE.

Scheduling a real-time federation (i.e., an HLA simulation) consists of two complementary facets that must be distinguished, namely, the temporal scheduling and the functional scheduling.

1) The temporal scheduling belongs to scheduling theory. Here, the objective is to formally verify that the distributed simulation will run while ensuring compliance with the real-time constraints according to its specifications. This aspect has been addressed in [2] and is not the purpose of the present study.

2) The functional scheduling ensures that the sequence of the communications and computations between each simulator will provide global simulation results that will always be fair and relevant according to simulation model semantics.

HLA time management mechanisms are providing a good baseline for this purpose. In the SDSE simulation, each federate implements a dedicated specific algorithm and some of these models are using Ordinary Differential Equations (ODEs) [START_REF] Cooke | Ordinary Differential Equations[END_REF]. Therefore, the distributed resolution of these ODEs is quite complex and needs to be clearly analyzed; this is the purpose of the paper.

Also, the notion of "time" within a real-time simulation is an important notion which needs to be clearly stated for consistent temporal and functional scheduling study. Three different concepts can be associated to the notion of time, each one describing a specific aspect [START_REF] Fujimoto | Parallel Discrete Event Simulation[END_REF]:

1) The physical time, denoted by t PHY , is the reference time of the physical system under study (the source system that we want to model and to simulate).

2) The simulated time, denoted by t SIM , is a representation of the physical time within the simulation. It is a set of ordered values representing different time instants for the modeled physical system. 3) The absolute time or wall-clock time, denoted by t WCT , is the time elapsed while executing the simulation. It can be measured by a hardware clock, like the CPU clock. In our real-time simulations, with the usage of conservative time management methods, these times might be equal (with respect to a certain unit and a certain precision). Therefore, during the run of a simulation, we have:

t PHY = t SIM = t WCT (1)
For example, if the millisecond is considered as the reference unit, 1 millisecond of the real physical system is equal to 1 unit of logical HLA simulation time, which is also equal to 1 millisecond of execution time.

To resume, this paper focuses on the study of functional scheduling for HLA simulations (i.e., federations) using HLA time management semantics where some simulators (i.e., federates) are modeled by ODEs. Based on our experience with the SDSE aircraft simulation development, we present here a detailed analysis from initial ODEs characterization to distributed resolution description. The paper is organized as follows:

• Section 2 describes ODEs characterization for our real-time simulation;

• Section 3 outlines the resolution of ODEs for federation with a communication sequence;

• Section 4 presents the resolution of ODEs for federation with a communication loop;

• Section 5 summarizes these concepts and shows an illustration for the SDSE aircraft simulation.

ODEs Characterization 2.1 State-Space Description

The dynamic behavior of a system can usually be modeled by a set of Ordinary Differential Equations (ODEs) and/or Partial Derivative Equations (PDEs) derived from different laws of physics. Under the action of external excitation (inputs) and/or initial conditions, the system evolves accordingly, and one can then monitor signals of interest (outputs), representative of the system behavior. In this paper, we suppose that the federates simulate systems governed by ODEs and described by state-space models of the form:

ẋ(t) = f (x(t), u(t), t) (2) y(t) = g(x(t), u(t), t) (3) x(0) = x 0 (4)
where x = [x 1 , . . . , x n ] ∈ R n denotes the state vector composed of n state variables x i , u = [u 1 , . . . , u m ] ∈ R m the input vector, y = [y 1 , . . . , y p ] ∈ R p the output vector and ẋ(t) := dx(t)/dt. The state equation (2) is a set of n coupled first-order ODEs, in which each time derivative ẋi is expressed in terms of the state variables x 1 , . . . , x n and the system inputs u 1 , . . . , u m . The state equation represents the dynamic behavior of the system under the action of the inputs and the initial conditions defined by equation ( 4). The output equation (3) yields any signal of interest depending on state variables and/or input variables. In the general case, f (.) et g(.) are non-linear vector functions.

From an HLA point of view, a federate subscribes to object attributes that will constitute the input vector that drives the state equation, and publishes object attributes yielded by the output equation.

Numerical methods for ODEs

Numerical solutions of ODEs can be performed by a large variety of methods depending on the nature of the problem and the desired accuracy [START_REF] Cooke | Ordinary Differential Equations[END_REF]. Starting from an initial point, a numerical method takes a short step forward in time to find the next point, and so on. Among all the available numerical methods, one can find single-step methods (e.g., Euler method), Runge-Kutta methods, and multistep methods (e.g., Adams-Bashforth and Adams-Moulton methods). These methods imply the choice of an adequate time step ∆t, which can have a significant impact on the accuracy of the solution; they all lead to recurrence equations. For instance, the forward Euler method applied to original state equation (2) leads to the recurrence equation ( 5):

x k+1 = x k + ∆t f (x n , u n , t n ) (5)
where n denotes here the iteration, t n = n∆t the corresponding time, and

x n = x(t n ), u n = u(t n
) the values of x and u at t n . The Euler method is an explicit method in the sense that it calculates the next state value x n+1 from previous values x n and u n . By extension, explicit Adams-Bashforth methods, can be written using s previous values as in equation ( 6):

x n+1 = x n + ∆t s k=0 b k f (x n-k , u n-k , t n-k ) (6)
where the b k are chosen such that the method has order s. As an example, the second-order and third-order Adams-Bashforth methods are respectively given by equations ( 7)-( 8):

(AB2)

x n+1 = x n + ∆t 3 2 f (x n , u n , t n ) - 1 2 f (x n-1 , u n-1 , t n-1 ) (7)
(AB3)

x n+1 = x n + ∆t 23 12 f (x n , u n , t n ) - 4 3 f (x n-1 , u n-1 , t n-1 ) + 5 12 f (x n-2 , u n-2 , t n-2 ) (8)
On the opposite, there are implicit methods such as the backward Euler method described here in equation ( 9):

x n+1 = x n + ∆t f (x n+1 , u n+1 , t n+1 ) (9)
where one needs to solve a non-linear algebraic equation at every time step in order to produce the next value x n+1 ; this is usually achieved by using fixed point iteration or Newton-Raphson method. The computational cost is obviously higher, but implicit methods are usually more stable than their explicit counterparts. Nevertheless, the implicit methods are rejected in this HLA simulation context as they necessitate the value of the input u n+1 at the next step n + 1, i.e., in the future.

Direct feedthrough

The output equation ( 3) is an algebraic equation that produces the output vector y(t) containing the signals of interest. Generally, the output vector y(t) may depend on both the state vector x(t) and the input vector u(t), and eventually, explicitly on the time t, as in equation ( 3). Contrary to the state equation ( 2) which requires numeric solvers as shown above, the discretization of the output equation, described in [START_REF] Noulard | Running Real Time Distributed Simulations under Linux and CERTI[END_REF], is straightforward:

∀t, y(t) = g (x(t), u(t), t) ⇒ y n = g (x n , u n , t n ) (10)
Nevertheless, special attention must be paid to the presence, or not, of the input vector u(t) (or u n ) in the output equation.

When the system output y(t) depends on the input vector u(t) (or on some of the inputs u i ), it is called direct feedthrough.

As it will be presented in Sections 3 and 4, this can reveal particularly intricate when simulating a distributed system. In the case where there is no direct feedthrough, the output equation ( 11) is:

∀t, y(t) = g (x(t), t) ⇒ y n = g (x n , t n ) (11)

Typical execution within a federate

The choice of the numerical method and the time step ∆t depends on the characteristics of the system but is not discussed here. The discretization of a state-space model with an explicit linear multi step method leads to the following discreet system (12):

   ẋ(t) = f (x(t), u(t), t) y(t) = g(x(t), u(t), t) x(0) = x 0 ⇒    x n+1 = x n + ∆t s k=0 b k f (x n-k , u n-k , t n-k ) y n = g(x n , u n , t n ) x(0) = x 0 (12)
At iteration n, knowing x n and u n , the output equation is first executed to provide y n , then the state equation yields the state vector update x n+1 for the next iteration. Note that depending on the numerical methods, previous values x n-k and u n-k must be stored for the recurrence equation. For the sake of simplicity but without loss of generality, we will consider in the following linear time-invariant (LTI) systems, i.e, systems described by linear differential and algebraic equations with constant coefficients. The original equations ( 2)-( 3) can be rewritten as linear equations ( 13)-( 14):

ẋ(t) = Ax(t) + Bu(t) (13) y(t) = Cx(t) + Du(t) (14)
where A, B, C and D are respectively n × n, n × m, p × n and p × m constant matrices. For example, the application of the forward Euler method yields equations ( 15)-( 16):

x n+1 = x n + ∆t (Ax n + Bu n ) (15) y n = Cx n + Du n (16)
In the following, the influence of direct feedthrough on simulation, is illustrated on a first example with a federation involving federates with a communication sequence. Thereafter, we pursue the study on a federation that presents communication loop between federates. These studies will lead to the establishment of simple rules for the automatic implementation of consistent functional scheduling for these types of federation.

Federation with a Communication Sequence

Let an HLA federation composed of 4 federates which are communicating in series as illustrated in Figure 2. We study the impact of direct feedthrough in some of the federates. Federate 0 is a federate which provides a signal u 1,n (every cycle n). For the sake of simplicity but without loss of generality, it is assumed that the numerical solution of the state-space models representing the dynamic behavior of the next 3 federates is described in equations ( 17)-( 18):

x i,n+1 = A i x i,n + B i u i,n (17)
y i,n = C i x i,n (18)
where i ∈ 1, 2, 3 denotes the federate index. The federates are linked in series such that the output of Federate i is entirely the input of Federate i + 1, i.e., ∀i ∈ 2, 3, u i,n = y i-1,n . We can see that these federates do not exhibit direct feedthrough in their respective output equation, and therefore no special care is needed for the execution of the federates. This gives the following global system with equations ( 19)-(20):

  x 1,n+1 x 2,n+1 x 3,n+1   =   A 1 0 0 B 2 C 1 A 2 0 0 B 3 C 2 A 3     x 1,n x 2,n x 3,n   +   B 1 0 0   u 1,n (19)
  y 1,n y 2,n y 3,n   =   C 1 0 0 0 C 2 0 0 0 C 3     x 1,n x 2,n x 3,n   (20)
From the initial conditions [x 1,0 x 2,0 x 3,0 ] , we can calculate the state and output for every iteration step [y 1,n y 2,n y 3,n ] . We assume that, from initial moment t = 0, all states from Federates 1, 2 and 3 are x i,0 , ∆t is the time step between cycle and lk ≤ ∆t is the lookahead. Federate 0 can be implemented according the pseudo-code described in Table 1 and the other Federates as in Table 2.

In the remainder of the paper, for clarity sake, acronyms will be used to denote well-known HLA services and callbacks:

• UAV stands for updateAttributesValues(...);

• RAV for reflectAttributesValues(...);

• TAR for timeAdvanceRequest(...);

• NER for nextEventRequest(...);

• and TAG for timeAdvanceGrant(...). Also, we consider the case where each federate is regulator and constrained for HLA time management mechanisms and we do not consider the zero-lookahead case (i.e., lk = 0).

O u 1,0 u 1,0 O u 1,1 u 1,1 O u 1,2 u 1,2 O u 1,3 u 1,3 LT=0 LT=10 LT=20 LT=30 LT=40 UA V( 1) TA R( 10 ) TA G ( 10 )

Direct feedthrough case

We still consider the test case described above but we are assuming now that Federates 2 and 3 contain direct feedthrough relation. The equation system for Federates 1, 2 and 3 is depicted in (21)-( 22)-( 23)-( 24):

x 1,n+1 = A 1 x 1,n + B 1 u 1,n x 2,n+1 = A 2 x 2,n + B 2 u 2,n x 3,n+1 = A 3 x 3,n + B 3 u 3,n (21) y 1,n = C 1 x 1,n y 2,n = C 2 x 2,n + D 2 u 2,n y 3,n = C 3 x 3,n + D 3 u 3,n (22)
  x 1,n+1 x 2,n+1 x 3,n+1   =   A 1 0 0 B 2 C 1 A 2 0 B 3 D 2 C 1 B 3 C 2 A 3     x 1,n x 2,n x 3,n   +   B 1 0 0   u 1,n (23)
  y 1,n y 2,n y 3,n   =   C 1 0 0 D 2 C 1 C 2 0 D 3 D 2 C 1 D 3 C 3 C 3     x 1,n x 2,n x 3,n   (24)
Reminding that ∀i ∈ 1, 2, 3, u i,n = y i-1,n we intuitively feel that direct feedthrough will cause problems. Let's consider Federate 2 implemented as described in Table 2. At phase 2.(a), the federate should calculate its output value y 2,n from x 2,n and u 2,n = y 1,n . However for n > 0, according to this pseudo-code implementation, Federate 1 doesn't know about y 1,n but only about y 1,n-1 . This means that it knows values from the previous step but not the one from the same step which has just been calculated per federate 1. Therefore for the same cycle, Federate 2 needs to know its input at instant n to calculate y 2,n .

As Federates are communicating in sequence, this problem is even worth for Federate 3 which must wait for the update of the output from Federate 2, which is itself depending on Federate 1. So there is a functional schedule principle to build in order to ensure consistency between publications and acquisitions.

We solve this problem by combining the usage of NER with the usage of TAR and we also define a dedicated and correct value for the lookahead. NER service allows to request an advance in the HLA logical time but unlike the TAR service, this time advance will be granted either for the required time or for the time of the first available RAV if it is smaller than the requested time. Here the choice of a correct lookahead is a key issue to ensure a correct number of UAV/RAV for the same simulation cycle. Implementation strategies For federate 0 and Federate 1 don't change but, as described in Tables 3 and4, Federates 2 and 3 contain NER usage in their implementation strategies. Figure 4 illustrates the temporal behavior of this new strategy with ∆t = 10 et lk = 1.

1) Initialization phase a) n ← 0 b) Initial choice for x i,

Federation with a Communication Loop

We consider now a federation containing a communication loop between Federates as depicted in Figure 5 where u 2,n = y 1,n and u 1,n = y 2,n .

Standard case

As before, Federate 0 delivers signal u 1,n at every cycle n. We assume that the numerical solutions of state-space models representing the dynamic behaviors of Federate 1 and Federate 2 are given by equations ( 25)-( 26)-( 27)-(28):

x 1,n+1 = A 1 x 1,n + B 1 u 1,n + B 1 u 1,n (25) y 1,n = C 1 x 1,n (26) x 2,n+1 = A 2 x 2,n + B 2 u 2,n (27) y 2,n = C 2 x 2,n (28
O u 1,0 UA V( 11 ) TA R( 20 ) TA G( 20 ) O u 1,1 UA V( 1) TA R( 10 ) RA V ( 1)
TA G( 10) 29)-(30) system:

O y 1,0 S x 1,1 UA V( 11 ) TA R( 20 ) RA V( 11 ) TA G( 20 ) O y 1,1 S x 1,2 u 1,0 u 1,1 NE R( 10 ) TA G( 10 ) RA V( 1) TA G( 1) UA V( 2) TA R( 10 ) O y 2,0 S x 2,1 NE R( 20 ) TA G( 20 ) RA V( 11 ) TA G( 11 ) UA V( 12 ) TA R( 20 ) O y 2,1 S x 2,2 x 2,0 y 1,0 y 1,1 NE R( 10 ) TA G( 10 ) RA V( 2) TA G( 2) UA V( 10 ) TA R( 10 ) O y 3,0 S x 3,1 NE R( 20 ) TA G( 20 ) RA V( 12 ) TA G( 12 ) UA V( 20 ) TA R( 20 ) O y 3,1 S x 3,2 x 3,0 y 2,0 y 2,1 Figure 4. Execution sequence direct feedthrough case Federate 0 Federate 1 Federate 2 u 1,n y 1,n → u 2,n u 1,n ← y 2,n
x 1,n+1 x 2,n+1 = A 1 B 1 C 2 B 2 C 1 A 2 x 1,n x 2,n + B 1 0 u 1,n (29)
y 1,n y 2,n = C 1 0 0 C 2 x 1,n x 2,n (30)
Without direct feedthrough and despite the presence of a loop, the resolution takes place as it is done in the context of a federation with communication sequence. Thus, we can apply the implementation strategies described in Tables 1 and2.

Direct feedthrough case

Now, let's consider a direct feedthrough relation for federate 1 depending on its second input u 1,n , then the model of federate 1 can be described as in equations ( 31)-(32):

x 1,n+1 = A 1 x 1,n + B 1 u 1,n + B 1 u 1,n (31)
y 1,n = C 1 x 1,n + D 1 u 1,n (32)
This brings the global state equations ( 33)-( 34) system to:

x 1,n+1 x 2,n+1 = A 1 B 1 C 2 B 2 C 1 A 2 + B 2 D 1 C 2 x 1,n x 2,n + B 1 0 u 1,n (33)
y 1,n y 2,n = C 1 D 1 C 2 0 C 2 x 1,n x 2,n (34)
As for the federation with communication sequence, the presence of the direct feedthrough requires therefore that Federate 1 receives the output from Federate 2 for the same cycle. Therefore Federate 1 has to behave according to the implementation strategy described above in Table 3.

Algebraic loop case

MATLAB/SIMULINK users are usually familiar with the concept of algebraic loops, which occur when a signal loop exists with only direct feedthrough blocks within the loop. More details on algebraic loops and how they are handled in the SIMULINK tool can be found in the corresponding user's guide [START_REF] The Mathworks | SIMULINK User's guide (r2016a)[END_REF]. To illustrate this concept, from the direct feedthrough case, let's consider a second direct feedthrough relation for Federate 2 depending on its second input u 2,n . The model of federate 2 should now be described with equations ( 35)-(36):

x 2,n+1 = A 2 x 2,n + B 2 u 2,n (35)
y 2,n = C 2 x 2,n + D 2 u 2,n (36)
Obtaining the global system equation is much more complicated. Indeed, by replacing u 1,n et u 2,n per their respective values y 2,n et y 1,n , we obtain the following output equations (37)-(38):

y 1,n = C 1 x 1,n + D 1 y 2,n (37)
y 2,n = C 2 x 2,n + D 2 y 1,n (38)
In this specific case, to calculate the output for step n, every federate needs the output from the other federates at the same instant n; this is the so-called algebraic loop concept mentionned above. This loop shows an implicit equation requiring the resolution of a system at each step of integration and the solution, if existing, can be described per equations (39)-(40):

y 1,n = (I -D 1 D 2 ) -1 C 1 x 1,n + (I -D 1 D 2 ) -1 D 1 C 2 x 2,n (39)
y 2,n = (I -D 2 D 1 ) -1 D 2 C 1 x 1,n + (I -D 2 D 1 ) -1 C 2 x 2,n (40)
where I denotes the identity matrix with the correct dimension. Therefore, in the context of a linear description, if matrices I -D 1 D 2 and I -D 2 D 1 are non singular, a solution exists and the model for the global state is given per (41)-(42):

x 1,n+1 x 2,n+1 = A 1 + B 1 (I -D 2 D 1 ) -1 D 2 C 1 B 1 (I -D 2 D 1 ) -1 C 2 B 2 (I -D 1 D 2 ) -1 C 1 A 2 + B 2 (I -D 1 D 2 ) -1 D 1 C 2 x 1,n x 2,n + B 1 0 u 1,n (41)
y 1,n y 2,n = (I -D 1 D 2 ) -1 C 1 (I -D 1 D 2 ) -1 D 1 C 2 (I -D 2 D 1 ) -1 D 2 C 1 (I -D 2 D 1 ) -1 C 2 x 1,n x 2,n (42)
However, from the HLA point of view, this problem has simply no solution since it would imply an endless exchange of messages during the same iteration n. The distributed nature of the federation prevents solving the global equations ( 41)-( 42).

An approximate solution would be that, at iteration n, one of the federates involved in the algebraic loop takes the previous value of the input, e.g., u 1,n-1 . This would actually break the loop by introducing an artificial delay. This is also one of the solutions used by SIMULINK when a global solving is not possible.

Formalization and Application

Discussion

Through this non exhaustive list of simple examples, we raised the issue of direct feedthrough within HLA federations. If this problem is not taken into account, it usually introduces delays between federate models (state and output equation calculations).

Obviously, this put into question the validity of the global simulation semantics and, depending on the size of the simulation step, this can have a significant impact on the overall system response. For example, suppose 50 federates involved in a communication sequence (as depicted in Figure 2) where each one is related with the previous by a direct feedthrough relation. If we do not care of direct feedthrough, the 50 th federate will receive its input value with a delay of 49 time units, whereas it should have received it during its first cycle. For sure, this example is a little exaggerated but it illustrates quite well the issue introduced by the presence of direct feedthrough, and why, one should apply the proper functional scheduling strategy.

To resume, there are three main cases to consider: 1) Federation does not have direct feedthrough constraint ⇒ each federate can be executed following the classical schedule (see Figure 3). 2) Federation has direct feedthrough constraint without any algebraic loop ⇒ a scheduling can be found according to the strategy described in Subsection 5.2 below. 3) Federation has direct feedthrough constraint with algebraic loops ⇒ no exact solution can be found. If possible, federate implementations must be modified in order to break these algebraic loop; then you might be in the case 1) or 2). Or, one of the federates involved in the loop should take the previous input value, as it would break the loop.

Strategy for direct feedthrough

A generic strategy is presented to schedule executions in case 2). We first define a sequence of direct feedthrough, as a sequence of output equations, beginning with an output equation without direct feedthrough and followed exclusively by output equations with direct feedthrough, Furthermore, the output value of an output equation is completely or partly the input value of the next equation.

A correct strategy to properly ensure the functional scheduling of a HLA federation implying distributed ODEs could be: 1) Identify offline all sequences of direct feedthrough; 2) Set the lookahead lk of each federate depending on the size of the sequences of direct feedthrough; 3) Decide the adequate strategy for NER and TAR usage.

SDSE example

Taking the example of the SDSE federation, the offline analysis of the ODEs describing each federate shows some direct feedthrough (bold arrows) and a possible algebraic loop. This is illustrated in figure 6 where direct feedthrough relations are illustrated using bold arrows. on Earth (latitude, longitude and altitude) to compute the corresponding atmosphere relevant parameters (temperature, pressure, air density, speed of sound) for the current time step. In addition, Flight Dynamics needs these quantities for the same time step to calculate some other outputs (but not the latitude, longitude and altitude). Strictly speaking, it is not an algebraic loop. The situation is solved by executing the output equation of Federate 5 in two steps: first, the sub vector of y that yields the aircraft position is calculated, so it can be sent to Federate 9 to obtain the atmosphere characteristics; secondly the rest of the output equation of Federate 5 is executed once the atmosphere characteristics are acquired.

Conclusion

We have presented, in this paper, a complete analysis for distributed resolution of ODEs in real-time HLA simulations. Based on our experience with the SDSE simulation test case, we offered a theoretical study on functional scheduling and deduced a set of rules in order to ensure consistent simulations. Also, we applied these rules ourselves to ensure the functional scheduling of our SDSE simulator and we compared the results of this simulation with the similar Simulink implementation to verify the correct semantic behavior. In addition, some performance metrics and tests (e.g., worst case execution time, worst case transit time, time management algorithm performances) have been conducted to guarantee the correct execution of these simulations with respect to specified real-time constraints. Interested readers may find more details and descriptions about the SDSE aircraft simulation in references [1], [2].

Figure 1 .

 1 Figure 1. Overview of the SDSE aircraft simulation

Figure 2 .

 2 Figure 2. Federation with federates communicating in series

Figure 3 2 )

 32 illustrates the temporal behavior of this federation with ∆t = 10 and lk = 1. Letters O and S describe respectively output equation calculation and state equation calculation in every federate of interest. 1) Initialization phase a) n ← 0 Simulation phase (while loop) a) Update attributes: UAV(n∆t + lk) for u 1,n b) Ask for HLA time advance to cycle n + 1 for local time (n + 1)∆t: TAR((n + 1)∆t) c) HLA time advance is granted to cycle n + 1 for local time (n + 1)∆t: TAG((n + 1)∆t) d) n ← n + 1

UAFigure 3 .

 3 Figure 3. Execution sequence for standard case

) 1 )

 1 Initialization phase a) n ← 0 b) Initial choice for x i,0 value 2) Simulation phase (while loop) a) Receive attributes y 3,n = u 2,n : RAV(n∆t + 2lk) b) HLA time advance is granted to HLA simulation time n∆t + lk: TAG(n∆t + 2lk) c) Output calculation y 3,n partir de x 3,n et u 3,n d) Publish attributes: UAV(n∆t + 2lk) de y 3,n e) Request for HLA time advance to cycle n + 1 for local time (n + 1)∆t: TAR((n + 1)∆t) f) HLA time advance is granted to cycle n + 1 for local time (n + 1)∆t: TAG((n + 1)∆t) g) Local state calculation x 3,n+1 from x 3,n and u 3,n h) n ← n + 1

Figure 5 .

 5 Figure 5. Illustration of a federation with a communication loop

Figure 6 .

 6 Figure 6. Illustration of direct feedthrough within SDSE communications

Table 1

 1 Simulation phase (while loop) a) Output calculation y i,n from state value x i,n b) Publish attributes: UAV (n∆t + lk) for y i,n c) Request for HLA time advance to cycle n + 1 for local time (n + 1)∆t: TAR((n + 1)∆t)

	. Federate 0 pseudo-code

d) Receive attributes u i,n : RAV(n∆t + lk) e) HLA time advance is granted to cycle n + 1 for local time (n + 1)∆t: TAG((n + 1)∆t) f) Local state calculation x i,n+1 from x i,n and u i,n g) n ← n + 1

Table 2 .

 2 Federate 1-2-3 pseudo-code

	RTI
	Federate 0
	RTI
	Federate 1
	RTI
	Federate 2
	RTI
	Federate 3

Table 3 .

 3 Federate 2 pseudo-code

0 value 2) Simulation phase (while loop) a) Request for HLA time advance to cycle n + 1 de temps local (n + 1)∆t: NER((n + 1)∆t) b) Receive attributes y 1,n = u 2,n : RAV(n∆t + lk) c) HLA time advance is granted to HLA simulation time n∆t + lk: TAG(n∆t + lk) d) Output calculation y 2,n partir de x 2,n et u 2,n e) Publish attributes: UAV(n∆t + 2lk) de y 2,n f) Request for HLA time advance to cycle n + 1 for local time (n + 1)∆t: TAR((n + 1)∆t) g) HLA time advance is granted to cycle n + 1 for local time (n + 1)∆t: TAG((n + 1)∆t) h) Local state calculation x 2,n+1 from x 2,n and u 2,n i) n ← n + 1

Table 4 .

 4 Federate 3 pseudo-code

	RTI
	Federate 0
	RTI
	Federate 1
	RTI
	Federate 2
	RTI
	Federate 3
	LT=0

•

  Communications between Pilots Inputs and Flight Controller federates: the output equation for Flight Controller directly depends on the pilot commands. • Communications between Sensors and Flight Controller federates: the output equation for Flight Controller directly depends on the measured values from sensors. • Communications between Actuators (i.e. control surfaces) and Flight Dynamics federates: the output equation for Flight Dynamics directly depends on the control surfaces deflections.

			Federate 9
			Environment
		Federate 4
		Actuators
	Federate 1	Federate 2	Federate 5	Federate 8
	Pilot Inputs	Flight Controller	Flight Dynamics	PFD
		Federate 3
		Engines
		Federate 6	Federate 7
		Sensors	3D Vizualisation

y 1,0 y 2,0 y 1,1 y 2,1 y 1,2 y 2,2 y 1,3 y 2,3 O y 1,0 S x 1,1 O y 2,0 S x 2,1 O y 3,0 S x 3,1 O y 1,1 S x 1,2 O y 2,1 S x 2,2 O y 3,1 S x 3,2 O y 1,2 S x 1,3 O y 2,2 S x 2,3 O y 3,2 S x 3,3 O y 1,3 S x 1,4 O y 2,3 S x 2,4O y 3,3 S x 3,4

Acknowledgments

The authors would like to thank ERIC NOULARD from ONERA for his help with CERTI and REGINE LECONTE from ISAE for her help on the evolution and the maintenance of the SDSE simulation.