Wajeb Saab

Maaz Mohiuddin

Simon Bliudze

Jean-Yves Le Boudec

Ordering Events Based on Intentionality in Cyber-Physical Systems

We consider cyber-physical systems (CPSs) comprising a central controller that might be replicated for highreliability, and one or more process agents. The controller receives measurements from process agents, causing it to compute and issue setpoints that are sent back to process agents. The implementation of these setpoints causes a change in the state of the controlled physical process, and the new state is communicated to the controllers through resulting measurements. To ensure correct operation, the process agents must implement only those setpoints that were caused by their most recent measurements. However, in the presence of replication of the controller, network or computation delays, setpoints and measurements do not necessarily succeed in causing the intended behavior. To capture the dependencies among events associated with measurements and setpoints, we introduce the intentionality relation among such events in a CPS and illustrate its differences with respect to the happened-before relation. We propose a mechanism, intentionality clocks, and the design of controllers and process agents that can be used to guarantee the strong clock-consistency condition under the intentionality relation. Moreover, we prove that our design ensures correct operation despite crash, delay, and network faults. We also demonstrate the practical application of our abstraction through an illustration with a real-world CPS for electrical vehicles.

I. INTRODUCTION A. Motivation

We consider cyber-physical systems (CPSs) comprising software agents, namely controllers and process agents (PAs), that coordinate, as shown in Fig. 1, to maintain a physical process in a desirable state. A PA is a software agent that interfaces with sensors and actuators, and controls a subprocess. The controller receives measurements from PAs, computes and issues setpoints that are implemented by the PAs through actuators. The communication network might drop, delay, reorder or retransmit messages. The controller, being susceptible to crash and delay faults, is a single point of failure and is usually replicated for high-reliability. The same general architecture of CPSs is considered in recent papers [START_REF] Mohiuddin | Axo: Masking Delay Faults in Real-Time Control Systems[END_REF]- [START_REF] Bernstein | A Composable Method for Real-Time Control of Active Distribution Networks with Explicit Power Setpoints[END_REF].

On implementation of a setpoint by a PA, the state of the sub-process is altered and the new state is communicated to the controller via a measurement. The controller uses these measurements to recreate the new state of the entire process that it then uses to compute setpoints. The setpoints computed with this state are only valid as long as the recreated state reflects the actual state of the process. Any subsequent setpoint implementations change the process state, making the former setpoints unsafe for implementation. Therefore, to achieve the desired control, the state of the process at the time of setpoint implementation must be the same as that used by the controller for computing the corresponding setpoints. Thus, before implementing a setpoint, the PA must be able to ascertain whether the setpoint reflects the state it last advertised. In other words, it must be able to infer if a received setpoint was caused by the measurement it last advertised. Similarly, a controller must be able to ascertain if a measurement received from a PA represents the most recent-state of that sub-process or the state corresponding to earlier setpointimplementations. This causal relationship can be better understood using the notion of control rounds. Software agents must be able to attribute a round number to received messages, to compare it with the round number they are currently executing, and to treat the message appropriately.

Note that, we use the term "state of the process" as a proxy for the state of the PAs. While the state of the physical process is continuous and evolving, the state of the PAs is discrete and only changes upon a setpoint implementation. In realtime control, the CPS issues setpoints at a rate faster than the dynamics of the underlying process. Thus, the evolution of the state of the process between two setpoint implementations is minimal, thereby justifying our usage of the term.

B. Need for a New Relation

When the controller is replicated, assigning a consistent round number to events requires consensus between the replicas. Due to network losses and delays, and due to software faults, consensus might require unbounded time [START_REF] Michael | Impossibility of Distributed Consensus with One Faulty Process[END_REF], making it unsuitable for real-time systems such as CPSs.

In literature, this problem is circumvented by using message labels (that represent the causal order between the messages) to infer the round number. The causal order between the messages is derived using the happened-before relation [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF]. In the presence of replication of the controller, or random network or computation delays, messages that intend to cause

M 0 M 0 SP 1 M 1 M 1 SP 0 1 Computation Delay C 1 PA C 2
Fig. 2. Illustration of why the happened-before relation is not suitable for ordering events in a CPS in the presence of replication a certain effect do not necessarily succeed, due to competing messages. This is illustrated through the example in Fig. 2. In Fig. 2, the controller is replicated, and its two replicas C 1 and C 2 receive the measurement M 0 sent by the PA. This measurement is used by each replica in the computation of a setpoint, resulting in SP 1 and SP 1 in C 1 and C 2 , respectively. In such a scenario, SP 1 and SP 1 belong to the same "generation" or the same "control round", and are said to be equivalent. SP 1 is received by the PA, and its implementation results in M 1 . However, due to a delay at C 2 , M 1 is received before SP 1 is issued. Although M 1 can be said to have happened before SP 1 , M 1 is nonetheless caused by an equivalent of SP 1 , namely SP 1 .

In the previous example, we say that SP 1 intends to have caused M 1 , but it did not succeed because a competing equivalent event (SP 1 ) was received earlier. In order to formally capture this phenomenon, we introduce the intentionality relation. This relation enables an implementation that provides a solution to the ordering problem.

C. Contributions

First, in Section III, we describe an abstraction of CPSs with one, possibly replicated, central controller, and one or more local PAs. This abstraction models the execution trace of a CPS by using two inherent relations: (1) the networking relation used to represent message exchange and (2) the computation relation used to represent the computation by software agents.

Second, in Section IV, we formally define the intentionality relation (→, read intends) by using the networking relation and the computation relation. We also give an intuition for competing events and formally define this notion under the name, intentional equivalence. Moreover, we show that physical time, on its own, is inadequate in providing strong clock-consistency under the intentionality relation [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF]. In other words, if physical time is used to timestamp and to order events a and b, where T S(a) is the timestamp of event a, then T S(a) < T S(b) ⇔ a → b.

Third, in Section V, we describe a mechanism, intentionality clocks, that uses logical clocks to attribute a round number to messages and can be used to guarantee strong clockconsistency under the intentionality relation. We formalize the desired behavior in terms of two correctness properties: safety and optimal selection. These properties specify which messages can be used in computation by software agents, and which can be discarded. We present the design of a controller and a PA that use intentionality clocks to guarantee these correctness properties. All guarantees are formally proven.

Lastly, we demonstrate a practical case study of the concepts developed throughout the paper by using, as an example, a real-world CPS for real-time control of electric vehicles (EVs) [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF]. We analyze their design in the light of the intentionality relation and show that when the controller is replicated, it violates the optimal selection property. Consequently, the CPS can enter a deadlock situation which can be avoided by using the CPS design that we proposed.

II. RELATED WORK

Here, we summarize the different bodies of work that addressed the problem of ordering events in distributed systems, and we note their shortcomings with respect to the intentionality relation. The detailed discussion of the shortcomings, along with examples, can be found in Sections I-B, IV-E and V-A. We also relate the CPS design we presented to other CPS designs found in literature.

Ordering events in distributed systems is traditionally done through a causal order, captured by the happened-before relation [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF]. Providing a causal order adheres to what is referred to as the real-time causal consistency semantic [START_REF] Mahajan | Consistency, Availability, and Convergence[END_REF]. This is achieved through one of several mechanisms such as timestamps, Lamport clocks [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF], or vector clocks [START_REF] Fidge | Timestamps in Message-Passing Systems that Preserve the Partial Ordering[END_REF].

As previously mentioned (Section I-B) causal order and the happened-before relation fail to capture the inherent "round number" in CPSs. Therefore, we introduce the intentionality relation. We elaborate, in Section IV-E, on the issues that prevent physical time and timestamping-based labeling schemes from guaranteeing the strong clock-consistency condition under the intentionality relation.

Lamport clocks and vector clocks are complementary under the happened-before relation, as Lamport clocks describe this relation, and vector clocks infer it. Intentionality clocks use a scalar clock inspired from Lamport clocks. In Section V-A, we discuss the differences between intentionality clocks and Lamport clocks in greater detail. Vector clocks were not considered as an avenue because they only provide a partial ordering between events in general distributed systems, as not all labels are comparable.

We use the same model of the CPS as previous work [START_REF] Mohiuddin | Axo: Masking Delay Faults in Real-Time Control Systems[END_REF]- [START_REF] Bernstein | Real-Time Control of an Ensemble of Heterogeneous Resources[END_REF], with a central, possibly replicated, controller and several PAs. Some previous works that use the same model [START_REF] Saab | Quarts: Quick Agreement for Real-Time Control Systems[END_REF], [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] assume the existence of a labeling scheme to achieve their goal. By providing one such labeling scheme in Section V-A, this paper complements the existing work in CPS design. We explore the CPS design in [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] in greater detail in Section VI; and we show that our system model, formalism, and mechanism apply to such a system. We also expose the possible issues in a naive extension of its design, which arise when the controller is replicated, and that can be avoided with our CPS design.

Several industrial solutions circumvent the ordering problem by using frameworks, such as the time-triggered architecture (TTA) [START_REF] Kopetz | The Time-Triggered Architecture[END_REF], which provide synchrony guarantees. However, using such frameworks requires specialized hardware, in addition to a complete redesigning of the application to fit the framework. In contrast, we propose a solution that requires neither. This facilitates deployment in existing CPSs.

III. CPS MODEL

In this section, we describe the model of the CPS we consider for defining the intentionality relation and for designing intentionality clocks. CPSs consists of four types of software agents: sensors, actuators, controllers and PAs. We consider CPSs with one central controller and one or more PAs, as shown in Fig. 1. The sensors and actuators interact with the physical process, by reading and altering its state, respectively. The software agents, namely controller and PAs, together achieve the desired control of the physical process. This model of the CPS is in agreement with the general model of CPSs considered in literature [START_REF] Mohiuddin | Axo: Masking Delay Faults in Real-Time Control Systems[END_REF]- [START_REF] Bernstein | A Composable Method for Real-Time Control of Active Distribution Networks with Explicit Power Setpoints[END_REF], and applies to a wide range of applications such as real-time control of smart grids, autonomous vehicles, and manufacturing processes.

PAs are low-level software agents responsible for controlling a sub-process, i.e., one part of the controlled process. PAs interface with the sensors and actuators as shown in Fig. 1. They implement the setpoints received from the controller through their actuators, read the state of the resources through sensors and send them as measurements to the controller. The number of PAs in a CPS is a constant and each PA is denoted with a unique identifier.

A controller performs high-level control of the physical process by receiving measurements from PAs and by sending setpoints to PAs.

We consider crash and delay faults [START_REF] Mohiuddin | Axo: Masking Delay Faults in Real-Time Control Systems[END_REF] in all the software agents, namely controller and PAs. Byzantine faults are not considered. Being a single point of failure, the controller is often replicated for high-reliability. We assume that the PAs are not replicated. Hence, a software agent is either a controller replica or a (non-replicated) PA. All replicas of the controller have different identifiers. Furthermore, we assume that the communication network between the software agents might drop, delay, reorder or duplicate messages. Byzantine contamination of messages is not considered.

We abstract the execution of a CPS as a trace of events occurring on different software agents. We define three types of events that can occur on a software agent: sending event, reception event, and timeout event. When software agent A sends a message to agent B, we say that A experiences a sending event. Upon successful reception of the message by B, we say that B experiences a reception event. Whereas, if the message is lost or is delayed beyond a deadline, B experiences a timeout event. A timeout event could also be caused by the internal logic of a software agent, such as firing of a timer (as seen in Section VI) or a response to system state.

An event is represented by the 4-tuple (sa, pa, m, l) where (1) sa is the identifier of the agent on which the event occurs, (2) pa is the identifier of either the PA on which the event occurred, or the PA for which the event is intended, (3) m is the message encapsulated in the event, and (4) l is an event label given by the software agent on which the event occurs. For sending or reception events, the encapsulated message is the measurement or setpoint exchanged, and timeout events encapsulate ⊥, representing the absence of a message.

Let C, P be the set of identifiers of all controller replicas and PAs, respectively. Then, the set of identifiers of all software agents is S = C ∪ P. Let M be the set of all messages such that ⊥∈ M, and L be a partially ordered set of labels. Then, we denote the set of all events that occur in an execution trace by E ⊂ S × P × M × L. No two distinct events have the same 4-tuple (sa, pa, m, l). Furthermore, we require that the abstract labeling scheme used to obtain L ensures that labels of events occurring on the same software agent, for the same PA, are different. In practice, such a labeling of events is achieved through physical timestamps, a permanent sequence numbering scheme, Lamport clocks [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF], Vector clocks [START_REF] Fidge | Timestamps in Message-Passing Systems that Preserve the Partial Ordering[END_REF], etc. Also, for CPSs that do not implement any labeling mechanism on events, the model still applies by successively numbering all events of each software agent with increasing integers.

Sending events are considered as output events, as they are the output of a computation at a software agent. Reception and timeout events are considered input events1 . The set of input events E i , which includes reception and timeout events, and the set of output events E o , which includes sending events, are such that E = E i ∪ E o and E i ∩ E o = ∅. Note that, due to network retransmissions, a single output event can result in different input events at the same PA, as each of these input events will have different labels l.

We consider the following computation model of a controller. In each computation, the controller uses exactly one input event from each PA and produces exactly one output event for each PA. Moreover, when the controller computes by using timeout events for one or more PAs, it is able to appropriately account for the missing information that it would have received from the corresponding reception events. Else, the controller refrains from computation of setpoints until more reception events occur (measurements are received).

To bootstrap the CPS, we assume that a controller starts with p sending events, one for each PA. These events are called initial sending events. The set of all initial sending events is represented by I.

IV. RELATIONS BETWEEN EVENTS IN A CPS

In this section, we formalize the notion of intentionality, an intrinsic relation between events in a CPS which captures the order between measurements and setpoints. First, we define the sub-relations that constitute intentionality, namely, the network relation and the computation relation. Then, in Section IV-C, we define an equivalence relation between events, called as intentional equivalence. Finally, we define the intentionality relation in Section IV-D.

For a relation r -→ and an event a, we denote by r(a) and r -1 (a) the image and pre-image of a by r -→, respectively.

A. Network Relation

Software agents exchange messages using a communication network. Thus, a network relation ( n -→) exists between events at different agents. This relation maps an output event (sending event) at one agent to an input event (reception/timeout event) at another agent. Formally, we abstract the properties of a network relation as follows.

Definition 1 (Network Relation). n -→ is a network relation, iff n -→⊂ E o × E i and • for any a ∈ E o , there exists b ∈ E i s.t. 1) a n -→ b, b.pa = a.pa, and n -1 (b) = {a} 2) If a.sa ∈ C, then b.sa = a.pa 3) If a.sa ∈ P, then b.sa ∈ C • for any b ∈ E i , there exists a ∈ E o s.t. n -1 (b) = {a} Intuitively, a n -→ b if a is a
sending event and b is its corresponding reception event or the corresponding timeout event that occurs on the intended destination. Notice that for a sending event that occurs on a controller, the corresponding input event occurs on a PA, and for a sending event that occurs on a PA, the corresponding input event occurs on a controller.

B. Computation Relation

A computation performed by a software agent can be represented as a mapping from a set of input events to a set of output events. In each computation, a controller uses p measurements, one from each PA and computes p setpoints, one for each PA. Upon reception of a setpoint, a PA implements it through the actuator, then reads the state of the sensor and sends the new state as a measurement. The set of input events used by a PA for computation is a singleton set. We abstract the properties of a computation relation ( c -→) as follows.

Definition 2 (Computation Relation). c -→ is a computation relation, iff c -→⊂ E i × E o 1) for any a ∈ E i a) If a.sa ∈ P, then ∃b ∈ E o : a.sa = b.sa, a.pa = b.pa, c(a) = {b}, and c -1 (b) = {a} b) If a.sa ∈ C, then • ∀ i ∈ P, ∃! b ∈ c(a) : b.pa = i • ∀ b ∈ c(a), b.sa = a.sa • ∀ b, b ∈ c(a), c -1 (b) = c -1 (b ) • ∀ b ∈ c(a), ∀ i ∈ P, ∃! a ∈ c -1 (b) : a .pa = i 2) for any a ∈ E o \ I, there exists b ∈ E i s.t. b.pa = a.pa and b ∈ c -1 (a).

C. Intentional Equivalence Relation

In the presence of controller replication and message retransmission, certain events that occur in the same "control round" in a CPS are functionally the same, i.e. they steer the physical process to a similar state. The intentional equivalence relation (≡) captures this.

We list the properties that will be used to define this relation. The intentional equivalence relation is defined as the smallest relation satisfying the following properties.

n -→ c, then b ≡ c 3) If a ≡ ã, a n -→ b and ã n -→ c =⇒ b ≡ c 4) If a c -→ b, a c -→ c and b.pa = c.pa, then b ≡ c 5) If a ≡ ã, a c -→ b, ã c -→ c and b.pa = c.pa =⇒ b ≡ c
The intuition behind rule (1) is that initial sending events are computed without any knowledge of the state of the system. They are polling events and do not steer the system in any direction. Thus, those sent to the same PA are equivalent.

For rule [START_REF] Saab | Quarts: Quick Agreement for Real-Time Control Systems[END_REF], the underlying intuition is that reception events and their corresponding timeout events convey similar information to the controller, as do multiple reception events corresponding to the same sending events (i.e., retransmission). A reception event informs the controller of the state of the process, whereas the corresponding timeout event forces the controller to estimate the missing state before computation. Recall from Section III that, when a controller of a CPS computes using timeout events, it accounts for the missing information, in order to ensure correctness. Thus, reception events and corresponding timeout events are equivalent.

From Definition 2, there exists a single output event resulting from a computation relation for each PA. Rule [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] states that if an event causes, after computation, two events for the same PA, then the resulting events are equivalent. In fact, the resulting events are the same event, and are therefore equivalent by the reflexive property of the intentional equivalence relation (Theorem IV.1). Finally, rules (3) and ( 5) mean that equivalent events, when subject to the same relation, result in equivalent events.

To better understand the intuition behind equivalent events resulting in similar changes to the state of a CPS, consider a controller that receives partial information from its PAs. A well-designed controller would compute only if it can reconstruct the missing information from the partially received information, thus resulting in safe setpoints. For example, a controller with two PAs, one with a 1 Hz update rate and another with a 10 Hz update rate, would require a measurement from the latter every 100 ms, whereas from the former only once every 1 s as it knows that the state of the slow PA has not changed during that second. Thus, a timeout event on the measurement from that PA is equivalent to a reception event.

Events that result from retransmissions encapsulate the same message verbatim, and are thus deemed equivalent. In the presence of replication, however, different replicas might send different sets of setpoints in a given "control round". This might be due to different sets of measurements received, a different internal state, or a non-deterministic computation. In other words, CPS equivalent events can be different and can lead to different outcomes. Therefore, in practice, CPSs must be able to live with this. This can be done, for instance, with an agreement on input as shown in our companion work, Quarts [START_REF] Saab | Quarts: Quick Agreement for Real-Time Control Systems[END_REF]. Alternatively, this is achieved if the setpoints are idempotent and the CPS logic permits such deviations [START_REF] Panda | SCL: Simplifying Distributed SDN Control Planes[END_REF].

The following theorem states that the intentional equivalence is, indeed, an equivalence relation.

D. Intentionality Relation

We define the intentionality relation (→), where a → b is read as "a intends b", by using the relations defined in the previous sections. It is the smallest relation satisfying the following properties: The proofs of Theorems IV.1, IV.2, and IV.3 are not central to the discussion, and are thus not included for spatial consideration. We discuss them further in the technical report [START_REF] Saab | Ordering Events Based on Intentionality in Cyber-Physical Systems[END_REF].

1) If a n -→ b, then a → b 2) If a n -→ b and ã ≡ a, then ã → b 3) If a c -→ b, then a → b 4) If a c -
Let E p = {e ∈ E|e.pa = p}. Recall that E p /≡ represents the factorization of the set E p by the relation ≡. Then, as a consequence of Theorem IV.3, the intentionality relation induces a total order on E p /≡, for any p ∈ P. In other words, the intentionality relation induces a total order on any set of events concerning the same PA and belonging to different equivalence classes.

In Section V, we use the definition of intentionality to formally specify the desirable correctness properties of a CPS. We also present a CPS design and prove that it guarantees the said properties.

E. Intentionality and Physical Time

In this section, we answer the question "Is physical time sufficient to guarantee the strong clock-consistency condition under the intentionality relation?"

As CPSs are real-time systems, they generally keep track of physical time. To maintain synchronized global time on all software agents, their physical clocks are synchronized either using GPS-based clock-synchronization or network-based

a 2 P A 1 C P A 2 a 1 b 2 c 2 d2 h 1 d 1 e 1 g 1 f 1 c 1 b 1 d 2 e 2 f 2
Fig. 4. Difference between temporal order and intentionality due to delays clock-synchronization (e.g., PTP [START_REF]IEEE Standard for a Precision Clock Synchronization Protocol for Networked Measurement and Control Systems . IEEE Std 1588-2002[END_REF], NTP [START_REF] David | Internet Time Synchronization: the Network Time Protocol[END_REF]). These timesynchronization solutions provide a synchronization accuracy δ that ranges from sub-microsecond to one millisecond. CPSs often leverage the availability of synchronized physical clocks to reason about the temporal ordering of events. However, as we will see in the following examples, the temporal ordering of events does not coincide with intentionality.

Consider a CPS with two controller replicas and a single PA, as shown in Fig. 3. We will consider a perfect timesynchronization (δ = 0). Here, a → b, a ≡ ã =⇒ ã → b. However, the time of occurrence of the effect b, is less than that of the "cause" ã.

Another example emphasizing the difference between temporal order and intentionality, shown in Fig. 4, concerns a CPS with one non-replicated controller and two non-replicated PAs. Due to network delay, the reception event b 2 occurs on PA 2 much later than the reception event b 1 at PA 1 . As a result, the reception event d 1 from PA 1 occurs much earlier than the corresponding reception event d2 from P A 2 . Instead, the controller moves on with a timeout event d 2 for P A 2 . Consequently, the events e 1 , f 1 , g 1 and h 1 take place. Then, we have

a 2 → d 2 → e 1 → h 1 . But, d 2 ≡ d2 =⇒ d2 → h 1 .
However, on the controller, the time of occurrence of h 1 occurred is less than that of d2 . Therefore, as the correct temporal order is h 1 before d2 , it does not coincide with intentionality.

Hence, we conclude that on their own, physical clocks are not sufficient to reason about intentionality and require an additional mechanism to do so. In Section V-A, we present a mechanism that describes the intentionality relation by using logical clocks instead of physical clocks. Note that, such a mechanism does not require synchronized physical time.

V. CPS DESIGN

As discussed in Section I, the controller and PAs must ensure that the events used in their computation reflect the most recent state of the process. The state of the process changes with each setpoint implementation, i.e., a computation event by a PA. Thus, when a controller uses an input event e : (sa 1 , pa 1 , m 1 , l 1 ) and there exists another event e : (sa 1 , pa 1 , m 2 , l 2 ) such that e → e , then the message m 2 reflects a more recent state of the sub-process controlled by PA 1 than the message m 1 . If the controller uses the "old" event e in the computation of setpoints, then the resulting setpoints might not be compatible with the most-recent state of the process, thereby causing incorrect control of the CPS.

Definition 3 formally specifies the safety property that the software agents in a CPS must satisfy in order to guarantee correct control. Definition 3 (Safety). If a software agent uses an event a for computation, then the last sending event b that occurred on this software agent, where b.pa = a.pa, is such that b → a.

The safety property requires that (1) a PA must use only those events in computation that have accounted for its most recent state, captured by its last measurement sending event, and (2) a controller must use only those events in computation that reflect the state change caused by its last setpoint sending events.

Notice that discarding all messages can trivially satisfy the safety property. However, this will render the physical process uncontrolled by a complete loss of availability of the control. Therefore, the selection of events must be optimal, i.e., only those events that violate the safety property must be discarded. This property is formally defined as follows.

Definition 4 (Optimal Selection

). An event a must only be discarded if there exists another event b, such that b → a that: (1) occurred on this software agent, or (2) this software agent was informed that b occurred on its replica.

The first part of the optimal selection property, on its own, states that an event a must be accepted by a software agent if the last sending event b on this software agent intended to cause a. Recall that if b → a, then all events that occurred on this software agent before b also intend to cause a. Thus, a presents new information, i.e., information about the state of the process after the implementation of the last-sent setpoint. In the presence of replication, however, an event that was intended by the last sending event, does not necessarily present new information, as seen in the following example.

A controller replica that last computed in round l might receive two measurements from the same PA: one from round l+1 and another from round l+2, before it is ready to compute. This can occur due to another controller replica driving the CPS into round l+2. In this scenario, the first controller might ignore the message from round l+1 as the message from round l + 2 supersedes it. The controller will, therefore, compute using reception events from round l + 2, and declare timeout events for PAs from which it has not received measurements from that round. This condition is captured by the second part of the optimal selection property.

Next, we describe the design of a CPS in which the software agents satisfy the safety and optimal selection properties. Our design uses a label scheme, intentionality clocks, that is adapted from Lamport clocks to guarantee strong clockconsistency under the intentionality relation. We describe the design of intentionality clocks in Section V-A; and the design of controllers and PAs that use the intentionality clocks in Sections V-B and V-C, respectively. We present the formal guarantees in Section V-D.

A. Intentionality clocks

Intentionality clocks is a mechanism to maintain, update and synchronize logical clocks across all software agents in a CPS. It is adapted from the Lamport clocks abstraction that was designed for general distributed systems, to accommodate the specificities of events in CPSs.

Each agent maintains and updates a local logical clock that is used to set the event label of sending events on that agent. These labels are communicated along with the message encapsulated in the event and are used to obtain the event labels of the corresponding reception and timeout events.

The labels are assigned such that they guarantee the strong clock-consistency condition [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF]. In other words, for two events a and b, we say that a clock mechanism describes the intentionality relation if the event labels are such that a.l < b.l =⇒ a → b and a.l = b.l =⇒ a ≡ b. We say that the clock mechanism infers the intentionality relation if a → b =⇒ a.l < b.l and a ≡ b =⇒ a.l = b.l. The strong clock-consistency condition is satisfied if a clock mechanism both describes and infers the intentionality relation.

The design of intentionality clocks is given by the following rules (for details, see Sections V-B and V-C):

1) The event label of a sending event at an agent is the value of its logical clock, right before the sending event. 2) If a is a sending event and b is its corresponding reception or timeout event, then b.l = a.l + 1.

3) The logical clock of an agent is only incremented before a computation. It is never decremented. 4) The logical clock of delayed software agents is resynchronized using the labels of the received events.

There are two main distinctions in the design of intentionality clocks when compared to Lamport clocks [START_REF] Lamport | Time Clocks and the Ordering of Events in a Distributed System[END_REF]. First, in our solution, the logical clock at an agent is incremented only when the agent performs a computation. This enables the controller to infer the intentionality relation by using its local logical clock, and to have a notion of the "control round". In contrast, the Lamport clock at an agent is updated after every event that occurs at that software agent Thus, the agents lose the information required to infer the intentionality relation and to have a notion of the "control round".

Second, in intentionality clocks, the label of a reception event is one more than the label of the corresponding sending event. In contrast, in Lamport clocks, for a reception event b that occurs when the value of the logical clock of the agent is C is b.l = max(a.l, C) + 1, where a is the corresponding sending event. Due to the presence of delayed controller replicas or message retransmissions by the network, reception events from previous "control rounds" might have a higher label than reception events from current "control round". Consequently, it cannot describe or infer the intentionality relation, i.e., cannot guarantee the strong clock-consistency condition.

The formal guarantees of intentionality clocks are presented in Section V-D and proven in the Appendix.

B. Controller Design

Algorithms 1 and 2 describe the design of a controller with intentionality clocks; this design satisfies the safety and optimal selection properties. The model of the controller is the same as that used in our previous work [START_REF] Mohiuddin | Axo: Masking Delay Faults in Real-Time Control Systems[END_REF], [START_REF] Saab | Quarts: Quick Agreement for Real-Time Control Systems[END_REF]. The parts in red are our modifications for satisfying the aforementioned properties (together with the implementation of Algorithm 3).

Each controller maintains a logical clock C, a list of input events (i.e., reception and timeout events) Z, and a list of their corresponding event labels L.

Upon receiving a measurement from a PA, the controller declares a reception event with a label one more than the label of the received message (Algorithm 1, line 9). The controller adds the measurement to Z, and adds the label of the reception event to L.

The controller also occasionally checks if it has accumulated enough information about the state of the physical process, through the measurements, required to compute setpoints. To this end, it uses the ready_to_compute() function. We make no assumptions on this function or the frequency with which the controller invokes it. The controller can choose to start a computation of setpoints by considering any form of information provided by measurements from Z. Moreover, the labels in L provided by intentionality clocks expose additional information to the ready_to_compute() function, by giving insight into the intentionality relation between the events.

Then, the controller computes setpoints by calling the compute() function that uses Z, L and C as shown in Algorithm 2. In this function, we specify how the controller must use the logical clock and the labels in order to satisfy the safety and optimal selection properties. First, the controller computes the highest label of the events it has seen: C . This represents the most recent events the controller has encountered. Then, for all PAs from which the highest label of received events is less than C , the controller declares a timeout event (line 4), thereby explicitly acknowledging that it lacks the most recent information from this PA. The controller can then account for this missing information in the subsequent computations. The logical clock C is set as C + 1 to mark the computation operation and the resulting setpoints are issued by way of sending events with the label being the current logical clock.

Note that, in Algorithm 2, the computation of setpoints takes as input a set of input events comprising exactly one reception or timeout event from each PA. This is in accordance with the computation relation ( c -→) described in Section IV-B. The controller is designed to be soft state [START_REF] Lui | On the Robustness of Soft State Protocols[END_REF]. When a controller boots or reboots after a crash, its logical clock is set to zero, and the lists Z and L are reinitialized. Thus, it would use all subsequent reception or timeout events for computation, because their labels would be ≥ 0. This behavior is in accordance with the safety property, as a freshly rebooted controller de-facto has no last setpoint sending event. However, as described in Section V-C, these sending events would be disregarded by the PAs as they do not reflect their most recent state. Note that, upon booting or rebooting, the controller sends setpoints corresponding to the initial sending events, indicated by S 0 , with the label 0. From lines 1 and 9 in Algorithm 2, we see that the logical clock of a newly booted controller is re-synchronized with that of the other software agents before computation, by taking the maximum of the labels of the received measurements. This is discussed further, in Section V-C.

C. PA Design

Algorithm 3 describes the design of a PA with intentionality clocks; the design complements Algorithms 1 and 2 in satisfying the safety and optimal selection properties. Each PA maintains a clock C that is initialized with 0 upon booting.

Upon reception of a setpoint from a controller, the PA declares a reception event with a label one more than that received in the message. Then, the PA compares the label of the event with its local logical clock. If the reception event has a higher label, then the PA implements the setpoint, else it is discarded because it violates the safety property. In other words, a reception event with a label less than the logical clock of the PA means that the corresponding sending event was not computed with the most recent state of this PA. In this way, setpoints from delay-faulty controllers or freshly booted controllers are not implemented, thereby upholding the safety property.

After implementing the setpoint through an actuator, the PA increments its logical clock to mark the computation of a new measurement. The PA computes a new measurement through a sensor and sends the measurement to the controller by a sending event labeled with the current logical clock.

Each PA stores its logical clock before computing measurements. When a PA recovers after a crash, it initializes its logical clock to the last stored value 2 . In this way, the PA keeps track of the last state it advertised to the controller.

Controllers and PAs update their local logical clocks to reflect the labels they observe (Algorithm 2 line 1, Algorithm 3 line 10). Notice that, upon receiving a setpoint with a label lower than its logical clock, a PA also sends both the latest computed measurement and the current value of the local clock. This serves to re-synchronize the software agents that miss some control rounds due to messages losses, crashes and recoveries, or delays.

D. Formal Guarantees

We formally prove that our mechanism of intentionality clocks and Algorithms 1, 2, and 3 guarantees safety and optimal selection. The first step lies in proving that intentionality clocks under our mechanism infer and describe the intentionality relation. Proof. The proof can be found in Appendix B.

VI. CASE STUDY: CPS FOR SCHEDULING ELECTRIC VEHICLE CHARGING

We present, via a case study, a practical application of the intentionality relation and the intentionality clocks mechanism. We take the example of a CPS for scheduling the charging of a fleet of electric vehicles (EVs), which provides a schedule that accounts for both the vehicles' demand and the vehicle-to-grid regulation services [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF]. One of the goals of the CPS is to charge each vehicle based on the vehicle's demand, and the other goal is to provide frequency support for the grid. Modulating the charging schedule of the EVs, by charging at a higher or lower rate, or even discharging into the grid for some time in cases of downward excursion of frequency, can provide frequency support. In practice, the EVs must respond to regulation requests every two to four seconds [START_REF] Kamboj | Deploying Power Grid-Integrated Electric Vehicles as a Multi-Agent System[END_REF].

The paper [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] presents a solution to the problem assuming an ideal communication, and without considering the failure of the software agents. However, such a mission-critical control requires high levels of reliability in a real deployment. It is, therefore, desirable to replicate the controller.

In Section VI-B, we use the intentionality relation to analyze their CPS design for possible issues that could arise when the controller is replicated. We find that a naive extension to their design violates the optimal selection property. Consequently, due to software and network faults, the CPS can encounter a deadlock situation, whereby frequency support will no longer be achieved. This might result in the instability of the underlying electrical grid.

A. CPS Model

In this section, we analyze the system model of [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] and show how our system model, presented in Section III, applies.

Figure 5 presents an architectural view of the system consisting of EVs and the aggregator, as shown in [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF]. It comprises EVs that represent the controlled sub-process from our model, the EV agents labeled 1 through n that correspond to the PAs from our model and the central aggregator that is the controller of the CPS. The controller receives a charging/discharging (henceforth referred to as charging) schedule from each PA, which represents the measurements in our model. The controller computes control signals (setpoints) for each of the EVs, using the received charging schedules.

The algorithms of the controller and PA used for achieving the desired vehicle-to-grid regulatory service are described in [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] in Algorithms 2A and 2B. We abstract this process and summarize it, as shown in Algorithm 4. The controller periodically starts an iterative process that consists of several control rounds, subject to convergence of the algorithm. Each iterative process is independent from the previous one, as if the controller had a fresh start. The triggering of a new iterative process (line 2) is a timeout event in our model, which further causes n sending events, one for each PA (line 4). These represent the initial sending events as they are computed without using measurements.

From Algorithm 4, we see that the CPS in [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] uses a labeling scheme with label m. In Section VI-B, we show that this labeling scheme violates optimal selection (Definition 4), consequently the system can enter a deadlock situation. In line 6, we see that the controller waits for schedules from each PA before beginning computation in line 7. Thus, the ready_to_compute() function of this CPS is the presence of one reception event from each PA, corresponding to the current round with label m. Timeout events do not occur within an iteration.

In lines 7-13, we see the computation relation takes as input one schedule from each PA and produces one control signal for each PA. In other words, the computation relation takes as input one reception event from each PA and outputs one sending event for each PA. This is same as our computation relation (Definition 2).

The iterative process (lines 5-13) continues until the control has converged. It terminates with the sending of control and stop signals in the setpoints of the last computation.

Each EV agent (PA) also keeps track of the on-going control round by using the indicator m. It only accepts charging schedules (setpoints) that belong to the current round (line 20). Upon receiving a setpoint from the current round, the PA checks if it is accompanied with a stop signal (line 22). The presence of a stop signal is an indication of the termination of the iterative process and results in implementation of the setpoint. Alternatively, when the stop signal is absent, the PA sends its new charging schedule as a measurement to the controller and increments m by one.

B. Deadlock due to Violation of Optimal Selection

Here, we describe a scenario in which the CPS in [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] can enter a deadlock situation when the controller is replicated, due to the controller replicas having different round indicators.

Consider a scenario with two replicas of the controller C 1 and C 2 , with indicator m 1 and m 2 , respectively. Consider a PA, PA 0 with indicator m 0 . Consider a situation when both controllers sent out setpoints to PA 0 with label 5. Then, m 1 = m 2 = 5. PA 0 receives this setpoint, implements it, sends a measurement with label m 0 = 5 and increments its indicator to 6. Thus, after this round, we have m 1 = m 2 = 5 and m 0 = 6. Now, if C 2 does not receive the measurement from PA 0 and C 1 received measurements from all PAs, C 1 will begin computation and C 2 will be stalled. After this computation, we have m 1 = 6 and m 2 = 5. Moreover, PA 0 implements the setpoint and increments its indicator m 0 to 7.

Next, let C 1 crash due to a software failure and reboot. Then, it requests for schedules with m 1 = 0. However, as P A 0 is expecting messages from round 7, it ignores these requests. P A 0 also discards any messages from C 2 , as they will have a label 5. Moreover, C 2 discards the received measurements with label 6 because its label is m 2 = 5. This is a violation of Algorithm 4: Abstraction of the iterative process of scheduling EV charging [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] 1 At the Aggregator (Controller) 2 repeat periodically the optimal selection property by C 2 . Hence, the CPS enters a deadlock state because there is no mechanism to resynchronize the counters.

Our design (Section V) satisfies the optimal selection property even in the presence of controller replication, and software and network faults (as shown in Theorem V.2). Therefore, the problem encountered by the design in [START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF] can be avoided by applying intentionality clocks and tuning the design of the software agents according to Algorithms 1, 2, and 3.

VII. CONCLUSION AND FUTURE WORK

We address the problem of enabling software agents in a CPS, namely controller and PAs, to provide a notion of rounds of computation in presence of network losses or delays, or replication of the controller. We show, that in such settings, the causal (or happened-before) relation, traditionally used in distributed systems literature, does not enable capturing the control rounds. Instead, we introduce a new relation that we call "intentionality relation". We formally define this relation for CPSs with one, possibly replicated, central controller and one or more PAs. A possible avenue for future work is to extend the intentionality relation to a wider range of CPS, specifically, comprising a hierarchy of controllers, or spontaneous sensors that send out-of-band measurements to the controllers.

We present a clock mechanism, intentionality clocks, that can be used to both describe and infer the intentionality relation. We also formalize the correctness properties, namely safety and optimal selection, that describe how the agents must treat events in order to respect intentionality. We present the design of a controller and a PA that guarantees the correctness properties. Lastly, through a case study of a real-world CPS for charging EVs, we demonstrate the practical relevance of the introduced concepts.

In future work, we intend to extend our design to more generic CPSs as mentioned above. Moreover, we are currently in the process of implementing our controller and PA design in a CPS for real-time control of electric grids [START_REF] Bernstein | A Composable Method for Real-Time Control of Active Distribution Networks with Explicit Power Setpoints[END_REF]. We will study the impact of violations of the correctness properties (safety and optimal selection) on the physical process, through experiments in a virtual commissioning environment [START_REF] Achara | T-RECS: A Software Testbed for Multi-Agent Real-Time Control of Electric Grids[END_REF]. This enterprise is a preparation for the deployment of the said CPS for real-time control of a medium-voltage grid.

VIII. ACKNOWLEDGMENTS

This research was supported by the "SNSF -NRP 70" Energy Turnaround project.

1 )

 1 If a and b are the initial sending events at controller replicas, and a.pa = b.pa then a ≡ b 2) If a n -→ b and a

Fig. 3 .

 3 Fig. 3. Difference between temporal order and intentionality due to replication of the controller Theorem IV.1. "Intentional equivalence" is reflexive, symmetric and transitive.

  → b and ã ≡ a, then ã → b 5) If a → b and b → c, then a → c Based on these properties, the following theorems hold. Theorem IV.2. For any two events a, b: a → b =⇒ b → a. Theorem IV.3. For any two events a, b, such that a.pa = b.pa: (a → b ∧ b → a) ⇐⇒ a ≡ b.

Algorithm 1 : 2 C 3 Z 4 S 5 L 9 15 if decision then 16 SAlgorithm 2 : 9 C

 123459151629 Abstract model of a controller 1 on boot or reboot ← 0; // CPS clock on this controller ← {}; // List of input events ← S 0 ; // List of setpoints to be sent ← {}; // List of labels of input events 6 Issue S with label C 7 end; 8 on reception of a message m with label l from a PA i Declare reception event a : (sa, i, m, l + 1); 10 Add m to Z; 11 Add a.l to L; 12 end; 13 repeat 14 decision ← ready to compute(C, Z, L); , C, Z, L ← compute(C, Z, L); 17 Issue S with label C; // Sending events Function: compute(C, Z, L) 1 C ← max(C + 3, max(L)); 2 for each PA i do 3 if the maximum label in L from PA i is not equal to C then 4 Declare timeout event a : (sa, i, ⊥, C ); 5 Add ⊥ to Z; 6 Add a.l to L; ← C + 1; 10 S ← setpoints computed using measurements with label C ; 11 Return S, C, Z, L;

Algorithm 3 : 5 C

 35 Model of PA with causal clocks. ←stored C; 6 end; 7 on reception of a message m with a label l from a controller 8 Declare reception event a : (pa, pa, m, l + 1); the controller with label C; 17 end;

Theorem V. 1 (

 1 Strong Clock-Consistency). In a CPS that implements Algorithms 1, 2, 3: for any two events a and b, C (a) < C (b) ⇐⇒ a → b C (a) = C (b) and a.pa = b.pa ⇐⇒ a ≡ b Proof. The proof can be found in Appendix A. Theorem V.2 (Safety & Optimal Selection). A CPS that implements Algorithms 1, 2 and 3 guarantees safety and optimal selection.

Fig. 5 .

 5 Fig.5. Architecture and information flow of a CPS for EV charging[START_REF] Lin | Optimal Scheduling with Vehicle-to-Grid Regulation Service[END_REF]

3 m ← 0; 4 5 repeat 6 if 7 Perform computation of control signals; 8 m ← m + 1 ; 9 if 14 Send 16 17

 345678191416 Send a request for schedules with label m to each PA; schedules labeled m from each PA are received then control has not converged then 10 Send new control signals with label m to each PA; control signals and stop signals with label m to each PA; 15 forever; At the EV Agent (PA) 18 m ← 0; 19 on reception of a control signal with label k from aggregator 20 if k == m then 21 Compute new charging schedule; 22 if stop signal received then 23 m ← 0; 24 Implement charging schedule until next control signal; 25 else 26 Send schedule to aggregator with label m; 27 m ← m + 1;

This dichotomy of input-output events is similar to that of sendingreceiving events used in classic distributed systems literature. We use a different name because we also have timeout events in our model.

Constantly increasing counters might cause a counter overflow. However, a 64-bit counter incremented once every millisecond takes much longer than the lifetime of any CPS to wrap-around.

Inductive hypothesis: C (c) = C (c). (For cases 2 and 4, c = c and this hypothesis holds trivially).

Inductive step: We show that the statement of the lemma also holds for a and b. In cases 2 and 3, by Lemma A.6:

In cases 4 and 5, by Lemma A.5: Thus, e.pa = f.pa. By Lemma A.6:

Thus, by the inductive hypothesis: e ≡ f . Thus, by rule 3 of intentional equivalence: a ≡ b.

Case 2: a and b are output events. By Definition 2: ∃ g, h: g c -→ a, h c -→ b, and g.pa = h.pa. Then, from Lemma A.5:

Thus, by inductive hypothesis: g ≡ h. Thus, by rule 5 of intentional equivalence: a ≡ b.

Case 3: a is an input event and b is an output event. We prove that this is an impossible case. By Definition 1: there exists an input event e, such that e.pa = a.pa and e n -→ a. By Definition 2: there exists an output event f , such that f.pa = b.pa and f c -→ b. Then, by Lemma A.6:

Then, C (e) = C (f ) and e.pa = f.pa. From the induction hypothesis: e ≡ f . By the properties of the intentional equivalence relation, e and f are either both input events or both output events. Contradiction.

We prove the lemma by induction on k. Inductive hypothesis: Let, for some k > 1 ∀ e, f :

Base case: For

From the inductive hypothesis: a → c. Also, from the base case: c → b. Therefore, from rule 5 of intentionality: a → b.

Proof. In Algorithm 2, a is either from a reception event from a PA with label C or a timeout event with a label C (line 6).

In Algorithm 2 line 4: we declare timeout events with label C , for each PA i, such that the maximum label in L of the events corresponding to i is different from C . We trace the events that caused one such timeout event b. Since a n -→ b, it follows that a is the sending event at a PA, such that a.pa = b.pa, which was lost or delayed thus causing the timeout event b. Let g be the last sending event at a controller, such that there exists a chain of events g n

, where c is a reception event at the controller on which b occurred, and c.l = C . (Definitions 1 and 2). Thus, from result of this lemma at the PA: d.l = C -1. From Lemma A.5: f.l = C -2. From the earlier statement at the controller: g.l = C -3. The controller declares timeout events to acknowledge that it lacks information from some PAs, that it has from others. This information is the state of the sub-processes after the implementation of the setpoints encapsulated in the last sending events that resulted from the same computation relation at some controller replica. Thus, from Definitions 1 and 2: there exists another chain of events, g n -→ h c -→ a, such that h.pa = a.pa. From result of this lemma at the PA: h.l = C -2. From Lemma A.5: a.l = C -1.

APPENDIX B PROOF OF THEOREM V.2

A CPS that implements Algorithms 1, 2, 3 guarantees safety and optimal selection.

Proof. The proof has two parts: safety and optimal selection.

Safety

At the controller: Algorithm 2 line 10, the controller computes setpoints using measurements with label C . C is greater than the label of the last setpoint issued at line 17 in Algorithm 1 with label C, as C = max(C +3, max(L)) (Algorithm 2, line 1). From Theorem V.1:

Thus, when an event a with label C is used by the controller for computation, and another event b was the last event issued by the controller with label C, then C > C =⇒ b → a.

At the PA: The computation at PA, i.e., implementation of a setpoint followed by subsequent computation of a measurement, is triggered only if Algorithm 3 line 9 is true. Therefore, if an event a with label a.l is used in a computation by a PA when its CPS causal clock is C, then C < a.l. However, C is the label of the sending event b, corresponding to the last measurement issued. Therefore, by Theorem V.1:

Optimal Selection

At the controller: Consider a reception event a with label a.l at a controller, when its clock is C.

Case 1: a.l ≤ C. Let b is a sending event corresponding to the last setpoint sent to the same PA. From Theorem V.1:

From Algorithm 2 line 1, we have a.l < C + 3 < C . In line 10, the controller computes with events of label C . Thus, the event a is not used in computation, i.e., discarded. Hence, a must be discarded. From Algorithm 2 line 1, we have a.l < max(L) < C . In line 10, the controller computes with events of label C . Thus, the event a is not used in computation, i.e., discarded.

Case 3: a.l > C and a.l = max(L) From Lemma B.1: a.l > C + 3. From Algorithm 2 line 1: a.l = C . In line 10, the controller computes with events of label C . Thus, the event a is used in computation, i.e., not discarded.

At the PA: An event a is discarded by the PA only if a.l ≤ C (Algorithm 3, line 9), where C is the label of the event b corresponding to the last measurement sent. Thus, from Theorem V.1, if a is discarded, b → a. Proof. We start by proving the statement for output events at the controller, by using induction. , where e and a are both output events. Thus, the statement is true for the first two sending events on all controllers.

Inductive hypothesis: Consider an instance of the execution trace at time t 0 at which there is no event in the CPS b, such that C (b) ≥ l. At this instant, let A be the set of all sending events on all controller. Clearly, ∀a ∈ A, C (a) < l. We assume the hypothesis be true for all a ∈ A.

Inductive step:

Consider a 0 ∈ A, be such that ∀ a ∈ A, C (a 0 ) ≥ C (a). By the induction hypothesis: C (a 0 ) = 4k, where k ∈ N. We show that the hypothesis is also true for the first sending event b that occurs at a controller after t 0 such that a 0 → b. f is a sending event on a controller and C (f ) > C (a 0 ). Thus, f / ∈ A. Also, f is the first sending event that occurred after t 0 . Thus, the hypothesis holds for all sending events at controllers. Sending events at controllers have label of the form 4k. By Lemma A.6: input events at PAs have label of the form 4k + 1. From Algorithm 3: output events at the PAs have label of the form 4k + 2. By Lemma A.6: input events at the controller have label of the form 4k + 3.


Anne Bouillard

Céline Comte

Élie De Panafieu

Fabien Mathieu

Of Kernels and Queues: when network calculus meets analytic combinatorics NetCal 2018 (

Stochastic network calculus is a tool for computing error bounds on the performance of queueing systems. However, deriving accurate bounds for networks consisting of several queues or subject to non-independent traffic inputs is challenging. In this paper, we investigate the relevance of the tools from analytic combinatorics, especially the kernel method, to tackle this problem. Applying the kernel method allows us to compute the generating functions of the queue state distributions in the stationary regime of the network. As a consequence, error bounds with an arbitrary precision can be computed. In this preliminary work, we focus on simple examples which are representative of the difficulties that the kernel method allows us to overcome.

Introduction

The development of new wireless communication technologies (5G) shed a new light on queueing theory, as the strong requirements on buffer occupancy, latencies, and reliability, bring the need for accurate dimensioning rules. In many scenarios, data packets arrive by batches and are processed by a server that can deal with a fixed number of packets per time slot [START_REF] Sehier | Transport Network Design for FrontHaul[END_REF]. The G/D/1 queue is thus a natural model.

A powerful tool to analyze such queues is Stochastic Network Calculus (SNC) [START_REF] Fidler | A guide to the stochastic network calculus[END_REF]. The aim of SNC is to derive precise error bounds on the performance of systems, combining deterministic network calculus and probabilistic tools.

Among the techniques developed so far, the Tailbounded approach [START_REF] Jiang | Stochastic Network Calculus[END_REF] introduces a violation probability in the parameters of the deterministic setting. It makes possible the computation of error bounds in networks, like in [START_REF] Ciucu | Scaling properties of statistical end-toend bounds in the network calculus[END_REF], but these bounds are usually loose. A second technique, introduced in [START_REF] Chang | Performance Guarantees in Communication Networks[END_REF], relies on moment generating functions. It can be very accurate for one queue. For example, in [START_REF] Poloczek | Scheduling analysis with martingales[END_REF][START_REF] Ciucu | On multiplexing flows: Does it hurt or not?[END_REF], the authors obtain tight upper and lower bounds for the single-server case under various service policies and arrival processes, using martingales and Doob's inequality. However, for more general topologies, the method becomes non applicable due to interdependencies between the processes. Recently, some (looser) bounds have been computed using Hölder's inequality [START_REF] Beck | Advances in theory and applicability of stochastic network calculus[END_REF][START_REF] Nikolaus | On per-flow delay bounds in tandem queues under (in)dependent arrivals[END_REF].

The use of generating functions to investigate random processes is the core principle of analytic combinatorics, a subfield of combinatorics (see [START_REF] Flajolet | Analytic Combinatorics[END_REF]). This community developed mathematical tools to study random walks [START_REF] Flajolet | The formal theory of birth-and-death processes, lattice path combinatorics and continued fractions[END_REF], such as the kernel method [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF][START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF], described later. The link between random walks and queueing theory is known and results on the former were transferred to the latter [START_REF] Fayolle | Random Walks in the Quarter Plane: Algebraic Methods, Boundary Value Problems, Applications to Queueing Systems and Analytic Combinatorics[END_REF].

In this article, we show how generating functions and the kernel method can be applied to derive precise results on queueing systems.

In §2, we first recall the main definitions and notations of generating functions. The main contribution of the paper is given in §3, where we show in detail how to apply the kernel method to study the GI/D/1 queue. Although the result itself is well-known (we retrieve the Pollaczek-Khinchine formula), the interest of the analysis is that it contains all the pieces for further extensions, such as several flows of packets, several queues, or non i.i.d. arrivals. Some of these extensions are developed in §4: random service, multi-flow and multi-queue. Finally, we confront our results with simulations in §5.

Generating functions

In this section, we recall some basics of generating functions. Let (a n ) n≥0 be a sequence of non-negative numbers. Its generating function is the formal series

A(u) = n≥0 a n u n .
The n-th monomial a n will also be denoted by [u n ]A(u). In combinatorics, a n is often the number of objects of size n within a given family. In probability, a n is usually the probability that a random variable A with values in N is equal to n:

A(u) = n≥0 P(A = n)u n .
In that case, we write A ∼ A; the convergence radius ρ of the function A is at least 1,

A(1) = 1, A ′ (1) = E[A],
and we assume lim u→ρ A(u) = +∞ to simplify the asymptotic analyses.

Two elementary operations can be performed on generating functions. Suppose that A and B are the generating functions of two random variables A and B, respectively.

1. If the events {A = n} and {B = n} are disjoint for each n ∈ N, then

A(u) + B(u) = n≥0 P({A = n} ∪ {B = n})u n .
2. If A and B are independent, then A(u)B(u) is the generating function of the r.v.

A + B.

Consider the example of a Galton-Watson tree, which is a branching process where the number of children of each node is i.i.d. with distribution given by the generating function A. The number of nodes of the tree is

X = 1 + A k=1 X k ,
where A ∼ A is the number of children of the root and X k is the number of nodes in the subtree rooted at the k-th child of the root. X k has the same distribution as X, hence the same generating function, denoted by T A . Therefore, we obtain

T A (u) = uA(T A (u)).
(

) 1
This equation characterizes T A . Figure 1 shows how T A (u) is computed. Being the generating function of a probability distribution, T A (u) must be a solution of Eq. ( 1) that is analytic at 0, also known as a small root of the equation. T A (u) is then the abscissa coordinate of the first intersection of A(x) with the line x/u. There is a maximal value ρ T A > 1 of u for which a root exists.

x x x/u A(x) β 1 T A (u) x/ρ T A T A (ρ T A ) Figure 1: Equation T A (u) = uA(T A (u)).
By deriving both sides at u = 1, one gets that E

[X] = 1 + E[X]E[A], so that E[X] = (1 -E[A]) -1 if E[A] < 1.
Adopting a combinatorics viewpoint allows us to consider generating functions that do not represent a distribution. For example, T k A is the generating function of the distribution of the total size of k independent Galton-Watson trees, and

1 1-T A = k≥0 T k A
is the sum of the distributions for all possible k. Let (X j ) j≥0 denote the sequence of the sizes of i.i.d. Galton-Watson trees, then

[u n ] 1 1 -T A (u) = P(∃k, X 1 + • • • + X k = n).
Studying the behavior of this series will prove useful to derive the asymptotic probability that an arbitrary number of trees has a given total size. As T A (u) < 1 for all 0 ≤ u < 1 and T A (1) = 1, we can apply the result of [10, p. 294, Th. V.1]:

[u n ] 1 1 -T A (u) ∼ n→∞ 1 T ′ A (1) = 1 -E[A]. (2)
All the definitions can be extended to the multivariate case.

3 The single-server queue

In this section, we present the simple example of a single-server queue with one flow of packets, as depicted in Figure 2. The results presented here are not new (we eventually rediscover the Pollaczek-Khinchine formula, and apply tools developed by [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF]), but our aim is to present the method that will be generalized later.

Queueing model

The queue is initially empty. At each time slot t ≥ 1, one packet (if any) is served and then A t packets arrive. The sequence (A t ) t≥1 is i.i.d. with generating function A and mean λ < 1. We let X t denote the number of packets in the queue at the end of time slot t. The system is driven by the equations

X 0 = 0 and X t+1 = (X t -1) + + A t+1 , ∀t ≥ 0, (3)
where (•) + = max(•, 0).

X t ∼ Φ A t ∼ A 1
Figure 2: A single-server queue crossed by a single flow. The server processes one packet during each time slot.

Generating function

We define the generating function of the queue state as

Φ(u, z) = n≥0 t≥0 P(X t = n)u n z t . (4)
For each t ∈ N, taking the coefficient of z t yields

z t Φ(u, z) = n≥0 P(X t = n)u n ,
which is the generating function of X t . We will show that

Lemma 1. Φ(u, z) = 1 + zA(u) (Φ(u, z) -Φ(0, z))u -1 + Φ(0, z) . (5)
Sketch of proof. Eq. ( 3) implies X 0 = 0 and for any t ≥ 1, with the notation a n = P(A = n),

P(X t = n -1) = a n-1 (P(X t-1 = 0) + P(X t-1 = 1)) + n m=2 a n-m P(X t-1 = m).
Multiplying this relation by u n z t , summing over n, t and dividing both sides by u leads to the equation of the lemma. The formula can also be directly derived using the Symbolic Method [START_REF] Flajolet | Analytic Combinatorics[END_REF].

Eq. ( 5) can be rewritten as

Φ(u, z) 1 -zA(u)u -1 = 1 -Φ(0, z)zA(u) u -1 -1 . (6)
Although Eq. ( 6) completely characterizes Φ(u, z), it is not straightforward to derive an explicit formula for Φ(u, z) from it, as we would need an expression for Φ(0, z). This expression will be obtained with the kernel method.

Kernel method

When the left-hand side of Eq. ( 6) cancels, so does the right-hand side. The kernel method [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF][START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF] consists in taking u = U (z) such that the second factor of the left-hand side cancels. Here, U (z) is implicitly defined by the equality U (z) = zA(U (z)), and we recognize from Eq. ( 1) the size distribution of a Galton-Watson tree where the offspring distribution has the generating function A. Therefore, we have U = T A .

Injecting T A (z) in Eq. ( 6) cancels its left-hand side, and its right-hand side can be rewritten as

Φ(0, z) = 1 1 -T A (z)
.

Going back to Eq. ( 6), we finally obtain

Φ(u, z) = 1 + 1 1-T A (z) zA(u) 1 -u -1 1 -zA(u)u -1 . ( 7
)
The kernel method has the following interpretation in terms of the queue sample paths. The generating function Φ(0, z) = t≥0 P(X t = 0)z t is associated with the probability of having an empty queue. Consider the duration between two consecutive instants when the queue is empty, which we call an inter-empty period. It was showed in [START_REF] Kendall | Some Problems in the Theory of Queues[END_REF] that we can build a Galton-Watson tree with offspring distribution A from an inter-empty period: each node represents a time slot; its children are the time slots when the packets arrived during this time slot are served. Having an empty queue at time t means that the realization between times 0 and t is made up of an arbitrary number of inter-empty periods. This corresponds exactly to

1 1-T A (z)
, where T A is as defined in Eq. (1).

Asymptotic performance

In this paragraph, our aim is to bound the probability that X t exceeds some value R in stationary regime. Note that, by monotony, this will also be an upper bound for the initially empty queue. We proceed in two steps. We first compute Π, the generating function of the stationary distribution of (X t ), and then we derive the asymptotic behavior of Π.

Computation of Π We know that, under the stability condition A ′ (1) = λ < 1, the distribution of X t converges to a stationary distribution π as t tends to +∞. The first step of our analysis consists in finding the generating function Π of this distribution π. Recall that, for each t ∈ N, the generating function of X t is Π t (u) = [z t ]Φ(u, z). By [10, p. 624], it suffices to study the limit of Π t (u) as t tends to +∞, when u is fixed. The obtained limit is exactly Π(u).

Let us fix u = u 0 . We see in Eq. ( 7) that Φ(u 0 , z) has two potential poles, 1 and u 0 /A(u 0 ). It can be checked that T A ( u 0 A(u 0 ) ) = u 0 , so that u 0 /A(u 0 ) is actually not a pole. In order to derive the asymptotic behavior of Π t (u 0 ) as t tends to +∞, we first compute a simpler equivalent of Φ(u 0 , z) in the neighborhood of its pole z = 1. After some rewriting, we obtain

Φ(u 0 , z) = u 0 u 0 -zA(u 0 ) + 1 1 -T A (z) zA(u 0 )(u 0 -1) u 0 -zA(u 0 ) .
As a consequence,

Φ(u 0 , z) ∼ z→1 u 0 u 0 -A(u 0 ) + 1 1 -T A (z) A(u 0 )(u 0 -1) u 0 -A(u 0 ) ,
and from Eq. ( 2), the terms are equivalent to

[z t ]Φ(u 0 , z) ∼ t→∞ (1 -λ) A(u 0 )(u 0 -1) u 0 -A(u 0 ) .
Therefore, the generating function of π is equal to the one given by the Pollaczek-Khinchine formula

Π(u) = (1 -λ) A(u)(u -1) u -A(u) .
Performance The second solution β of the equation u = A(u) is the convergence radius of the function Π (with β = +∞ in the degenerate case where A(u) is linear). The error bound, i.e. the probability that the buffer occupancy is at least R, is n≥R π(n). Its generating function is

E(u) = R≥0 n≥R π(n) u R = 1 -uΠ(u) 1 -u . ( 8
)
The asymptotic analysis of this generating function yields

Theorem 1. With X ∼ Π, P(X ≥ R) ∼ R→∞ (1 -λ) β A ′ (β) -1 β -R . ( 9
)
4 Extensions of the single-server queue

The analysis in the previous section shows that deriving an equation satisfied by the generating function from the system dynamics is the easy step; solving this equation is harder. We now consider a few simple extensions of the model of §3.1, where the kernel method allows us to perform the analysis and derive explicit formulas for the performance metrics.

Random service

We consider a first extension of the model of §3.1 where the service is random. Specifically, at each time slot t ≥ 1, the server processes one packet (if any) with some probability p > λ, and zero packet otherwise. The system is driven by the equations X 0 = 0 and X t+1 = (X t -S t ) + + A t+1 , ∀t ≥ 0, where (S t ) t∈N is a sequence of independent, Bernoulli distributed random variables with parameter p. The corresponding generating function is S(u) = 1 -p + pu.

The generating function Φ of the system state is as defined in (4). The equation satisfied by Φ is a rewriting of Eq. ( 6), where u -1 is replaced by S(u -1 ):

Φ(u, z)[1 -zA(u)S(u -1 )] = 1 -Φ(0, z)zA(u)[S(u -1 ) -1].
Applying the kernel method consists in choosing u = U (z) such that zA(U (z))S(U (z) -1 ) = 1. We can rewrite this as U (z) = G(zA(U (z))), where G is the generating series of the geometric distribution (defined on the set of positive integers) with parameter p:

G(s) = ps 1 -(1 -p)s .
In much the same way as in §3.3, we obtain

Φ(0, z) = 1 1 -zA(U (z)) = 1 1 -T A•G (z)
.

The second equality holds because Ū (z) = zA(U (z)) satisfies the equation Ū (z) = zA(G( Ū (z))), so that Ū is also the generating function of the size of a Galton-Watson tree with offspring distribution A • G. Finally, we obtain

Φ(u, z) = 1 + 1 1-T A•G (z) zA(u)(1 -S(u -1 )) 1 -zA(u)S(u -1
) .

The interpretation is similar to that of §3.3, except that the number of time slots dedicated to a given customer is now geometrically distributed with parameter p.

The stationary distribution has the generating function

Π(u) = 1 - λ p A(u) 1 -S(u -1 ) 1 -A(u)S(u -1
) .

Let γ be the largest solution of the equation A(u)S(u -1 ) = 1, or, equivalently, u = G(A(u)). Similarly to §3.4, we obtain Theorem 2. For X ∼ Π,

P(X ≥ R) ∼ R→∞ (1 -λ p )(A(γ) -1) (γ -1) A ′ (γ)S(γ -1 ) -A(γ)pγ -2 γ -R .
The kernel method is also applicable to the case where the server processes up to c packets at each time slot, for some integer c ≥ 1. Assume that the distribution of the number of served packets has generating function S(u) if the queue contains at least c packets, and S k (u) if the queue contains exactly k packets, for each 0 ≤ k < c. The equivalent of Eq. ( 6) is now [START_REF] Flajolet | Analytic Combinatorics[END_REF] We refer the reader to [10, p. 508] or [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF] for a detailed analysis of this equation, and provide here a short version. There are c independent functions (U k (z)) 0≤k<c , analytic at 0, that cancel the second term of the left hand-side, because S is a degree c polynomial. Thus, we obtain c equations for the c unknowns ∂ k (∂u) k Φ(0, z) for 0 ≤ k < c. Solving this system of equations and injecting the solution in Eq. ( 10) leads to the expression of Φ(u, z).

Φ(u, z)[1 -zA(u)S(u -1 )] = 1 - c-1 k=0 zA(u)(S(u -1 ) -S k (u -1 )) u k k! ∂ k (∂u) k Φ(0, z).

Several flows with priorities

Consider the system in Figure 3. As in §3.1, the queue is initially empty and the server processes one packet at each time slot. There are two flows of packets. Flow 1 has priority over flow 2, so that a packet from flow 1 is served whenever the queue contains at least one packet from this flow at the beginning of this time slot. At each time slot t ≥ 1, A t packets from flow 1 and B t packets from flow 2 arrive. The sequences ( and (B t ) t≥1 are independent and i.i.d. with generating function A and B and mean λ A and λ B , respectively, such that λ A + λ B < 1. We denote by X t and Y t the respective numbers of packets from flows 1 and 2 in the queue at the end of time slot t. The system is then driven by the equations X 0 = 0, Y 0 = 0 and

A t ) t≥1 X t , Y t ∼ Φ A t ∼ A B t ∼ B 1
X t+1 = (X t -1) + + A t+1 , Y t+1 = (Y t -1 {Xt=0} ) + + B t+1 , ∀t ≥ 0. ( 11
)
We define a generating function for the state (X t , Y t ) t≥0 , with three variables u, v, and z, respectively representing the numbers of packets from flows 1 and 2 and the time:

Φ(u, v, z) = n≥0 m≥0 t≥0 P(X t = n, Y t = m)u n v m z t .
The equation satisfied by Φ follows from Eq. ( 11):

Φ(u, v, z) = 1 + zA(u)B(v) (Φ(u, v, z) -Φ(0, v, z))u -1 + (Φ(0, v, z) -Φ(0, 0, z))v -1 + Φ(0, 0, z) .
Here (Φ(u, v, z) -Φ(0, v, z))u -1 represents the service of a packet from flow 1 (if any) and (Φ(0, v, z) -Φ(0, 0, z))v -1 the service of a packet from flow 2 (if any, and when there is no packet from flow 1). This equation can be rewritten as

Φ(u, v, z) 1 -zA(u)B(v)u -1 = 1 -zA(u)B(v)× Φ(0, v, z) u -1 -v -1 + Φ(0, 0, z) v -1 -1 . (12)
We could find an expression for Φ(0, 0, z) by applying twice the kernel method on this equation, but we prefer a more intuitive approach. The generating function Φ(0, 0, z) is associated with the probability that the queue is empty. This probability only depends on the global arrival process (A t + B t ) t≥1 , regardless of the division of packets into flows. Therefore, we know from §3.3 that

Φ(0, 0, z) = 1 1 -T AB (z)
,

where T AB is the generating function of the size of a Galton-Watson tree with offspring distribution AB.

We apply the kernel method to derive the expression of Φ(0, v, z). Let us take u = U (v, z) such that U (v, z) = zB(v)A(U (v, z)), in order to cancel the left-hand side of Eq. [START_REF] Kendall | Some Problems in the Theory of Queues[END_REF]. We obtain U (v, z) = T A (zB(v)), which is again strongly related to a Galton-Watson tree. The interpretation is similar to that of §3.3, except that U (v, z) is now the generating function of the number of time slots passed and flow-2 packets arrived during an inter-empty period of flow 1. The priority of flow 1 ensures that no packet from flow 2 is served in the meantime. After simplifications, we get

Φ(0, v, z) = v + 1 1-T AB (z) T A (zB(v))(v -1) v -T A (zB(v)) ,
and the expression for Φ(u, v, z) immediately follows. It is not difficult to see that this method can be generalized to queues with more than two flows with a total order on priority levels. Suppose that we focus on the number of packets from flow 2 in the stationary state. We are then interested in the generating function Φ (1, v, z). The same approach as in §3.4 can be used to obtain the following result.

Theorem 3. The limit distribution of the number of packets of flow 2 is given by the generating series

Π(v) = (1 -λ A -λ B ) B(v)(1 -v)(T A (B(v)) -1) (1 -B(v))(v -T A (B(v)) .
Let δ be the largest solution of the equation v = T A (B(v)) (as T A and B are convex, there are exactly 2 solutions, and the smallest is 1). Similarly to §3.4,

Theorem 4. For Y ∼ Π, P(Y ≥ R) ∼ R→∞ (1 -λ A -λ B )B(δ)(δ -1) (1 -B(δ))(1 -(T A • B) ′ (δ)) δ -R . ( 13
)

Several queues

We can also use generating functions to describe the dynamics of networks of queues.

As an example, consider the network of Figure 4, consisting of two single-server queues. At each time slot t ≥ 1, A t packets arrive at queue 1 and B t packets arrive at queue 2.

As before, the sequences (A t ) t≥1 and (B t ) t≥1 are independent and i.i.d. with generating function A and B and mean λ A and λ B , respectively, such that λ A + λ B < 1. Additionally, the packets served at queue 1 are subsequently forwarded to queue 2 for service. We let X t and Y t denote the numbers of packets at queues 1 and 2, respectively, at time t.

The dynamics of the system, which is initially empty, are driven by the equations X 0 = 0, Y 0 = 0 and We define a generating function for the state (X t , Y t ) t≥0 , with three variables u, v, and z, respectively representing the numbers of packets at queues 1 and 2 and the time:

X t+1 = (X t -1) + + A t+1 , Y t+1 = (Y t -1) + + B t+1 + 1 {Xt>0} , ∀t ≥ 0. ( 14
) X t ∼ Φ 1 Y t ∼ Φ 1 A t ∼ A B t ∼ B
Φ(u, v, z) = n≥0 m≥0 t≥0 P(X t = n, Y t = m)u n v m z t .
This function Φ satisfies the following equation:

Φ(u, v, z)[1 -zA(u)B(v)u -1 ] = 1-zA(u)B(v) [Φ(u, 0, z)(u -1 -vu -1 ) + Φ(0, v, z)(u -1 -v -1 ) + Φ(0, 0, z)(vu -1 + v -1 -u -1 -1)]. (15)
The kernel method cannot be applied directly. Indeed, we need to compute three generating functions (Φ(u, 0, z), Φ(0, v, z), and Φ(0, 0, z)), while we can only apply the kernel method (at most) twice. It is, however, possible to find an additional relation between these functions:

Φ(u, 0, z) = 1 + zA(u)B(0) Φ(0, 0, z) + [v 1 ]Φ(0, v, z) .
This relation can be obtained in two different ways. We can derive Eq. ( 15) according to v at v = 0. Alternatively, we can go back to the system dynamics: queue 2 is empty at the end of some time slot t ≥ 1 if it does not receive any external arrival during this time slot and, at the end of time slot t -1, queue 1 was empty and queue 2 contained at most one packet.

This equation gives the relation

Φ(u, 0, z) = 1 + A(u) A(0) [Φ(0, 0, z) -1].
We are now in a position to apply the kernel method, by defining first

U (v, z) = zA(U (v, z))B(v) and then V (z) = zA(V (z))B(V (z)).
The generating function Π of the stationary distribution of the number of packets in queue 2 can be computed similarly to the previous cases. For simplicity, we only give its asymptotic behavior: Buffer occupancy R δ -R Eq. ( 16) Simulation (c) Two servers Although we obtain the exact asymptotic in those two cases, it seems that these are lower bounds of the error. Indeed, we only computed the first term. But, once again, as we were able to compute an exact formula for the error bounds, more terms are derivable using Taylor expansions.

Theorem 5. With δ previously defined, for Y ∼ Π, P(Y ≥ R) ∼ R→∞ (1 -λ A -λ B )δ(δ -1) (1 -B(δ))(1 -(T A • B) ′ (δ)) δ -R . ( 16

Conclusion

In this paper, we have demonstrated on simple examples that methods from analytic combinatorics can be successfully applied to the analysis of queueing systems. We have focused on computing backlog bounds, but we believe delay bounds can be derived by using the same techniques as in [START_REF] Poloczek | Scheduling analysis with martingales[END_REF], for the FIFO, EDF (earliest-deadline-first) and priorities policies. Moreover, combining the computations described in §4 would allow other service policies to enter our framework, in particular some discrete version of GPS (generalized processor sharing). Following the approach of [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF], we could also consider a continuous-time extension of our work based on Laplace transforms. The greatest challenge is to cope with networks of queues. A simple example with two queues has been analyzed. The same analysis can be extended for more than two queues, but this analysis is still partial since the packets are aggregated at each queue. Further investigation needs to be done, in particular in view of the techniques presented in [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF].

Figure 3 :

 3 Figure 3: A single-server queue crossed by two flows.

Figure 4 :

 4 Figure 4: A network of two queues crossed by two flows.

Figure 5 :

 5 Figure 5: Numerical evaluation of the kernel method. Parameters A = D 2/30,6 and B = D 2/5,1 .

Non-independent arrivals

The analysis can be extended to networks with more generic arrival processes. We take the network of Figure 3 as an example.

• The arrivals of flows 1 and 2 may be dependent. The global arrival process is then described by a generating function A(u, v) that cannot be written as a product A(u)B(v).

• Within each flow, the numbers of arrivals at different time slots may not be i.i.d. anymore. Instead, they may be described by a modulated process (which includes modulated Markov On-Off processes). The modulation is described by a finite Markov chain. The system dynamics are then described by a system of equations on generating functions (one per state of the Markov chain).

Numerical evaluation

We tested our formulas against simulations in three different scenarios: the single-server case of §3, the multiflow single-server case of §4.2 and the tandem network of two singleserver queues of §4.3. The service is deterministic. Performing simulations consists in computing the stationary distribution of the truncated processes (the number of packets in each queue never exceeds 200) whose dynamics are described by Eqs. ( 3), [START_REF] Jiang | Stochastic Network Calculus[END_REF] or [START_REF] Poloczek | Scheduling analysis with martingales[END_REF]. The approximation of the stationary distribution is obtained by iteratively computing the distribution after t steps for a large enough t (the stopping criterion is when the distance in total variation between the t-th and the t + 1-th distribution is less than 10 -12 ).

Each arrival process has a bimodal distribution with generating function D p,M (u) = (1 -p) + pu M , for some p and M : at each time slot, either M packets arrive, which occurs with probability p, or no packet arrives. With this distribution, the arrival rate is D ′ p,M (1) = M p. In the numerical results of Figure 5, we take A = D 2/30,6 and B = D 2/5,1 . This choice of the functions A and B is arbitrary (other distributions lead to similar observations) Fig. 5a illustrates the case of a single queue. The curve β 1-R is the one that would be obtained using Doob's inequality from [START_REF] Poloczek | Scheduling analysis with martingales[END_REF] (we do not use β -R because we consider the size of the queue before service and not after as in [START_REF] Poloczek | Scheduling analysis with martingales[END_REF]). The simulation confirms that we obtain the exact asymptotic behavior, and shows that we improve Doob's inequality by a factor 1.5. We remark that the simulation curve has some irregularities for small values of R. This is explained by the arrival of packets in batches of 6. It is possible to compute the exact error bound from Eq. ( 8) by deriving the first terms explicitly: 5b illustrates the case of a single-server queue with two flows. We focus on the buffer occupancy of flow 2. Since flow 1 has priority, its buffer occupancy is still given by Fig. 5a. Again, the simulation validates our theoretical results. Up to our knowledge, there is no formula similar to Doob's inequality, so the curve δ -R only mimics a Doob-like inequality.